weipa/src/DataVar.cpp


/*******************************************************
*
* Copyright (c) 2003-2010 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/

#include <weipa/DataVar.h>
#include <weipa/DomainChunk.h>
#include <weipa/ElementData.h>
#include <weipa/NodeData.h>
#ifndef VISIT_PLUGIN
#include <escript/Data.h>
#endif

#if USE_NETCDF
#include <netcdfcpp.h>
#endif

#if USE_SILO
#include <silo.h>
#endif

#include <numeric> // for accumulate
#include <iostream> // for cerr
#include <stdio.h>

using namespace std;

namespace weipa {
    
//
// Constructor
//
DataVar::DataVar(const string& name) :
    initialized(false), varName(name),
    numSamples(0), rank(0), ptsPerSample(0)
{
}

//
// Copy constructor
//
DataVar::DataVar(const DataVar& d) :
    initialized(d.initialized), domain(d.domain),
    varName(d.varName), numSamples(d.numSamples),
    rank(d.rank), ptsPerSample(d.ptsPerSample), funcSpace(d.funcSpace),
    centering(d.centering), shape(d.shape), sampleID(d.sampleID)
{
    if (numSamples > 0) {
        CoordArray::const_iterator it;
        for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
            float* c = new float[numSamples];
            copy(*it, (*it)+numSamples, c);
            dataArray.push_back(c);
        }
    }
}

//
// Destructor
//
DataVar::~DataVar()
{
    cleanup();
}

//
//
//
void DataVar::cleanup()
{
    CoordArray::iterator it;
    for (it = dataArray.begin(); it != dataArray.end(); it++)
        delete[] *it;
    dataArray.clear();
    shape.clear();
    sampleID.clear();
    numSamples = 0;
    initialized = false;
}

//
//
//
bool DataVar::initFromEscript(escript::Data& escriptData, const_DomainChunk_ptr dom)
{
#ifndef VISIT_PLUGIN
    cleanup();

    if (!escriptData.isConstant() && !escriptData.actsExpanded()) {
        cerr << "WARNING: Weipa only supports constant & expanded data, "
            << "not initializing " << varName << endl;
        return false;
    }

    domain = dom;
    rank = escriptData.getDataPointRank();
    ptsPerSample = escriptData.getNumDataPointsPerSample();
    shape = escriptData.getDataPointShape();
    funcSpace = escriptData.getFunctionSpace().getTypeCode();
    numSamples = escriptData.getNumSamples();
    centering = domain->getCenteringForFunctionSpace(funcSpace);

#ifdef _DEBUG
    cout << varName << ":\t" << numSamples << " samples,  "
        << ptsPerSample << " pts/s,  rank: " << rank << endl;
#endif

    NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
    if (nodes == NULL)
        return false;

    meshName = nodes->getName();
    siloMeshName = nodes->getFullSiloName();
    initialized = true;

    // no samples? Nothing more to do.
    if (numSamples == 0)
        return true;

    const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
    sampleID.insert(sampleID.end(), numSamples, 0);
    copy(iPtr, iPtr+numSamples, sampleID.begin());

    size_t dimSize = 1;
    if (rank > 0)
        dimSize *= shape[0];
    if (rank > 1)
        dimSize *= shape[1];
    if (rank > 2) {
        cerr << "WARNING: Rank " << rank << " data is not supported!\n";
        initialized = false;
    }

    if (initialized) {
        size_t dataSize = dimSize * ptsPerSample;
        float* tempData = new float[dataSize*numSamples];
        float* destPtr = tempData;
        if (escriptData.isConstant()) {
            const escript::DataAbstract::ValueType::value_type* values =
                escriptData.getSampleDataRO(0);
            for (int pointNo=0; pointNo<numSamples*ptsPerSample; pointNo++) {
                copy(values, values+dimSize, destPtr);
                destPtr += dimSize;
            }
        } else {
            for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
                const escript::DataAbstract::ValueType::value_type* values =
                    escriptData.getSampleDataRO(sampleNo);
                copy(values, values+dataSize, destPtr);
                destPtr += dataSize;
            }
        }

        const float* srcPtr = tempData;
        for (int i=0; i < dimSize; i++, srcPtr++) {
            float* c = averageData(srcPtr, dimSize);
            dataArray.push_back(c);
        }
        delete[] tempData;

        initialized = reorderSamples();
    }

    return initialized;

#else // VISIT_PLUGIN
    return false;
#endif
}

//
// Initialise with domain data
//
bool DataVar::initFromMeshData(const_DomainChunk_ptr dom, const IntVec& data,
        int fsCode, Centering c, NodeData_ptr nodes, const IntVec& id)
{
    cleanup();
    
    domain = dom;
    rank = 0;
    ptsPerSample = 1;
    centering = c;
    sampleID = id;
    meshName = nodes->getName();
    siloMeshName = nodes->getFullSiloName();
    numSamples = data.size();

    if (numSamples > 0) {
        float* c = new float[numSamples];
        dataArray.push_back(c);
        IntVec::const_iterator it;
        for (it=data.begin(); it != data.end(); it++)
            *c++ = static_cast<float>(*it);
    }
    initialized = true;

    return initialized;
}

//
// Reads variable data from dump file
//
bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
{
    cleanup();
    
#if USE_NETCDF
    NcError ncerr(NcError::silent_nonfatal);    
    NcFile* input = new NcFile(filename.c_str());
    if (!input->is_valid()) {
        cerr << "Could not open input file " << filename << "." << endl;
        delete input;
        return false;
    }

    NcDim* dim;
    NcAtt* att;

    att = input->get_att("type_id");
    int typeID = att->as_int(0);
    if (typeID != 2) {
        cerr << "WARNING: Only expanded data supported!" << endl;
        delete input;
        return false;
    }

    att = input->get_att("rank");
    rank = att->as_int(0);

    dim = input->get_dim("num_data_points_per_sample");
    ptsPerSample = dim->size();

    att = input->get_att("function_space_type");
    funcSpace = att->as_int(0);

    centering = dom->getCenteringForFunctionSpace(funcSpace);

    dim = input->get_dim("num_samples");
    numSamples = dim->size();

#ifdef _DEBUG
    cout << varName << ":\t" << numSamples << " samples,  "
        << ptsPerSample << " pts/s,  rank: " << rank << endl;
#endif

    domain = dom;
    NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
    if (nodes == NULL) {
        delete input;
        return false;
    }

    meshName = nodes->getName();
    siloMeshName = nodes->getFullSiloName();
    initialized = true;

    size_t dimSize = 1;
    vector<long> counts;

    if (rank > 0) {
        dim = input->get_dim("d0");
        int d = dim->size();
        shape.push_back(d);
        counts.push_back(d);
        dimSize *= d;
    }
    if (rank > 1) {
        dim = input->get_dim("d1");
        int d = dim->size();
        shape.push_back(d);
        counts.push_back(d);
        dimSize *= d;
    }
    if (rank > 2) {
        cerr << "WARNING: Rank " << rank << " data is not supported!\n";
        initialized = false;
    }
 
    if (initialized && numSamples > 0) {
        sampleID.insert(sampleID.end(), numSamples, 0);
        NcVar* var = input->get_var("id");
        var->get(&sampleID[0], numSamples);

        size_t dataSize = dimSize*numSamples*ptsPerSample;
        counts.push_back(ptsPerSample);
        counts.push_back(numSamples);
        float* tempData = new float[dataSize];
        var = input->get_var("data");
        var->get(tempData, &counts[0]);

        const float* srcPtr = tempData;
        for (int i=0; i < dimSize; i++, srcPtr++) {
            float* c = averageData(srcPtr, dimSize);
            dataArray.push_back(c);
        }
        delete[] tempData;

        initialized = reorderSamples();
    }

    delete input;
#endif // USE_NETCDF

    return initialized;
}

//
// Returns true if the data values are nodal, false if they are zonal.
//
bool DataVar::isNodeCentered() const
{
    return (centering == NODE_CENTERED);
}

//
// Returns a subset of the src array according to stride parameter.
// If samples consist of multiple values they are averaged beforehand.
// Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
// (z0,z1,...)
//
float* DataVar::averageData(const float* src, size_t stride)
{
    float* res;

    if (ptsPerSample == 1) {
        res = new float[numSamples];
        float* dest = res;
        for (int i=0; i<numSamples; i++, src+=stride)
            *dest++ = *src;
    } else {
        ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
        int cellFactor = cells->getElementFactor();
        res = new float[cellFactor * numSamples];
        float* dest = res;
        QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
        if (!qmi.mask.empty()) {
            const float* tmpSrc = src;
            for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
                for (int l=0; l<cellFactor; l++) {
                    double tmpVal = 0.0;
                    for (int j=0; j<ptsPerSample; j++) {
                        if (qmi.mask[l][j] != 0) {
                            tmpVal += *(tmpSrc+stride*j);
                        }
                    }
                    *dest++ = (float)(tmpVal / qmi.factor[l]);
                }
            }
        } else {
            for (int i=0; i<numSamples; i++) {
                double tmpVal = 0.0;
                for (int j=0; j<ptsPerSample; j++, src+=stride) {
                    tmpVal += *src;
                }
                tmpVal /= ptsPerSample;
                for (int l=0; l<cellFactor; l++) {
                    *dest++ = static_cast<float>(tmpVal);
                }
            }
        }
    }
    return res;
}

//
// Filters and reorders the raw sample values according to the node/element
// IDs. This is used to have data arrays ordered according to the underlying
// mesh (i.e. DataID[i]==MeshNodeID[i])
//
bool DataVar::reorderSamples()
{
    if (numSamples == 0)
        return true;

    const IntVec* requiredIDs = NULL;
    int requiredNumSamples = 0;
    int cellFactor = 1;

    if (centering == NODE_CENTERED) {
        NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
        requiredIDs = &nodes->getNodeIDs();
        requiredNumSamples = nodes->getNumNodes();
    } else {
        ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
        if (cells == NULL)
            return false;

        requiredIDs = &cells->getIDs();
        requiredNumSamples = cells->getNumElements();
        cellFactor = cells->getElementFactor();
        if (cellFactor > 1) {
            numSamples *= cellFactor;
            // update sample IDs
            IntVec newSampleID(numSamples);
            IntVec::const_iterator idIt = sampleID.begin();
            IntVec::iterator newIDit = newSampleID.begin();
            for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
                fill(newIDit, newIDit+cellFactor, *idIt);
            }
            sampleID.swap(newSampleID);
        }
    }

    if (requiredNumSamples > numSamples) {
        cerr << "ERROR: " << varName << " has " << numSamples
            << " instead of " << requiredNumSamples << " samples!" << endl;
        return false;
    }

    IndexMap sampleID2idx = buildIndexMap();
    numSamples = requiredNumSamples;

    // now filter the data
    for (size_t i=0; i < dataArray.size(); i++) {
        float* c = new float[numSamples];
        const float* src = dataArray[i];
        IntVec::const_iterator idIt = requiredIDs->begin();
        size_t destIdx = 0;
        for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
            size_t srcIdx = sampleID2idx.find(*idIt)->second;
            copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
        }
        delete[] dataArray[i];
        dataArray[i] = c;
    }

    // sample IDs now = mesh node/element IDs
    sampleID = *requiredIDs;

    return true;
}

//
//
//
int DataVar::getNumberOfComponents() const
{
    return (rank == 0 ? 1 : accumulate(shape.begin(), shape.end(), 0));
}

//
//
//
float* DataVar::getDataFlat() const
{
    int totalSize = numSamples * getNumberOfComponents();
    float* res = new float[totalSize];
    if (rank == 0) {
        copy(dataArray[0], dataArray[0]+numSamples, res);
    } else if (rank == 1) {
        float *dest = res;
        for (size_t c=0; c<numSamples; c++) {
            for (size_t i=0; i<shape[0]; i++) {
                *dest++ = dataArray[i][c];
            }
        }
    } else if (rank == 2) {
        float *dest = res;
        for (size_t c=0; c<numSamples; c++) {
            for (int i=0; i<shape[1]; i++) {
                for (int j=0; j<shape[0]; j++) {
                    *dest++ = dataArray[i*shape[0]+j][c];
                }
            }
        }
    }

    return res;
}

//
//
//
void DataVar::sampleToStream(ostream& os, int index)
{
    // index is -1 for dummy samples, i.e. if writing the full mesh but
    // only a reduced number of samples is required
    if (rank == 0) {
        if (index < 0) {
            os << 0.;
        } else {
            os << dataArray[0][index];
        }
    } else if (rank == 1) {
        if (index < 0) {
            os << 0. << " " << 0.  << " " << 0.;
        } else if (shape[0] < 3) {
            os << dataArray[0][index] << " " << dataArray[1][index]
                << " " << 0.;
        } else {
            os << dataArray[0][index] << " " << dataArray[1][index]
                << " " << dataArray[2][index];
        }
    } else if (rank == 2) {
        if (index < 0) {
            os << 0. << " " << 0. << " " << 0. << " ";
            os << 0. << " " << 0. << " " << 0. << " ";
            os << 0. << " " << 0. << " " << 0.;
        } else if (shape[1] < 3) {
            os << dataArray[0][index] << " " << dataArray[1][index]
                << " " << 0. << " ";
            os << dataArray[2][index] << " " << dataArray[3][index]
                << " " << 0. << " ";
            os << 0. << " " << 0. << " " << 0.;
        } else {
            os << dataArray[0][index] << " " << dataArray[1][index]
                << " " << dataArray[2][index] << " ";
            os << dataArray[3][index] << " " << dataArray[4][index]
                << " " << dataArray[5][index] << " ";
            os << dataArray[6][index] << " " << dataArray[7][index]
                << " " << dataArray[8][index];
        }
    }
    os << endl;
}

//
//
//
void DataVar::writeToVTK(ostream& os, int ownIndex)
{
    if (numSamples == 0)
        return;

    if (isNodeCentered()) {
        // data was reordered in reorderSamples() but for VTK we write the
        // original node mesh and thus need the original ordering.
        // Note, that this also means we may not have samples for all nodes
        // in which case we set idx to -1 and write a dummy sample
        const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
        const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
        const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
        int firstId = nodeDist[ownIndex];
        int lastId = nodeDist[ownIndex+1];
        IndexMap sampleID2idx = buildIndexMap();
        for (int i=0; i<nodeGNI.size(); i++) {
            if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
                IndexMap::const_iterator it = sampleID2idx.find(requiredIDs[i]);
                int idx = (it==sampleID2idx.end() ? -1 : (int)it->second);
                sampleToStream(os, idx);
            }
        }
    } else {
        // cell data: ghost cells have been removed so do not write ghost
        // samples (which are the last elements in the arrays)
        int toWrite = domain->getElementsByName(meshName)->getNumElements();
        for (int i=0; i<toWrite; i++) {
            sampleToStream(os, i);
        }
    }
}

///////////////////////////////
// SILO related methods follow
///////////////////////////////

//
// If the data is tensor data then the components of the tensor are stored
// separately in the Silo file. This method then returns a string that
// contains the proper Silo expression to put the tensor together again.
// For non-tensor data this method returns an empty string.
//
string DataVar::getTensorDef() const
{
    if (rank < 2 || !initialized)
        return string();
    
    /// Format string for Silo 2x2 tensor
    const string tensor2DefFmt =
        "{{ <%sa_00>, <%sa_01> },"
        " { <%sa_10>, <%sa_11> }}";

    /// Format string for Silo 3x3 tensor
    const string tensor3DefFmt =
        "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
        " { <%sa_10>, <%sa_11>, <%sa_12> },"
        " { <%sa_20>, <%sa_21>, <%sa_22> }}";

    string tensorDef;
    string tensorDir = varName+string("_comps/");
    if (shape[1] == 3) {
        char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
        sprintf(tDef, tensor3DefFmt.c_str(),
                tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
                tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
                tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
        tensorDef = tDef;
        delete[] tDef;
    } else {
        char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
        sprintf(tDef, tensor2DefFmt.c_str(),
                tensorDir.c_str(), tensorDir.c_str(),
                tensorDir.c_str(), tensorDir.c_str(),
                tensorDir.c_str(), tensorDir.c_str());
        tensorDef = tDef;
        delete[] tDef;
    }
    return tensorDef;
}

//
// Writes the data to given Silo file under the virtual path provided.
// The corresponding mesh must have been written already.
//
bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
                          const string& units)
{
#if USE_SILO
    if (!initialized)
        return false;

    if (numSamples == 0)
        return true;

    int ret;

    if (siloPath != "") {
        ret = DBSetDir(dbfile, siloPath.c_str());
        if (ret != 0)
            return false;
    }
 
    char* siloMesh = const_cast<char*>(siloMeshName.c_str());
    int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
    DBoptlist* optList = DBMakeOptlist(2);
    if (units.length()>0) {
        DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
    }

    if (rank == 0) {
        ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
                numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
    }
    else if (rank == 1) {
        const string comps[3] = {
            varName+string("_x"), varName+string("_y"), varName+string("_z")
        };
        const char* varnames[3] = {
            comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
        };

        ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
                (char**)varnames, &dataArray[0], numSamples, NULL,
                0, DB_FLOAT, dcenter, optList);
    }
    else {
        string tensorDir = varName+string("_comps/");
        ret = DBMkdir(dbfile, tensorDir.c_str());
        if (ret == 0) {
            int one = 1;
            DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);

            for (int i=0; i<shape[1]; i++) {
                for (int j=0; j<shape[0]; j++) {
                    ostringstream varname;
                    varname << tensorDir << "a_" << i << j;
                    ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
                            dataArray[i*shape[0]+j], numSamples,
                            NULL, 0, DB_FLOAT, dcenter, optList);
                    if (ret != 0) break;
                }
                if (ret != 0) break;
            }
        } // ret==0
    } // rank

    DBFreeOptlist(optList);
    DBSetDir(dbfile, "/");
    return (ret == 0);

#else // !USE_SILO
    return false;
#endif
}

} // namespace weipa

1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2010 by University of Queensland
5	* Earth Systems Science Computational Center (ESSCC)
6	* http://www.uq.edu.au/esscc
7	*
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14	#include <weipa/DataVar.h>
15	#include <weipa/DomainChunk.h>
16	#include <weipa/ElementData.h>
17	#include <weipa/NodeData.h>
18	#ifndef VISIT_PLUGIN
19	#include <escript/Data.h>
20	#endif
21
22	#if USE_NETCDF
23	#include <netcdfcpp.h>
24	#endif
25
26	#if USE_SILO
27	#include <silo.h>
28	#endif
29
30	#include <numeric> // for accumulate
31	#include <iostream> // for cerr
32	#include <stdio.h>
33
34	using namespace std;
35
36	namespace weipa {
37
38	//
39	// Constructor
40	//
41	DataVar::DataVar(const string& name) :
42	initialized(false), varName(name),
43	numSamples(0), rank(0), ptsPerSample(0)
44	{
45	}
46
47	//
48	// Copy constructor
49	//
50	DataVar::DataVar(const DataVar& d) :
51	initialized(d.initialized), domain(d.domain),
52	varName(d.varName), numSamples(d.numSamples),
53	rank(d.rank), ptsPerSample(d.ptsPerSample), funcSpace(d.funcSpace),
54	centering(d.centering), shape(d.shape), sampleID(d.sampleID)
55	{
56	if (numSamples > 0) {
57	CoordArray::const_iterator it;
58	for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
59	float* c = new float[numSamples];
60	copy(it, (it)+numSamples, c);
61	dataArray.push_back(c);
62	}
63	}
64	}
65
66	//
67	// Destructor
68	//
69	DataVar::~DataVar()
70	{
71	cleanup();
72	}
73
74	//
75	//
76	//
77	void DataVar::cleanup()
78	{
79	CoordArray::iterator it;
80	for (it = dataArray.begin(); it != dataArray.end(); it++)
81	delete[] *it;
82	dataArray.clear();
83	shape.clear();
84	sampleID.clear();
85	numSamples = 0;
86	initialized = false;
87	}
88
89	//
90	//
91	//
92	bool DataVar::initFromEscript(escript::Data& escriptData, const_DomainChunk_ptr dom)
93	{
94	#ifndef VISIT_PLUGIN
95	cleanup();
96
97	if (!escriptData.isConstant() && !escriptData.actsExpanded()) {
98	cerr << "WARNING: Weipa only supports constant & expanded data, "
99	<< "not initializing " << varName << endl;
100	return false;
101	}
102
103	domain = dom;
104	rank = escriptData.getDataPointRank();
105	ptsPerSample = escriptData.getNumDataPointsPerSample();
106	shape = escriptData.getDataPointShape();
107	funcSpace = escriptData.getFunctionSpace().getTypeCode();
108	numSamples = escriptData.getNumSamples();
109	centering = domain->getCenteringForFunctionSpace(funcSpace);
110
111	#ifdef _DEBUG
112	cout << varName << ":\t" << numSamples << " samples, "
113	<< ptsPerSample << " pts/s, rank: " << rank << endl;
114	#endif
115
116	NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
117	if (nodes == NULL)
118	return false;
119
120	meshName = nodes->getName();
121	siloMeshName = nodes->getFullSiloName();
122	initialized = true;
123
124	// no samples? Nothing more to do.
125	if (numSamples == 0)
126	return true;
127
128	const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
129	sampleID.insert(sampleID.end(), numSamples, 0);
130	copy(iPtr, iPtr+numSamples, sampleID.begin());
131
132	size_t dimSize = 1;
133	if (rank > 0)
134	dimSize *= shape[0];
135	if (rank > 1)
136	dimSize *= shape[1];
137	if (rank > 2) {
138	cerr << "WARNING: Rank " << rank << " data is not supported!\n";
139	initialized = false;
140	}
141
142	if (initialized) {
143	size_t dataSize = dimSize * ptsPerSample;
144	float* tempData = new float[dataSize*numSamples];
145	float* destPtr = tempData;
146	if (escriptData.isConstant()) {
147	const escript::DataAbstract::ValueType::value_type* values =
148	escriptData.getSampleDataRO(0);
149	for (int pointNo=0; pointNo<numSamples*ptsPerSample; pointNo++) {
150	copy(values, values+dimSize, destPtr);
151	destPtr += dimSize;
152	}
153	} else {
154	for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
155	const escript::DataAbstract::ValueType::value_type* values =
156	escriptData.getSampleDataRO(sampleNo);
157	copy(values, values+dataSize, destPtr);
158	destPtr += dataSize;
159	}
160	}
161
162	const float* srcPtr = tempData;
163	for (int i=0; i < dimSize; i++, srcPtr++) {
164	float* c = averageData(srcPtr, dimSize);
165	dataArray.push_back(c);
166	}
167	delete[] tempData;
168
169	initialized = reorderSamples();
170	}
171
172	return initialized;
173
174	#else // VISIT_PLUGIN
175	return false;
176	#endif
177	}
178
179	//
180	// Initialise with domain data
181	//
182	bool DataVar::initFromMeshData(const_DomainChunk_ptr dom, const IntVec& data,
183	int fsCode, Centering c, NodeData_ptr nodes, const IntVec& id)
184	{
185	cleanup();
186
187	domain = dom;
188	rank = 0;
189	ptsPerSample = 1;
190	centering = c;
191	sampleID = id;
192	meshName = nodes->getName();
193	siloMeshName = nodes->getFullSiloName();
194	numSamples = data.size();
195
196	if (numSamples > 0) {
197	float* c = new float[numSamples];
198	dataArray.push_back(c);
199	IntVec::const_iterator it;
200	for (it=data.begin(); it != data.end(); it++)
201	c++ = static_cast<float>(it);
202	}
203	initialized = true;
204
205	return initialized;
206	}
207
208	//
209	// Reads variable data from dump file
210	//
211	bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
212	{
213	cleanup();
214
215	#if USE_NETCDF
216	NcError ncerr(NcError::silent_nonfatal);
217	NcFile* input = new NcFile(filename.c_str());
218	if (!input->is_valid()) {
219	cerr << "Could not open input file " << filename << "." << endl;
220	delete input;
221	return false;
222	}
223
224	NcDim* dim;
225	NcAtt* att;
226
227	att = input->get_att("type_id");
228	int typeID = att->as_int(0);
229	if (typeID != 2) {
230	cerr << "WARNING: Only expanded data supported!" << endl;
231	delete input;
232	return false;
233	}
234
235	att = input->get_att("rank");
236	rank = att->as_int(0);
237
238	dim = input->get_dim("num_data_points_per_sample");
239	ptsPerSample = dim->size();
240
241	att = input->get_att("function_space_type");
242	funcSpace = att->as_int(0);
243
244	centering = dom->getCenteringForFunctionSpace(funcSpace);
245
246	dim = input->get_dim("num_samples");
247	numSamples = dim->size();
248
249	#ifdef _DEBUG
250	cout << varName << ":\t" << numSamples << " samples, "
251	<< ptsPerSample << " pts/s, rank: " << rank << endl;
252	#endif
253
254	domain = dom;
255	NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
256	if (nodes == NULL) {
257	delete input;
258	return false;
259	}
260
261	meshName = nodes->getName();
262	siloMeshName = nodes->getFullSiloName();
263	initialized = true;
264
265	size_t dimSize = 1;
266	vector<long> counts;
267
268	if (rank > 0) {
269	dim = input->get_dim("d0");
270	int d = dim->size();
271	shape.push_back(d);
272	counts.push_back(d);
273	dimSize *= d;
274	}
275	if (rank > 1) {
276	dim = input->get_dim("d1");
277	int d = dim->size();
278	shape.push_back(d);
279	counts.push_back(d);
280	dimSize *= d;
281	}
282	if (rank > 2) {
283	cerr << "WARNING: Rank " << rank << " data is not supported!\n";
284	initialized = false;
285	}
286
287	if (initialized && numSamples > 0) {
288	sampleID.insert(sampleID.end(), numSamples, 0);
289	NcVar* var = input->get_var("id");
290	var->get(&sampleID[0], numSamples);
291
292	size_t dataSize = dimSizenumSamplesptsPerSample;
293	counts.push_back(ptsPerSample);
294	counts.push_back(numSamples);
295	float* tempData = new float[dataSize];
296	var = input->get_var("data");
297	var->get(tempData, &counts[0]);
298
299	const float* srcPtr = tempData;
300	for (int i=0; i < dimSize; i++, srcPtr++) {
301	float* c = averageData(srcPtr, dimSize);
302	dataArray.push_back(c);
303	}
304	delete[] tempData;
305
306	initialized = reorderSamples();
307	}
308
309	delete input;
310	#endif // USE_NETCDF
311
312	return initialized;
313	}
314
315	//
316	// Returns true if the data values are nodal, false if they are zonal.
317	//
318	bool DataVar::isNodeCentered() const
319	{
320	return (centering == NODE_CENTERED);
321	}
322
323	//
324	// Returns a subset of the src array according to stride parameter.
325	// If samples consist of multiple values they are averaged beforehand.
326	// Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
327	// (z0,z1,...)
328	//
329	float* DataVar::averageData(const float* src, size_t stride)
330	{
331	float* res;
332
333	if (ptsPerSample == 1) {
334	res = new float[numSamples];
335	float* dest = res;
336	for (int i=0; i<numSamples; i++, src+=stride)
337	dest++ = src;
338	} else {
339	ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
340	int cellFactor = cells->getElementFactor();
341	res = new float[cellFactor * numSamples];
342	float* dest = res;
343	QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
344	if (!qmi.mask.empty()) {
345	const float* tmpSrc = src;
346	for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
347	for (int l=0; l<cellFactor; l++) {
348	double tmpVal = 0.0;
349	for (int j=0; j<ptsPerSample; j++) {
350	if (qmi.mask[l][j] != 0) {
351	tmpVal += (tmpSrc+stridej);
352	}
353	}
354	*dest++ = (float)(tmpVal / qmi.factor[l]);
355	}
356	}
357	} else {
358	for (int i=0; i<numSamples; i++) {
359	double tmpVal = 0.0;
360	for (int j=0; j<ptsPerSample; j++, src+=stride) {
361	tmpVal += *src;
362	}
363	tmpVal /= ptsPerSample;
364	for (int l=0; l<cellFactor; l++) {
365	*dest++ = static_cast<float>(tmpVal);
366	}
367	}
368	}
369	}
370	return res;
371	}
372
373	//
374	// Filters and reorders the raw sample values according to the node/element
375	// IDs. This is used to have data arrays ordered according to the underlying
376	// mesh (i.e. DataID[i]==MeshNodeID[i])
377	//
378	bool DataVar::reorderSamples()
379	{
380	if (numSamples == 0)
381	return true;
382
383	const IntVec* requiredIDs = NULL;
384	int requiredNumSamples = 0;
385	int cellFactor = 1;
386
387	if (centering == NODE_CENTERED) {
388	NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
389	requiredIDs = &nodes->getNodeIDs();
390	requiredNumSamples = nodes->getNumNodes();
391	} else {
392	ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
393	if (cells == NULL)
394	return false;
395
396	requiredIDs = &cells->getIDs();
397	requiredNumSamples = cells->getNumElements();
398	cellFactor = cells->getElementFactor();
399	if (cellFactor > 1) {
400	numSamples *= cellFactor;
401	// update sample IDs
402	IntVec newSampleID(numSamples);
403	IntVec::const_iterator idIt = sampleID.begin();
404	IntVec::iterator newIDit = newSampleID.begin();
405	for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
406	fill(newIDit, newIDit+cellFactor, *idIt);
407	}
408	sampleID.swap(newSampleID);
409	}
410	}
411
412	if (requiredNumSamples > numSamples) {
413	cerr << "ERROR: " << varName << " has " << numSamples
414	<< " instead of " << requiredNumSamples << " samples!" << endl;
415	return false;
416	}
417
418	IndexMap sampleID2idx = buildIndexMap();
419	numSamples = requiredNumSamples;
420
421	// now filter the data
422	for (size_t i=0; i < dataArray.size(); i++) {
423	float* c = new float[numSamples];
424	const float* src = dataArray[i];
425	IntVec::const_iterator idIt = requiredIDs->begin();
426	size_t destIdx = 0;
427	for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
428	size_t srcIdx = sampleID2idx.find(*idIt)->second;
429	copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
430	}
431	delete[] dataArray[i];
432	dataArray[i] = c;
433	}
434
435	// sample IDs now = mesh node/element IDs
436	sampleID = *requiredIDs;
437
438	return true;
439	}
440
441	//
442	//
443	//
444	int DataVar::getNumberOfComponents() const
445	{
446	return (rank == 0 ? 1 : accumulate(shape.begin(), shape.end(), 0));
447	}
448
449	//
450	//
451	//
452	float* DataVar::getDataFlat() const
453	{
454	int totalSize = numSamples * getNumberOfComponents();
455	float* res = new float[totalSize];
456	if (rank == 0) {
457	copy(dataArray[0], dataArray[0]+numSamples, res);
458	} else if (rank == 1) {
459	float *dest = res;
460	for (size_t c=0; c<numSamples; c++) {
461	for (size_t i=0; i<shape[0]; i++) {
462	*dest++ = dataArray[i][c];
463	}
464	}
465	} else if (rank == 2) {
466	float *dest = res;
467	for (size_t c=0; c<numSamples; c++) {
468	for (int i=0; i<shape[1]; i++) {
469	for (int j=0; j<shape[0]; j++) {
470	dest++ = dataArray[ishape[0]+j][c];
471	}
472	}
473	}
474	}
475
476	return res;
477	}
478
479	//
480	//
481	//
482	void DataVar::sampleToStream(ostream& os, int index)
483	{
484	// index is -1 for dummy samples, i.e. if writing the full mesh but
485	// only a reduced number of samples is required
486	if (rank == 0) {
487	if (index < 0) {
488	os << 0.;
489	} else {
490	os << dataArray[0][index];
491	}
492	} else if (rank == 1) {
493	if (index < 0) {
494	os << 0. << " " << 0. << " " << 0.;
495	} else if (shape[0] < 3) {
496	os << dataArray[0][index] << " " << dataArray[1][index]
497	<< " " << 0.;
498	} else {
499	os << dataArray[0][index] << " " << dataArray[1][index]
500	<< " " << dataArray[2][index];
501	}
502	} else if (rank == 2) {
503	if (index < 0) {
504	os << 0. << " " << 0. << " " << 0. << " ";
505	os << 0. << " " << 0. << " " << 0. << " ";
506	os << 0. << " " << 0. << " " << 0.;
507	} else if (shape[1] < 3) {
508	os << dataArray[0][index] << " " << dataArray[1][index]
509	<< " " << 0. << " ";
510	os << dataArray[2][index] << " " << dataArray[3][index]
511	<< " " << 0. << " ";
512	os << 0. << " " << 0. << " " << 0.;
513	} else {
514	os << dataArray[0][index] << " " << dataArray[1][index]
515	<< " " << dataArray[2][index] << " ";
516	os << dataArray[3][index] << " " << dataArray[4][index]
517	<< " " << dataArray[5][index] << " ";
518	os << dataArray[6][index] << " " << dataArray[7][index]
519	<< " " << dataArray[8][index];
520	}
521	}
522	os << endl;
523	}
524
525	//
526	//
527	//
528	void DataVar::writeToVTK(ostream& os, int ownIndex)
529	{
530	if (numSamples == 0)
531	return;
532
533	if (isNodeCentered()) {
534	// data was reordered in reorderSamples() but for VTK we write the
535	// original node mesh and thus need the original ordering.
536	// Note, that this also means we may not have samples for all nodes
537	// in which case we set idx to -1 and write a dummy sample
538	const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
539	const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
540	const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
541	int firstId = nodeDist[ownIndex];
542	int lastId = nodeDist[ownIndex+1];
543	IndexMap sampleID2idx = buildIndexMap();
544	for (int i=0; i<nodeGNI.size(); i++) {
545	if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
546	IndexMap::const_iterator it = sampleID2idx.find(requiredIDs[i]);
547	int idx = (it==sampleID2idx.end() ? -1 : (int)it->second);
548	sampleToStream(os, idx);
549	}
550	}
551	} else {
552	// cell data: ghost cells have been removed so do not write ghost
553	// samples (which are the last elements in the arrays)
554	int toWrite = domain->getElementsByName(meshName)->getNumElements();
555	for (int i=0; i<toWrite; i++) {
556	sampleToStream(os, i);
557	}
558	}
559	}
560
561	///////////////////////////////
562	// SILO related methods follow
563	///////////////////////////////
564
565	//
566	// If the data is tensor data then the components of the tensor are stored
567	// separately in the Silo file. This method then returns a string that
568	// contains the proper Silo expression to put the tensor together again.
569	// For non-tensor data this method returns an empty string.
570	//
571	string DataVar::getTensorDef() const
572	{
573	if (rank < 2 \|\| !initialized)
574	return string();
575
576	/// Format string for Silo 2x2 tensor
577	const string tensor2DefFmt =
578	"{{ <%sa_00>, <%sa_01> },"
579	" { <%sa_10>, <%sa_11> }}";
580
581	/// Format string for Silo 3x3 tensor
582	const string tensor3DefFmt =
583	"{{ <%sa_00>, <%sa_01>, <%sa_02> },"
584	" { <%sa_10>, <%sa_11>, <%sa_12> },"
585	" { <%sa_20>, <%sa_21>, <%sa_22> }}";
586
587	string tensorDef;
588	string tensorDir = varName+string("_comps/");
589	if (shape[1] == 3) {
590	char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
591	sprintf(tDef, tensor3DefFmt.c_str(),
592	tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
593	tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
594	tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
595	tensorDef = tDef;
596	delete[] tDef;
597	} else {
598	char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
599	sprintf(tDef, tensor2DefFmt.c_str(),
600	tensorDir.c_str(), tensorDir.c_str(),
601	tensorDir.c_str(), tensorDir.c_str(),
602	tensorDir.c_str(), tensorDir.c_str());
603	tensorDef = tDef;
604	delete[] tDef;
605	}
606	return tensorDef;
607	}
608
609	//
610	// Writes the data to given Silo file under the virtual path provided.
611	// The corresponding mesh must have been written already.
612	//
613	bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
614	const string& units)
615	{
616	#if USE_SILO
617	if (!initialized)
618	return false;
619
620	if (numSamples == 0)
621	return true;
622
623	int ret;
624
625	if (siloPath != "") {
626	ret = DBSetDir(dbfile, siloPath.c_str());
627	if (ret != 0)
628	return false;
629	}
630
631	char* siloMesh = const_cast<char*>(siloMeshName.c_str());
632	int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
633	DBoptlist* optList = DBMakeOptlist(2);
634	if (units.length()>0) {
635	DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
636	}
637
638	if (rank == 0) {
639	ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
640	numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
641	}
642	else if (rank == 1) {
643	const string comps[3] = {
644	varName+string("_x"), varName+string("_y"), varName+string("_z")
645	};
646	const char* varnames[3] = {
647	comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
648	};
649
650	ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
651	(char**)varnames, &dataArray[0], numSamples, NULL,
652	0, DB_FLOAT, dcenter, optList);
653	}
654	else {
655	string tensorDir = varName+string("_comps/");
656	ret = DBMkdir(dbfile, tensorDir.c_str());
657	if (ret == 0) {
658	int one = 1;
659	DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
660
661	for (int i=0; i<shape[1]; i++) {
662	for (int j=0; j<shape[0]; j++) {
663	ostringstream varname;
664	varname << tensorDir << "a_" << i << j;
665	ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
666	dataArray[i*shape[0]+j], numSamples,
667	NULL, 0, DB_FLOAT, dcenter, optList);
668	if (ret != 0) break;
669	}
670	if (ret != 0) break;
671	}
672	} // ret==0
673	} // rank
674
675	DBFreeOptlist(optList);
676	DBSetDir(dbfile, "/");
677	return (ret == 0);
678
679	#else // !USE_SILO
680	return false;
681	#endif
682	}
683
684	} // namespace weipa
685