escript/src/DataFactory.cpp


/* $Id$ */

/*******************************************************
 *
 *           Copyright 2003-2007 by ACceSS MNRF
 *       Copyright 2007 by University of Queensland
 *
 *                http://esscc.uq.edu.au
 *        Primary Business: Queensland, Australia
 *  Licensed under the Open Software License version 3.0
 *     http://www.opensource.org/licenses/osl-3.0.php
 *
 *******************************************************/

#include "DataFactory.h"
#include "esysUtils/esys_malloc.h"

#include <boost/python/extract.hpp>
#include <iostream>
#include <exception>
#ifdef USE_NETCDF
#include <netcdfcpp.h>
#endif

using namespace boost::python;

namespace escript {

Data
Scalar(double value,
       const FunctionSpace& what,
       bool expanded)
{
    //
    // an empty shape is a scalar
    DataTypes::ShapeType shape;
    return Data(value,shape,what,expanded);
}

Data
Vector(double value,
       const FunctionSpace& what,
       bool expanded)
{
    DataTypes::ShapeType shape(1,what.getDomain().getDim());
    return Data(value,shape,what,expanded);
}

Data
Tensor(double value,
       const FunctionSpace& what,
       bool expanded)
{
    DataTypes::ShapeType shape(2,what.getDomain().getDim());
    return Data(value,shape,what,expanded);
}

Data
Tensor3(double value,
        const FunctionSpace& what,
        bool expanded)
{
    DataTypes::ShapeType shape(3,what.getDomain().getDim());
    return Data(value,shape,what,expanded);
}

Data
Tensor4(double value,
        const FunctionSpace& what,
        bool expanded)
{
    DataTypes::ShapeType shape(4,what.getDomain().getDim());
    return Data(value,shape,what,expanded);
}

Data 
load(const std::string fileName,
     const AbstractDomain& domain)
{
   #ifdef PASO_MPI
   throw DataException("Error - load :: not implemented for MPI yet.");
   #endif

   #ifdef USE_NETCDF
   NcAtt *type_att, *rank_att, *function_space_type_att;
   // netCDF error handler
   NcError err(NcError::silent_nonfatal);
   // Create the file.
   NcFile dataFile(fileName.c_str(), NcFile::ReadOnly);
   if (!dataFile.is_valid())
        throw DataException("Error - load:: opening of netCDF file for input failed.");
   /* recover function space */
   if (! (function_space_type_att=dataFile.get_att("function_space_type")) )
        throw DataException("Error - load:: cannot recover function_space_type attribute from escript netCDF file.");
   int function_space_type = function_space_type_att->as_int(0);
   delete function_space_type_att;
   /* test if function space id is valid and create function space instance */
   if (! domain.isValidFunctionSpaceType(function_space_type) ) 
        throw DataException("Error - load:: function space type code in netCDF file is invalid for given domain.");
   FunctionSpace function_space=FunctionSpace(domain, function_space_type);
   /* recover rank */
   if (! (rank_att=dataFile.get_att("rank")) )
        throw DataException("Error - load:: cannot recover rank attribute from escript netCDF file.");
   int rank = rank_att->as_int(0);
   delete rank_att;
   if (rank<0 || rank>DataTypes::maxRank)
        throw DataException("Error - load:: rank in escript netCDF file is greater than maximum rank.");
   /* recover type attribute */
   int type=-1;
   if ((type_att=dataFile.get_att("type")) ) {
       char* type_str = type_att->as_string(0);
       if (strncmp(type_str, "constant", strlen("constant")) == 0 ) {
          type =0;
       } else if (strncmp(type_str, "tagged", strlen("tagged")) == 0 ) {
           type =1;
       } else if (strncmp(type_str, "expanded", strlen("expanded")) == 0 ) {
           type =2;
       }
       esysUtils::free(type_str);
   } else {
      if (! (type_att=dataFile.get_att("type_id")) )
    throw DataException("Error - load:: cannot recover type attribute from escript netCDF file.");
      type=type_att->as_int(0);
   }
   delete type_att;

   /* recover dimension */
   int ndims=dataFile.num_dims();
   int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
   NcDim *d_dim, *tags_dim, *num_samples_dim, *num_data_points_per_sample_dim;
   /* recover shape */
   DataTypes::ShapeType shape;
   long dims[DataTypes::maxRank+2];
   if (rank>0) {
     if (! (d_dim=dataFile.get_dim("d0")) )
          throw DataException("Error - load:: unable to recover d0 from netCDF file.");
      d=d_dim->size();
      shape.push_back(d);
      dims[0]=d;
      len_data_point*=d;
   }
   if (rank>1) {
     if (! (d_dim=dataFile.get_dim("d1")) )
          throw DataException("Error - load:: unable to recover d1 from netCDF file.");
      d=d_dim->size();
      shape.push_back(d);
      dims[1]=d;
      len_data_point*=d;
   }
   if (rank>2) {
     if (! (d_dim=dataFile.get_dim("d2")) )
          throw DataException("Error - load:: unable to recover d2 from netCDF file.");
      d=d_dim->size();
      shape.push_back(d);
      dims[2]=d;
      len_data_point*=d;
   }
   if (rank>3) {
     if (! (d_dim=dataFile.get_dim("d3")) )
          throw DataException("Error - load:: unable to recover d3 from netCDF file.");
      d=d_dim->size();
      shape.push_back(d);
      dims[3]=d;
      len_data_point*=d;
   }
   /* recover stuff */
   Data out;
   NcVar *var, *ids_var, *tags_var;
   if (type == 0) {
      /* constant data */
      if ( ! ( (ndims == rank && rank >0) || ( ndims ==1 && rank == 0 ) ) )
          throw DataException("Error - load:: illegal number of dimensions for constant data in netCDF file.");
      if (rank == 0) {
          if (! (d_dim=dataFile.get_dim("l")) )
              throw DataException("Error - load:: unable to recover d0 for scalar constant data in netCDF file.");
          int d0 = d_dim->size();
          if (! d0 == 1) 
              throw DataException("Error - load:: d0 is expected to be one for scalar constant data in netCDF file.");
          dims[0]=1;
      }
      out=Data(0,shape,function_space);
      if (!(var = dataFile.get_var("data")))
              throw DataException("Error - load:: unable to find data in netCDF file.");
      if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) ) 
              throw DataException("Error - load:: unable to recover data from netCDF file.");
   } else if (type == 1) { 
      /* tagged data */
      if ( ! (ndims == rank + 1) )
         throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
      if (! (tags_dim=dataFile.get_dim("num_tags")) )
         throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
      ntags=tags_dim->size();
      dims[rank]=ntags;
      int *tags = (int *) esysUtils::malloc(ntags*sizeof(int));
      if (! ( tags_var = dataFile.get_var("tags")) )
      {
         esysUtils::free(tags);
         throw DataException("Error - load:: unable to find tags in netCDF file.");
      }
      if (! tags_var->get(tags, ntags) ) 
      {
         esysUtils::free(tags);
         throw DataException("Error - load:: unable to recover tags from netCDF file.");
      }

      DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
      if (!(var = dataFile.get_var("data")))
      {
         esysUtils::free(tags);
         throw DataException("Error - load:: unable to find data in netCDF file.");
      }
      if (! var->get(&(data[0]), dims) ) 
      {
         esysUtils::free(tags);
         throw DataException("Error - load:: unable to recover data from netCDF file.");
      }
      out=Data(DataArrayView(data,shape,0),function_space);
      for (int t=1; t<ntags; ++t) {
     out.setTaggedValueFromCPP(tags[t],shape, data, t*len_data_point);
//         out.setTaggedValueFromCPP(tags[t],DataArrayView(data,shape,t*len_data_point));
      }
      esysUtils::free(tags);
   } else if (type == 2) {
      /* expanded data */
      if ( ! (ndims == rank + 2) )
          throw DataException("Error - load:: illegal number of dimensions for exanded data in netCDF file.");
      if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
          throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
      num_samples = num_samples_dim->size();
      if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
          throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
      num_data_points_per_sample=num_data_points_per_sample_dim->size();
      // check shape:
      if ( ! (num_samples == function_space.getNumSamples() && num_data_points_per_sample == function_space.getNumDataPointsPerSample()) )
          throw DataException("Error - load:: data sample layout of file does not match data layout of function space.");
      // get ids
      if (! ( ids_var = dataFile.get_var("id")) )
         throw DataException("Error - load:: unable to find reference ids in netCDF file.");
      const int* ids_p=function_space.borrowSampleReferenceIDs();
      int *ids_of_nc = (int *)esysUtils::malloc(num_samples*sizeof(int));
      if (! ids_var->get(ids_of_nc, (long) num_samples) ) 
      {
         esysUtils::free(ids_of_nc);
         throw DataException("Error - load:: unable to recover ids from netCDF file.");
      }
      // check order:
      int failed=-1, local_failed=-1, i;
      #pragma omp parallel private(local_failed)
      {
          local_failed=-1;
          #pragma omp for private(i) schedule(static)
          for (i=0;i < num_samples; ++i) {
              if (ids_of_nc[i]!=ids_p[i]) local_failed=i;
          }
          #pragma omp critical
          if (local_failed>=0) failed = local_failed;
      }
      /* if (failed>=0) 
      {
         esysUtils::free(ids_of_nc);
         throw DataException("Error - load:: data ordering in netCDF file does not match ordering of FunctionSpace.");
      } */
      // get the data:
      dims[rank]=num_data_points_per_sample;
      dims[rank+1]=num_samples;
      out=Data(0,shape,function_space,true);
      if (!(var = dataFile.get_var("data")))
      {
         esysUtils::free(ids_of_nc);
         throw DataException("Error - load:: unable to find data in netCDF file.");
      }
      if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) ) 
      {
         esysUtils::free(ids_of_nc);
         throw DataException("Error - load:: unable to recover data from netCDF file.");
      }
      if (failed>=0) {
        try {
           std::cout << "Information - load: start reordering data from netCDF file " << fileName << std::endl;
           out.borrowData()->reorderByReferenceIDs(ids_of_nc);
        } 
        catch (std::exception&) {
           esysUtils::free(ids_of_nc);
           throw DataException("Error - load:: unable to reorder data in netCDF file.");
        }
      }
   } else {
       throw DataException("Error - load:: unknown escript data type in netCDF file.");
   }
   return out;
   #else
   throw DataException("Error - load:: is not compiled with netCFD. Please contact your insatllation manager.");
   #endif
}

