ripley/src/Rectangle.cpp


/*******************************************************
*
* Copyright (c) 2003-2011 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 <ripley/Rectangle.h>
extern "C" {
#include "paso/SystemMatrixPattern.h"
}

#if USE_SILO
#include <silo.h>
#ifdef ESYS_MPI
#include <pmpio.h>
#endif
#endif

#include <iomanip>

using namespace std;

namespace ripley {

Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :
    RipleyDomain(2),
    m_gNE0(n0),
    m_gNE1(n1),
    m_l0(l0),
    m_l1(l1),
    m_NX(d0),
    m_NY(d1)
{
    // ensure number of subdivisions is valid and nodes can be distributed
    // among number of ranks
    if (m_NX*m_NY != m_mpiInfo->size)
        throw RipleyException("Invalid number of spatial subdivisions");

    if (n0%m_NX > 0 || n1%m_NY > 0)
        throw RipleyException("Number of elements must be separable into number of ranks in each dimension");

    // local number of elements
    m_NE0 = n0/m_NX;
    m_NE1 = n1/m_NY;
    // local number of nodes (not necessarily owned)
    m_N0 = m_NE0+1;
    m_N1 = m_NE1+1;
    // bottom-left node is at (offset0,offset1) in global mesh
    m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
    m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);
    populateSampleIds();
}

Rectangle::~Rectangle()
{
}

string Rectangle::getDescription() const
{
    return "ripley::Rectangle";
}

bool Rectangle::operator==(const AbstractDomain& other) const
{
    if (dynamic_cast<const Rectangle*>(&other))
        return this==&other;

    return false;
}

void Rectangle::dump(const string& fileName) const
{
#if USE_SILO
    string fn(fileName);
    if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
        fn+=".silo";
    }

    const int NUM_SILO_FILES = 1;
    const char* blockDirFmt = "/block%04d";
    int driver=DB_HDF5;    
    string siloPath;
    DBfile* dbfile = NULL;

#ifdef ESYS_MPI
    PMPIO_baton_t* baton = NULL;
#endif

    if (m_mpiInfo->size > 1) {
#ifdef ESYS_MPI
        baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
                    0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
                    PMPIO_DefaultClose, (void*)&driver);
        // try the fallback driver in case of error
        if (!baton && driver != DB_PDB) {
            driver = DB_PDB;
            baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
                        0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
                        PMPIO_DefaultClose, (void*)&driver);
        }
        if (baton) {
            char str[64];
            snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
            siloPath = str;
            dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
        }
#endif
    } else {
        dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
                getDescription().c_str(), driver);
        // try the fallback driver in case of error
        if (!dbfile && driver != DB_PDB) {
            driver = DB_PDB;
            dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
                    getDescription().c_str(), driver);
        }
    }

    if (!dbfile)
        throw RipleyException("dump: Could not create Silo file");

    /*
    if (driver==DB_HDF5) {
        // gzip level 1 already provides good compression with minimal
        // performance penalty. Some tests showed that gzip levels >3 performed
        // rather badly on escript data both in terms of time and space
        DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
    }
    */

    boost::scoped_ptr<double> x(new double[m_N0]);
    boost::scoped_ptr<double> y(new double[m_N1]);
    double* coords[2] = { x.get(), y.get() };
#pragma omp parallel
    {
#pragma omp for
        for (dim_t i0 = 0; i0 < m_N0; i0++) {
            coords[0][i0]=(m_l0*(i0+m_offset0))/m_gNE0;
        }
#pragma omp for
        for (dim_t i1 = 0; i1 < m_N1; i1++) {
            coords[1][i1]=(m_l1*(i1+m_offset1))/m_gNE1;
        }
    }
    int dims[] = { m_N0, m_N1 };
    DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
            DB_COLLINEAR, NULL);

    if (m_mpiInfo->rank == 0) {
        vector<string> tempstrings;
        vector<char*> meshnames;
        for (dim_t i=0; i<m_mpiInfo->size; i++) {
            stringstream path;
            path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
            tempstrings.push_back(path.str());
            meshnames.push_back((char*)tempstrings.back().c_str());
        }
        vector<int> meshtypes(m_mpiInfo->size, DB_QUAD_RECT);
        DBSetDir(dbfile, "/");
        DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &meshnames[0],
               &meshtypes[0], NULL);
    }

    if (m_mpiInfo->size > 1) {
#ifdef ESYS_MPI
        PMPIO_HandOffBaton(baton, dbfile);
        PMPIO_Finish(baton);
#endif
    } else {
        DBClose(dbfile);
    }

#else // USE_SILO
    throw RipleyException("dump(): no Silo support");
#endif
}

const int* Rectangle::borrowSampleReferenceIDs(int functionSpaceType) const
{
    switch (functionSpaceType) {
        case Nodes:
            return &m_nodeId[0];
        case Elements:
            return &m_elementId[0];
        case FaceElements:
            return &m_faceId[0];
        default:
            break;
    }

    stringstream msg;
    msg << "borrowSampleReferenceIDs() not implemented for function space type "
        << functionSpaceType;
    throw RipleyException(msg.str());
}

bool Rectangle::ownSample(int fsCode, index_t id) const
{
#ifdef ESYS_MPI
    if (fsCode == Nodes) {
        const index_t myFirst=getNumNodes()*m_mpiInfo->rank;
        const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;
        return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
    } else
        throw RipleyException("ownSample() only implemented for Nodes");
#else
    return true;
#endif
}

void Rectangle::interpolateOnDomain(escript::Data& target,
                                    const escript::Data& in) const
{
    const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));
    const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));
    if (inDomain != *this)
        throw RipleyException("Illegal domain of interpolant");
    if (targetDomain != *this)
        throw RipleyException("Illegal domain of interpolation target");

    throw RipleyException("interpolateOnDomain() not implemented");
}

Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
                                                bool reducedColOrder) const
{
    if (reducedRowOrder || reducedColOrder)
        throw RipleyException("getPattern() not implemented for reduced order");

/*
    // create distribution
    IndexVector dist;
    for (index_t i=0; i<m_mpiInfo->size+1; i++)
        dist.push_back(i*getNumNodes());
    Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
            &dist[0], 1, 0);

    // connectors
    dim_t numNeighbours = ...;
    RankVector neighbour(numNeighbours);
    IndexVector offsetInShared(numNeighbours+1);
    IndexVector shared(offsetInShared[numNeighbours]);

    Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
            getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,
            1, 0, m_mpiInfo);
    Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
            getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,
            1, 0, m_mpiInfo);
    Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);

    // create patterns
    Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
            M, N, ptr, index);
    Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
            M, N, ...);
    Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
            ...);

    Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
            MATRIX_FORMAT_DEFAULT, distribution, distribution,
            mainPattern, colCouplePattern, rowCouplePattern,
            connector, connector);
    Paso_Pattern_free(mainPattern);
    Paso_Pattern_free(colCouplePattern);
    Paso_Pattern_free(rowCouplePattern);
    Paso_Distribution_free(distribution);
    Paso_SharedComponents_free(snd_shcomp);
    Paso_SharedComponents_free(rcv_shcomp);
*/
    throw RipleyException("getPattern() not implemented");
}

void Rectangle::Print_Mesh_Info(const bool full) const
{
    RipleyDomain::Print_Mesh_Info(full);
    if (full) {
        cout << "     Id  Coordinates" << endl;
        cout.precision(15);
        cout.setf(ios::scientific, ios::floatfield);
        for (index_t i=0; i < getNumNodes(); i++) {
            cout << "  " << setw(5) << m_nodeId[i]
                << "  " << (m_l0*(i%m_N0+m_offset0))/m_gNE0
                << "  " << (m_l1*(i/m_N0+m_offset1))/m_gNE1 << endl;
        }
    }
}

//protected
dim_t Rectangle::getNumFaceElements() const
{
    dim_t n=0;
    //left
    if (m_offset0==0)
        n+=m_NE1;
    //right
    if (m_mpiInfo->rank%m_NX==m_NX-1)
        n+=m_NE1;
    //bottom
    if (m_offset1==0)
        n+=m_NE0;
    //top
    if (m_mpiInfo->rank/m_NX==m_NY-1)
        n+=m_NE0;

    return n;
}

//protected
void Rectangle::assembleCoordinates(escript::Data& arg) const
{
    escriptDataC x = arg.getDataC();
    int numDim = m_numDim;
    if (!isDataPointShapeEqual(&x, 1, &numDim))
        throw RipleyException("setToX: Invalid Data object shape");
    if (!numSamplesEqual(&x, 1, getNumNodes()))
        throw RipleyException("setToX: Illegal number of samples in Data object");

    arg.requireWrite();
#pragma omp parallel for
    for (dim_t i1 = 0; i1 < m_N1; i1++) {
        for (dim_t i0 = 0; i0 < m_N0; i0++) {
            double* point = arg.getSampleDataRW(i0+m_N0*i1);
            point[0] = (m_l0*(i0+m_offset0))/m_gNE0;
            point[1] = (m_l1*(i1+m_offset1))/m_gNE1;
        }
    }
}

//private
void Rectangle::populateSampleIds()
{
    const index_t firstId = getNumNodes()*m_mpiInfo->rank;
    const index_t diff0 = m_N0*(m_N1-1)+1;
    const index_t diff1 = m_N0*((m_NX-1)*m_N1+1);
    m_nodeId.resize(getNumNodes());
#pragma omp parallel for
    for (dim_t k=0; k<getNumNodes(); k++) {
        index_t id = firstId+k;
        if (m_offset0 > 0 && k%m_N0==0)
            id -= diff0;
        if (m_offset1 > 0 && k/m_N0==0)
            id -= diff1;
        m_nodeId[k]=id;
    }
}

} // end of namespace ripley

1	caltinay	3691
2			/*******************************************************
3			*
4			* Copyright (c) 2003-2011 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 <ripley/Rectangle.h>
15			extern "C" {
16			#include "paso/SystemMatrixPattern.h"
17			}
18
19			#if USE_SILO
20			#include <silo.h>
21			#ifdef ESYS_MPI
22			#include <pmpio.h>
23			#endif
24			#endif
25
26			#include <iomanip>
27
28			using namespace std;
29
30			namespace ripley {
31
32			Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :
33			RipleyDomain(2),
34			m_gNE0(n0),
35			m_gNE1(n1),
36			m_l0(l0),
37			m_l1(l1),
38			m_NX(d0),
39			m_NY(d1)
40			{
41			// ensure number of subdivisions is valid and nodes can be distributed
42			// among number of ranks
43			if (m_NX*m_NY != m_mpiInfo->size)
44			throw RipleyException("Invalid number of spatial subdivisions");
45
46			if (n0%m_NX > 0 \|\| n1%m_NY > 0)
47			throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
48
49			// local number of elements
50			m_NE0 = n0/m_NX;
51			m_NE1 = n1/m_NY;
52			// local number of nodes (not necessarily owned)
53			m_N0 = m_NE0+1;
54			m_N1 = m_NE1+1;
55			// bottom-left node is at (offset0,offset1) in global mesh
56			m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
57			m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);
58			populateSampleIds();
59			}
60
61			Rectangle::~Rectangle()
62			{
63			}
64
65			string Rectangle::getDescription() const
66			{
67			return "ripley::Rectangle";
68			}
69
70			bool Rectangle::operator==(const AbstractDomain& other) const
71			{
72			if (dynamic_cast<const Rectangle*>(&other))
73			return this==&other;
74
75			return false;
76			}
77
78			void Rectangle::dump(const string& fileName) const
79			{
80			#if USE_SILO
81			string fn(fileName);
82			if (fileName.length() < 6 \|\| fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
83			fn+=".silo";
84			}
85
86			const int NUM_SILO_FILES = 1;
87			const char* blockDirFmt = "/block%04d";
88			int driver=DB_HDF5;
89			string siloPath;
90			DBfile* dbfile = NULL;
91
92			#ifdef ESYS_MPI
93			PMPIO_baton_t* baton = NULL;
94			#endif
95
96			if (m_mpiInfo->size > 1) {
97			#ifdef ESYS_MPI
98			baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
99			0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
100			PMPIO_DefaultClose, (void*)&driver);
101			// try the fallback driver in case of error
102			if (!baton && driver != DB_PDB) {
103			driver = DB_PDB;
104			baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
105			0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
106			PMPIO_DefaultClose, (void*)&driver);
107			}
108			if (baton) {
109			char str[64];
110			snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
111			siloPath = str;
112			dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
113			}
114			#endif
115			} else {
116			dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
117			getDescription().c_str(), driver);
118			// try the fallback driver in case of error
119			if (!dbfile && driver != DB_PDB) {
120			driver = DB_PDB;
121			dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
122			getDescription().c_str(), driver);
123			}
124			}
125
126			if (!dbfile)
127			throw RipleyException("dump: Could not create Silo file");
128
129			/*
130			if (driver==DB_HDF5) {
131			// gzip level 1 already provides good compression with minimal
132			// performance penalty. Some tests showed that gzip levels >3 performed
133			// rather badly on escript data both in terms of time and space
134			DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
135			}
136			*/
137
138			boost::scoped_ptr<double> x(new double[m_N0]);
139			boost::scoped_ptr<double> y(new double[m_N1]);
140			double* coords[2] = { x.get(), y.get() };
141			#pragma omp parallel
142			{
143			#pragma omp for
144			for (dim_t i0 = 0; i0 < m_N0; i0++) {
145			coords[0][i0]=(m_l0*(i0+m_offset0))/m_gNE0;
146			}
147			#pragma omp for
148			for (dim_t i1 = 0; i1 < m_N1; i1++) {
149			coords[1][i1]=(m_l1*(i1+m_offset1))/m_gNE1;
150			}
151			}
152			int dims[] = { m_N0, m_N1 };
153			DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
154			DB_COLLINEAR, NULL);
155
156			if (m_mpiInfo->rank == 0) {
157			vector<string> tempstrings;
158			vector<char*> meshnames;
159			for (dim_t i=0; i<m_mpiInfo->size; i++) {
160			stringstream path;
161			path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
162			tempstrings.push_back(path.str());
163			meshnames.push_back((char*)tempstrings.back().c_str());
164			}
165			vector<int> meshtypes(m_mpiInfo->size, DB_QUAD_RECT);
166			DBSetDir(dbfile, "/");
167			DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &meshnames[0],
168			&meshtypes[0], NULL);
169			}
170
171			if (m_mpiInfo->size > 1) {
172			#ifdef ESYS_MPI
173			PMPIO_HandOffBaton(baton, dbfile);
174			PMPIO_Finish(baton);
175			#endif
176			} else {
177			DBClose(dbfile);
178			}
179
180			#else // USE_SILO
181			throw RipleyException("dump(): no Silo support");
182			#endif
183			}
184
185			const int* Rectangle::borrowSampleReferenceIDs(int functionSpaceType) const
186			{
187			switch (functionSpaceType) {
188			case Nodes:
189			return &m_nodeId[0];
190			case Elements:
191			return &m_elementId[0];
192			case FaceElements:
193			return &m_faceId[0];
194			default:
195			break;
196			}
197
198			stringstream msg;
199			msg << "borrowSampleReferenceIDs() not implemented for function space type "
200			<< functionSpaceType;
201			throw RipleyException(msg.str());
202			}
203
204			bool Rectangle::ownSample(int fsCode, index_t id) const
205			{
206			#ifdef ESYS_MPI
207			if (fsCode == Nodes) {
208			const index_t myFirst=getNumNodes()*m_mpiInfo->rank;
209			const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;
210			return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
211			} else
212			throw RipleyException("ownSample() only implemented for Nodes");
213			#else
214			return true;
215			#endif
216			}
217
218			void Rectangle::interpolateOnDomain(escript::Data& target,
219			const escript::Data& in) const
220			{
221			const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));
222			const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));
223			if (inDomain != *this)
224			throw RipleyException("Illegal domain of interpolant");
225			if (targetDomain != *this)
226			throw RipleyException("Illegal domain of interpolation target");
227
228			throw RipleyException("interpolateOnDomain() not implemented");
229			}
230
231			Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
232			bool reducedColOrder) const
233			{
234			if (reducedRowOrder \|\| reducedColOrder)
235			throw RipleyException("getPattern() not implemented for reduced order");
236
237			/*
238			// create distribution
239			IndexVector dist;
240			for (index_t i=0; i<m_mpiInfo->size+1; i++)
241			dist.push_back(i*getNumNodes());
242			Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
243			&dist[0], 1, 0);
244
245			// connectors
246			dim_t numNeighbours = ...;
247			RankVector neighbour(numNeighbours);
248			IndexVector offsetInShared(numNeighbours+1);
249			IndexVector shared(offsetInShared[numNeighbours]);
250
251			Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
252			getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,
253			1, 0, m_mpiInfo);
254			Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
255			getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,
256			1, 0, m_mpiInfo);
257			Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
258
259			// create patterns
260			Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
261			M, N, ptr, index);
262			Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
263			M, N, ...);
264			Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
265			...);
266
267			Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
268			MATRIX_FORMAT_DEFAULT, distribution, distribution,
269			mainPattern, colCouplePattern, rowCouplePattern,
270			connector, connector);
271			Paso_Pattern_free(mainPattern);
272			Paso_Pattern_free(colCouplePattern);
273			Paso_Pattern_free(rowCouplePattern);
274			Paso_Distribution_free(distribution);
275			Paso_SharedComponents_free(snd_shcomp);
276			Paso_SharedComponents_free(rcv_shcomp);
277			*/
278			throw RipleyException("getPattern() not implemented");
279			}
280
281			void Rectangle::Print_Mesh_Info(const bool full) const
282			{
283			RipleyDomain::Print_Mesh_Info(full);
284			if (full) {
285			cout << " Id Coordinates" << endl;
286			cout.precision(15);
287			cout.setf(ios::scientific, ios::floatfield);
288			for (index_t i=0; i < getNumNodes(); i++) {
289			cout << " " << setw(5) << m_nodeId[i]
290			<< " " << (m_l0*(i%m_N0+m_offset0))/m_gNE0
291			<< " " << (m_l1*(i/m_N0+m_offset1))/m_gNE1 << endl;
292			}
293			}
294			}
295
296			//protected
297			dim_t Rectangle::getNumFaceElements() const
298			{
299			dim_t n=0;
300			//left
301			if (m_offset0==0)
302			n+=m_NE1;
303			//right
304			if (m_mpiInfo->rank%m_NX==m_NX-1)
305			n+=m_NE1;
306			//bottom
307			if (m_offset1==0)
308			n+=m_NE0;
309			//top
310			if (m_mpiInfo->rank/m_NX==m_NY-1)
311			n+=m_NE0;
312
313			return n;
314			}
315
316			//protected
317			void Rectangle::assembleCoordinates(escript::Data& arg) const
318			{
319			escriptDataC x = arg.getDataC();
320			int numDim = m_numDim;
321			if (!isDataPointShapeEqual(&x, 1, &numDim))
322			throw RipleyException("setToX: Invalid Data object shape");
323			if (!numSamplesEqual(&x, 1, getNumNodes()))
324			throw RipleyException("setToX: Illegal number of samples in Data object");
325
326			arg.requireWrite();
327			#pragma omp parallel for
328			for (dim_t i1 = 0; i1 < m_N1; i1++) {
329			for (dim_t i0 = 0; i0 < m_N0; i0++) {
330			double* point = arg.getSampleDataRW(i0+m_N0*i1);
331			point[0] = (m_l0*(i0+m_offset0))/m_gNE0;
332			point[1] = (m_l1*(i1+m_offset1))/m_gNE1;
333			}
334			}
335			}
336
337			//private
338			void Rectangle::populateSampleIds()
339			{
340			const index_t firstId = getNumNodes()*m_mpiInfo->rank;
341			const index_t diff0 = m_N0*(m_N1-1)+1;
342			const index_t diff1 = m_N0((m_NX-1)m_N1+1);
343			m_nodeId.resize(getNumNodes());
344			#pragma omp parallel for
345			for (dim_t k=0; k<getNumNodes(); k++) {
346			index_t id = firstId+k;
347			if (m_offset0 > 0 && k%m_N0==0)
348			id -= diff0;
349			if (m_offset1 > 0 && k/m_N0==0)
350			id -= diff1;
351			m_nodeId[k]=id;
352			}
353			}
354
355			} // end of namespace ripley
356