bool 
loadConfigured()
{
   #ifdef USE_NETCDF
   return true;
   #else
   return false;
   #endif
}

Data
convertToData(const boost::python::object& value,
              const FunctionSpace& what) 
{
     // first we try to extract a Data object from value 
     extract<Data> value_data(value);
     if (value_data.check()) {
         Data extracted_data=value_data();
         if (extracted_data.isEmpty()) {
            return extracted_data;
         } else {
            return Data(extracted_data,what);
         }
     } else {
        return Data(value,what);
     }
}

}  // end of namespace
1
2	/* $Id$ */
3
4	/*******************************************************
5	*
6	* Copyright 2003-2007 by ACceSS MNRF
7	* Copyright 2007 by University of Queensland
8	*
9	* http://esscc.uq.edu.au
10	* Primary Business: Queensland, Australia
11	* Licensed under the Open Software License version 3.0
12	* http://www.opensource.org/licenses/osl-3.0.php
13	*
14	*******************************************************/
15
16	#include "DataFactory.h"
17	#include "esysUtils/esys_malloc.h"
18
19	#include <boost/python/extract.hpp>
20	#include <iostream>
21	#include <exception>
22	#ifdef USE_NETCDF
23	#include <netcdfcpp.h>
24	#endif
25
26	using namespace boost::python;
27
28	namespace escript {
29
30	Data
31	Scalar(double value,
32	const FunctionSpace& what,
33	bool expanded)
34	{
35	//
36	// an empty shape is a scalar
37	DataTypes::ShapeType shape;
38	return Data(value,shape,what,expanded);
39	}
40
41	Data
42	Vector(double value,
43	const FunctionSpace& what,
44	bool expanded)
45	{
46	DataTypes::ShapeType shape(1,what.getDomain().getDim());
47	return Data(value,shape,what,expanded);
48	}
49
50	Data
51	Tensor(double value,
52	const FunctionSpace& what,
53	bool expanded)
54	{
55	DataTypes::ShapeType shape(2,what.getDomain().getDim());
56	return Data(value,shape,what,expanded);
57	}
58
59	Data
60	Tensor3(double value,
61	const FunctionSpace& what,
62	bool expanded)
63	{
64	DataTypes::ShapeType shape(3,what.getDomain().getDim());
65	return Data(value,shape,what,expanded);
66	}
67
68	Data
69	Tensor4(double value,
70	const FunctionSpace& what,
71	bool expanded)
72	{
73	DataTypes::ShapeType shape(4,what.getDomain().getDim());
74	return Data(value,shape,what,expanded);
75	}
76
77	Data
78	load(const std::string fileName,
79	const AbstractDomain& domain)
80	{
81	#ifdef PASO_MPI
82	throw DataException("Error - load :: not implemented for MPI yet.");
83	#endif
84
85	#ifdef USE_NETCDF
86	NcAtt type_att, rank_att, *function_space_type_att;
87	// netCDF error handler
88	NcError err(NcError::silent_nonfatal);
89	// Create the file.
90	NcFile dataFile(fileName.c_str(), NcFile::ReadOnly);
91	if (!dataFile.is_valid())
92	throw DataException("Error - load:: opening of netCDF file for input failed.");
93	/* recover function space */
94	if (! (function_space_type_att=dataFile.get_att("function_space_type")) )
95	throw DataException("Error - load:: cannot recover function_space_type attribute from escript netCDF file.");
96	int function_space_type = function_space_type_att->as_int(0);
97	delete function_space_type_att;
98	/* test if function space id is valid and create function space instance */
99	if (! domain.isValidFunctionSpaceType(function_space_type) )
100	throw DataException("Error - load:: function space type code in netCDF file is invalid for given domain.");
101	FunctionSpace function_space=FunctionSpace(domain, function_space_type);
102	/* recover rank */
103	if (! (rank_att=dataFile.get_att("rank")) )
104	throw DataException("Error - load:: cannot recover rank attribute from escript netCDF file.");
105	int rank = rank_att->as_int(0);
106	delete rank_att;
107	if (rank<0 \|\| rank>DataTypes::maxRank)
108	throw DataException("Error - load:: rank in escript netCDF file is greater than maximum rank.");
109	/* recover type attribute */
110	int type=-1;
111	if ((type_att=dataFile.get_att("type")) ) {
112	char* type_str = type_att->as_string(0);
113	if (strncmp(type_str, "constant", strlen("constant")) == 0 ) {
114	type =0;
115	} else if (strncmp(type_str, "tagged", strlen("tagged")) == 0 ) {
116	type =1;
117	} else if (strncmp(type_str, "expanded", strlen("expanded")) == 0 ) {
118	type =2;
119	}
120	esysUtils::free(type_str);
121	} else {
122	if (! (type_att=dataFile.get_att("type_id")) )
123	throw DataException("Error - load:: cannot recover type attribute from escript netCDF file.");
124	type=type_att->as_int(0);
125	}
126	delete type_att;
127
128	/* recover dimension */
129	int ndims=dataFile.num_dims();
130	int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
131	NcDim d_dim, tags_dim, num_samples_dim, num_data_points_per_sample_dim;
132	/* recover shape */
133	DataTypes::ShapeType shape;
134	long dims[DataTypes::maxRank+2];
135	if (rank>0) {
136	if (! (d_dim=dataFile.get_dim("d0")) )
137	throw DataException("Error - load:: unable to recover d0 from netCDF file.");
138	d=d_dim->size();
139	shape.push_back(d);
140	dims[0]=d;
141	len_data_point*=d;
142	}
143	if (rank>1) {
144	if (! (d_dim=dataFile.get_dim("d1")) )
145	throw DataException("Error - load:: unable to recover d1 from netCDF file.");
146	d=d_dim->size();
147	shape.push_back(d);
148	dims[1]=d;
149	len_data_point*=d;
150	}
151	if (rank>2) {
152	if (! (d_dim=dataFile.get_dim("d2")) )
153	throw DataException("Error - load:: unable to recover d2 from netCDF file.");
154	d=d_dim->size();
155	shape.push_back(d);
156	dims[2]=d;
157	len_data_point*=d;
158	}
159	if (rank>3) {
160	if (! (d_dim=dataFile.get_dim("d3")) )
161	throw DataException("Error - load:: unable to recover d3 from netCDF file.");
162	d=d_dim->size();
163	shape.push_back(d);
164	dims[3]=d;
165	len_data_point*=d;
166	}
167	/* recover stuff */
168	Data out;
169	NcVar var, ids_var, *tags_var;
170	if (type == 0) {
171	/* constant data */
172	if ( ! ( (ndims == rank && rank >0) \|\| ( ndims ==1 && rank == 0 ) ) )
173	throw DataException("Error - load:: illegal number of dimensions for constant data in netCDF file.");
174	if (rank == 0) {
175	if (! (d_dim=dataFile.get_dim("l")) )
176	throw DataException("Error - load:: unable to recover d0 for scalar constant data in netCDF file.");
177	int d0 = d_dim->size();
178	if (! d0 == 1)
179	throw DataException("Error - load:: d0 is expected to be one for scalar constant data in netCDF file.");
180	dims[0]=1;
181	}
182	out=Data(0,shape,function_space);
183	if (!(var = dataFile.get_var("data")))
184	throw DataException("Error - load:: unable to find data in netCDF file.");
185	if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) )
186	throw DataException("Error - load:: unable to recover data from netCDF file.");
187	} else if (type == 1) {
188	/* tagged data */
189	if ( ! (ndims == rank + 1) )
190	throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
191	if (! (tags_dim=dataFile.get_dim("num_tags")) )
192	throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
193	ntags=tags_dim->size();
194	dims[rank]=ntags;
195	int tags = (int ) esysUtils::malloc(ntags*sizeof(int));
196	if (! ( tags_var = dataFile.get_var("tags")) )
197	{
198	esysUtils::free(tags);
199	throw DataException("Error - load:: unable to find tags in netCDF file.");
200	}
201	if (! tags_var->get(tags, ntags) )
202	{
203	esysUtils::free(tags);
204	throw DataException("Error - load:: unable to recover tags from netCDF file.");
205	}
206
207	DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
208	if (!(var = dataFile.get_var("data")))
209	{
210	esysUtils::free(tags);
211	throw DataException("Error - load:: unable to find data in netCDF file.");
212	}
213	if (! var->get(&(data[0]), dims) )
214	{
215	esysUtils::free(tags);
216	throw DataException("Error - load:: unable to recover data from netCDF file.");
217	}
218	out=Data(DataArrayView(data,shape,0),function_space);
219	for (int t=1; t<ntags; ++t) {
220	out.setTaggedValueFromCPP(tags[t],shape, data, t*len_data_point);
221	// out.setTaggedValueFromCPP(tags[t],DataArrayView(data,shape,t*len_data_point));
222	}
223	esysUtils::free(tags);
224	} else if (type == 2) {
225	/* expanded data */
226	if ( ! (ndims == rank + 2) )
227	throw DataException("Error - load:: illegal number of dimensions for exanded data in netCDF file.");
228	if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
229	throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
230	num_samples = num_samples_dim->size();
231	if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
232	throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
233	num_data_points_per_sample=num_data_points_per_sample_dim->size();
234	// check shape:
235	if ( ! (num_samples == function_space.getNumSamples() && num_data_points_per_sample == function_space.getNumDataPointsPerSample()) )
236	throw DataException("Error - load:: data sample layout of file does not match data layout of function space.");
237	// get ids
238	if (! ( ids_var = dataFile.get_var("id")) )
239	throw DataException("Error - load:: unable to find reference ids in netCDF file.");
240	const int* ids_p=function_space.borrowSampleReferenceIDs();
241	int ids_of_nc = (int )esysUtils::malloc(num_samples*sizeof(int));
242	if (! ids_var->get(ids_of_nc, (long) num_samples) )
243	{
244	esysUtils::free(ids_of_nc);
245	throw DataException("Error - load:: unable to recover ids from netCDF file.");
246	}
247	// check order:
248	int failed=-1, local_failed=-1, i;
249	#pragma omp parallel private(local_failed)
250	{
251	local_failed=-1;
252	#pragma omp for private(i) schedule(static)
253	for (i=0;i < num_samples; ++i) {
254	if (ids_of_nc[i]!=ids_p[i]) local_failed=i;
255	}
256	#pragma omp critical
257	if (local_failed>=0) failed = local_failed;
258	}
259	/* if (failed>=0)
260	{
261	esysUtils::free(ids_of_nc);
262	throw DataException("Error - load:: data ordering in netCDF file does not match ordering of FunctionSpace.");
263	} */
264	// get the data:
265	dims[rank]=num_data_points_per_sample;
266	dims[rank+1]=num_samples;
267	out=Data(0,shape,function_space,true);
268	if (!(var = dataFile.get_var("data")))
269	{
270	esysUtils::free(ids_of_nc);
271	throw DataException("Error - load:: unable to find data in netCDF file.");
272	}
273	if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) )
274	{
275	esysUtils::free(ids_of_nc);
276	throw DataException("Error - load:: unable to recover data from netCDF file.");
277	}
278	if (failed>=0) {
279	try {
280	std::cout << "Information - load: start reordering data from netCDF file " << fileName << std::endl;
281	out.borrowData()->reorderByReferenceIDs(ids_of_nc);
282	}
283	catch (std::exception&) {
284	esysUtils::free(ids_of_nc);
285	throw DataException("Error - load:: unable to reorder data in netCDF file.");
286	}
287	}
288	} else {
289	throw DataException("Error - load:: unknown escript data type in netCDF file.");
290	}
291	return out;
292	#else
293	throw DataException("Error - load:: is not compiled with netCFD. Please contact your insatllation manager.");
294	#endif
295	}
296
297	bool
298	loadConfigured()
299	{
300	#ifdef USE_NETCDF
301	return true;
302	#else
303	return false;
304	#endif
305	}
306
307	Data
308	convertToData(const boost::python::object& value,
309	const FunctionSpace& what)
310	{
311	// first we try to extract a Data object from value
312	extract<Data> value_data(value);
313	if (value_data.check()) {
314	Data extracted_data=value_data();
315	if (extracted_data.isEmpty()) {
316	return extracted_data;
317	} else {
318	return Data(extracted_data,what);
319	}
320	} else {
321	return Data(value,what);
322	}
323	}
324
325	} // end of namespace
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision