/[escript]/trunk/ripley/src/ripleycpp.cpp
ViewVC logotype

Annotation of /trunk/ripley/src/ripleycpp.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4616 - (hide annotations)
Tue Jan 14 22:57:47 2014 UTC (5 years, 3 months ago) by caltinay
File size: 17239 byte(s)
- Fixed unit test failure by populating missing parts of the struct
- prepared for reversing read order from grids
- worked around a segfault caused by python version in jessie (2.7.6) in
  combination with gdal,pyproj and RTLD_GLOBAL in dlopen flags.


1 caltinay 3670
2 jfenwick 3981 /*****************************************************************************
3 caltinay 3670 *
4 jfenwick 4154 * Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 caltinay 3670 *
7     * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10     *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 caltinay 3670
16 caltinay 3691 #include <ripley/Brick.h>
17     #include <ripley/Rectangle.h>
18 caltinay 3670 #include <esysUtils/esysExceptionTranslator.h>
19    
20     #include <boost/python.hpp>
21     #include <boost/python/module.hpp>
22     #include <boost/python/def.hpp>
23     #include <boost/python/detail/defaults_gen.hpp>
24     #include <boost/version.hpp>
25    
26     using namespace boost::python;
27    
28 caltinay 3691 namespace ripley {
29 caltinay 3670
30 caltinay 4616 template<typename T>
31     std::vector<T> extractPyArray(const object& obj, const std::string& name,
32     int expectedLength=0)
33     {
34     std::vector<T> result;
35     if (extract<tuple>(obj).check() || extract<list>(obj).check()) {
36     if (expectedLength==0 || len(obj)==expectedLength) {
37     for (int i=0; i<len(obj); i++) {
38     result.push_back(extract<T>(obj[i]));
39     }
40     } else {
41     std::stringstream ssError;
42     ssError << "argument '" << name << "' has wrong length";
43     std::string error(ssError.str());
44     throw RipleyException(error.c_str());
45     }
46     } else {
47     std::stringstream ssError;
48     ssError << "argument '" << name << "' must be a tuple or list";
49     std::string error(ssError.str());
50     throw RipleyException(error.c_str());
51     }
52    
53     return result;
54     }
55    
56 caltinay 3971 escript::Data readBinaryGrid(std::string filename, escript::FunctionSpace fs,
57 caltinay 4616 const object& pyShape, double fill, int byteOrder, int dataType,
58 caltinay 4277 const object& pyFirst, const object& pyNum, const object& pyMultiplier,
59 caltinay 4616 const object& pyReverse)
60 caltinay 3971 {
61     int dim=fs.getDim();
62 caltinay 4615 GridParameters params;
63 caltinay 3971
64 caltinay 4616 params.first = extractPyArray<int>(pyFirst, "first", dim);
65     params.numValues = extractPyArray<int>(pyNum, "numValues", dim);
66     params.multiplier = extractPyArray<int>(pyMultiplier, "multiplier", dim);
67     params.reverse = extractPyArray<int>(pyReverse, "reverse", dim);
68     params.byteOrder = byteOrder;
69     params.dataType = dataType;
70     std::vector<int> shape(extractPyArray<int>(pyShape, "shape"));
71 caltinay 3971
72     const RipleyDomain* dom=dynamic_cast<const RipleyDomain*>(fs.getDomain().get());
73     if (!dom)
74     throw RipleyException("Function space must be on a ripley domain");
75    
76     escript::Data res(fill, shape, fs, true);
77 caltinay 4615 dom->readBinaryGrid(res, filename, params);
78 caltinay 3971 return res;
79     }
80    
81 caltinay 4013 escript::Data readNcGrid(std::string filename, std::string varname,
82 caltinay 4616 escript::FunctionSpace fs, const object& pyShape, double fill,
83     const object& pyFirst, const object& pyNum, const object& pyMultiplier,
84     const object& pyReverse)
85 caltinay 4013 {
86     int dim=fs.getDim();
87 caltinay 4615 GridParameters params;
88 caltinay 4013
89 caltinay 4616 params.first = extractPyArray<int>(pyFirst, "first", dim);
90     params.numValues = extractPyArray<int>(pyNum, "numValues", dim);
91     params.multiplier = extractPyArray<int>(pyMultiplier, "multiplier", dim);
92     params.reverse = extractPyArray<int>(pyReverse, "reverse", dim);
93     std::vector<int> shape(extractPyArray<int>(pyShape, "shape"));
94 caltinay 4013
95     const RipleyDomain* dom=dynamic_cast<const RipleyDomain*>(fs.getDomain().get());
96     if (!dom)
97     throw RipleyException("Function space must be on a ripley domain");
98    
99     escript::Data res(fill, shape, fs, true);
100 caltinay 4615 dom->readNcGrid(res, filename, varname, params);
101 caltinay 4013 return res;
102     }
103    
104 caltinay 3691 // These wrappers are required to make the shared pointers work through the
105     // Python wrapper
106 caltinay 3670
107 jfenwick 3892 // The double for n? is just to keep python happy when people need to deal with
108 caltinay 3971 // truediv
109 jfenwick 3892 escript::Domain_ptr _brick(double _n0, double _n1, double _n2, const object& l0,
110 caltinay 3781 const object& l1, const object& l2, int d0, int d1, int d2)
111 caltinay 3691 {
112 jfenwick 3892 int n0=static_cast<int>(_n0), n1=static_cast<int>(_n1), n2=static_cast<int>(_n2);
113 caltinay 3781 double x0=0., x1=1., y0=0., y1=1., z0=0., z1=1.;
114     if (extract<tuple>(l0).check()) {
115     tuple x=extract<tuple>(l0);
116     if (len(x)==2) {
117     x0=extract<double>(x[0]);
118     x1=extract<double>(x[1]);
119     } else
120     throw RipleyException("Argument l0 must be a float or 2-tuple");
121     } else if (extract<double>(l0).check()) {
122     x1=extract<double>(l0);
123     } else
124     throw RipleyException("Argument l0 must be a float or 2-tuple");
125    
126     if (extract<tuple>(l1).check()) {
127     tuple y=extract<tuple>(l1);
128     if (len(y)==2) {
129     y0=extract<double>(y[0]);
130     y1=extract<double>(y[1]);
131     } else
132     throw RipleyException("Argument l1 must be a float or 2-tuple");
133     } else if (extract<double>(l1).check()) {
134     y1=extract<double>(l1);
135     } else
136     throw RipleyException("Argument l1 must be a float or 2-tuple");
137    
138     if (extract<tuple>(l2).check()) {
139     tuple z=extract<tuple>(l2);
140     if (len(z)==2) {
141     z0=extract<double>(z[0]);
142     z1=extract<double>(z[1]);
143     } else
144     throw RipleyException("Argument l2 must be a float or 2-tuple");
145     } else if (extract<double>(l2).check()) {
146     z1=extract<double>(l2);
147     } else
148     throw RipleyException("Argument l2 must be a float or 2-tuple");
149    
150 caltinay 3791 return escript::Domain_ptr(new Brick(n0,n1,n2, x0,y0,z0, x1,y1,z1, d0,d1,d2));
151 caltinay 3691 }
152 caltinay 3670
153 caltinay 3791 const int _q[]={0x61686969,0x746c4144,0x79616e43};
154 jfenwick 3892 escript::Domain_ptr _rectangle(double _n0, double _n1, const object& l0,
155 caltinay 3781 const object& l1, int d0, int d1)
156 caltinay 3691 {
157 jfenwick 3892 int n0=static_cast<int>(_n0), n1=static_cast<int>(_n1);
158 caltinay 3781 double x0=0., x1=1., y0=0., y1=1.;
159     if (extract<tuple>(l0).check()) {
160     tuple x=extract<tuple>(l0);
161     if (len(x)==2) {
162     x0=extract<double>(x[0]);
163     x1=extract<double>(x[1]);
164     } else
165     throw RipleyException("Argument l0 must be a float or 2-tuple");
166     } else if (extract<double>(l0).check()) {
167     x1=extract<double>(l0);
168     } else
169     throw RipleyException("Argument l0 must be a float or 2-tuple");
170    
171     if (extract<tuple>(l1).check()) {
172     tuple y=extract<tuple>(l1);
173     if (len(y)==2) {
174     y0=extract<double>(y[0]);
175     y1=extract<double>(y[1]);
176     } else
177     throw RipleyException("Argument l1 must be a float or 2-tuple");
178     } else if (extract<double>(l1).check()) {
179     y1=extract<double>(l1);
180     } else
181     throw RipleyException("Argument l1 must be a float or 2-tuple");
182    
183 caltinay 3791 return escript::Domain_ptr(new Rectangle(n0,n1, x0,y0, x1,y1, d0,d1));
184 caltinay 3691 }
185 caltinay 3791 std::string _who(){int a[]={_q[0]^42,_q[1]^42,_q[2]^42,0};return (char*)&a[0];}
186 caltinay 3670
187 caltinay 4357 } // end of namespace ripley
188 caltinay 3670
189    
190     BOOST_PYTHON_MODULE(ripleycpp)
191     {
192     // This feature was added in boost v1.34
193     #if ((BOOST_VERSION/100)%1000 > 34) || (BOOST_VERSION/100000 >1)
194 caltinay 3691 // params are: bool show_user_defined, bool show_py_signatures, bool show_cpp_signatures
195     docstring_options docopt(true, true, false);
196 caltinay 3670 #endif
197    
198 caltinay 3691 register_exception_translator<ripley::RipleyException>(&(esysUtils::esysExceptionTranslator));
199 caltinay 3670
200 caltinay 4174 scope().attr("__doc__") = "To use this module, please import esys.ripley";
201 caltinay 4357 scope().attr("BYTEORDER_NATIVE") = (int)ripley::BYTEORDER_NATIVE;
202     scope().attr("BYTEORDER_LITTLE_ENDIAN") = (int)ripley::BYTEORDER_LITTLE_ENDIAN;
203     scope().attr("BYTEORDER_BIG_ENDIAN") = (int)ripley::BYTEORDER_BIG_ENDIAN;
204     scope().attr("DATATYPE_INT32") = (int)ripley::DATATYPE_INT32;
205     scope().attr("DATATYPE_FLOAT32") = (int)ripley::DATATYPE_FLOAT32;
206     scope().attr("DATATYPE_FLOAT64") = (int)ripley::DATATYPE_FLOAT64;
207 caltinay 4174
208 caltinay 3943 def("Brick", ripley::_brick, (arg("n0"),arg("n1"),arg("n2"),arg("l0")=1.0,arg("l1")=1.0,arg("l2")=1.0,arg("d0")=-1,arg("d1")=-1,arg("d2")=-1),
209 caltinay 3691 "Creates a hexagonal mesh with n0 x n1 x n2 elements over the brick [0,l0] x [0,l1] x [0,l2].\n\n"
210     ":param n0: number of elements in direction 0\n:type n0: ``int``\n"
211     ":param n1: number of elements in direction 1\n:type n1: ``int``\n"
212     ":param n2: number of elements in direction 2\n:type n2: ``int``\n"
213 caltinay 3781 ":param l0: length of side 0 or coordinate range of side 0\n:type l0: ``float`` or ``tuple``\n"
214     ":param l1: length of side 1 or coordinate range of side 1\n:type l1: ``float`` or ``tuple``\n"
215     ":param l2: length of side 2 or coordinate range of side 2\n:type l2: ``float`` or ``tuple``\n"
216 caltinay 3691 ":param d0: number of subdivisions in direction 0\n:type d0: ``int``\n"
217     ":param d1: number of subdivisions in direction 1\n:type d1: ``int``\n"
218     ":param d2: number of subdivisions in direction 2\n:type d2: ``int``");
219 caltinay 3670
220 caltinay 3943 def("Rectangle", ripley::_rectangle, (arg("n0"),arg("n1"),arg("l0")=1.0,arg("l1")=1.0,arg("d0")=-1,arg("d1")=-1),
221 caltinay 3691 "Creates a rectangular mesh with n0 x n1 elements over the rectangle [0,l0] x [0,l1].\n\n"
222     ":param n0: number of elements in direction 0\n:type n0: ``int``\n"
223     ":param n1: number of elements in direction 1\n:type n1: ``int``\n"
224 caltinay 3781 ":param l0: length of side 0 or coordinate range of side 0\n:type l0: ``float`` or ``tuple``\n"
225     ":param l1: length of side 1 or coordinate range of side 1\n:type l1: ``float`` or ``tuple``\n"
226 caltinay 3691 ":param d0: number of subdivisions in direction 0\n:type d0: ``int``\n"
227     ":param d1: number of subdivisions in direction 1\n:type d1: ``int``");
228 caltinay 3971 def("_theculprit_", ripley::_who);
229 caltinay 3670
230 caltinay 4616 def("_readBinaryGrid", &ripley::readBinaryGrid, (arg("filename"),
231     arg("functionspace"), arg("shape"), arg("fill")=0.,
232     arg("byteOrder"), arg("dataType"), arg("first"),
233     arg("numValues"), arg("multiplier"), arg("reverse")),
234     "Reads a binary Grid");
235 caltinay 3971
236 caltinay 4616 def("_readNcGrid", &ripley::readNcGrid, (arg("filename"), arg("varname"),
237     arg("functionspace"), arg("shape"), arg("fill"), arg("first"),
238     arg("numValues"), arg("multiplier"), arg("reverse")),
239     "Reads a grid from a netCDF file");
240 caltinay 4013
241 jfenwick 4242 class_<ripley::RipleyDomain, bases<escript::AbstractContinuousDomain>, boost::noncopyable >
242 caltinay 3691 ("RipleyDomain", "", no_init)
243     .def("print_mesh_info", &ripley::RipleyDomain::Print_Mesh_Info, (arg("full")=false),
244     "Prints out a summary about the mesh.\n"
245     ":param full: whether to output additional data\n:type full: ``bool``")
246 caltinay 4357 .def("writeBinaryGrid", &ripley::RipleyDomain::writeBinaryGrid)
247    
248 caltinay 3691 .def("dump", &ripley::RipleyDomain::dump, args("filename"),
249     "Dumps the mesh to a file with the given name.")
250 caltinay 4340 .def("getGridParameters", &ripley::RipleyDomain::getGridParameters,
251     "Returns the tuple (origin, spacing, elements) where the entries are tuples:\n"
252     "``origin``=the coordinates of the domain's global origin,\n"
253     "``spacing``=the element size (=node spacing) of the domain,\n"
254     "``elements``=the global number of elements in all dimensions\n\n"
255     ":rtype: ``tuple``")
256 caltinay 3691 .def("getDescription", &ripley::RipleyDomain::getDescription,
257 caltinay 3670 ":return: a description for this domain\n:rtype: ``string``")
258 caltinay 3691 .def("getDim", &ripley::RipleyDomain::getDim, ":rtype: ``int``")
259     .def("getDataShape", &ripley::RipleyDomain::getDataShape, args("functionSpaceCode"),
260 caltinay 3670 ":return: a pair (dps, ns) where dps=the number of data points per sample, and ns=the number of samples\n:rtype: ``tuple``")
261 caltinay 3691 .def("getNumDataPointsGlobal", &ripley::RipleyDomain::getNumDataPointsGlobal,
262 caltinay 3670 ":return: the number of data points summed across all MPI processes\n"
263     ":rtype: ``int``")
264 caltinay 3691 .def("addPDEToSystem",&ripley::RipleyDomain::addPDEToSystem,
265     args("mat", "rhs", "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact"),
266 caltinay 3670 "adds a PDE onto the stiffness matrix mat and a rhs\n\n"
267     ":param mat:\n:type mat: `OperatorAdapter`\n:param rhs:\n:type rhs: `Data`\n"
268     ":param A:\n:type A: `Data`\n"
269     ":param B:\n:type B: `Data`\n"
270     ":param C:\n:type C: `Data`\n"
271     ":param D:\n:type D: `Data`\n"
272     ":param X:\n:type X: `Data`\n"
273     ":param Y:\n:type Y: `Data`\n"
274     ":param d:\n:type d: `Data`\n"
275     ":param d_contact:\n:type d_contact: `Data`\n"
276 caltinay 3791 ":param y_contact:\n:type y_contact: `Data`"
277 caltinay 3670 )
278 caltinay 3691 .def("addPDEToRHS",&ripley::RipleyDomain::addPDEToRHS,
279 caltinay 3670 args("rhs", "X", "Y", "y", "y_contact"),
280     "adds a PDE onto the stiffness matrix mat and a rhs\n\n"
281     ":param rhs:\n:type rhs: `Data`\n"
282     ":param X:\n:type X: `Data`\n"
283     ":param Y:\n:type Y: `Data`\n"
284     ":param y:\n:type y: `Data`\n"
285     ":param y_contact:\n:type y_contact: `Data`"
286     )
287 caltinay 3691 .def("addPDEToTransportProblem",&ripley::RipleyDomain::addPDEToTransportProblem,
288 caltinay 3670 args( "tp", "source", "M", "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact"),
289     ":param tp:\n:type tp: `TransportProblemAdapter`\n"
290     ":param source:\n:type source: `Data`\n"
291     ":param M:\n:type M: `Data`\n"
292     ":param A:\n:type A: `Data`\n"
293     ":param B:\n:type B: `Data`\n"
294     ":param C:\n:type C: `Data`\n"
295     ":param D:\n:type D: `Data`\n"
296     ":param X:\n:type X: `Data`\n"
297     ":param Y:\n:type Y: `Data`\n"
298     ":param d:\n:type d: `Data`\n"
299     ":param y:\n:type y: `Data`\n"
300     ":param d_contact:\n:type d_contact: `Data`\n"
301 caltinay 3791 ":param y_contact:\n:type y_contact: `Data`"
302 caltinay 3670 )
303 caltinay 3691 .def("newOperator",&ripley::RipleyDomain::newSystemMatrix,
304 caltinay 3670 args("row_blocksize", "row_functionspace", "column_blocksize", "column_functionspace", "type"),
305     "creates a SystemMatrixAdapter stiffness matrix and initializes it with zeros\n\n"
306     ":param row_blocksize:\n:type row_blocksize: ``int``\n"
307     ":param row_functionspace:\n:type row_functionspace: `FunctionSpace`\n"
308     ":param column_blocksize:\n:type column_blocksize: ``int``\n"
309     ":param column_functionspace:\n:type column_functionspace: `FunctionSpace`\n"
310 caltinay 3791 ":param type:\n:type type: ``int``"
311 caltinay 3670 )
312 caltinay 3691 .def("newTransportProblem",&ripley::RipleyDomain::newTransportProblem,
313 caltinay 3670 args("theta", "blocksize", "functionspace", "type"),
314     "creates a TransportProblemAdapter\n\n"
315     ":param theta:\n:type theta: ``float``\n"
316     ":param blocksize:\n:type blocksize: ``int``\n"
317     ":param functionspace:\n:type functionspace: `FunctionSpace`\n"
318 caltinay 3791 ":param type:\n:type type: ``int``"
319 caltinay 3670 )
320 caltinay 3691 .def("getSystemMatrixTypeId",&ripley::RipleyDomain::getSystemMatrixTypeId,
321 caltinay 3670 args("solver", "preconditioner", "package", "symmetry"),
322 caltinay 3697 ":return: the identifier of the matrix type to be used for the global stiffness matrix when a particular solver, package, preconditioner, and symmetric matrix is used.\n"
323 caltinay 3670 ":rtype: ``int``\n"
324     ":param solver:\n:type solver: ``int``\n"
325     ":param preconditioner:\n:type preconditioner: ``int``\n"
326     ":param package:\n:type package: ``int``\n"
327 caltinay 3791 ":param symmetry:\n:type symmetry: ``int``"
328 caltinay 3670 )
329 caltinay 3691 .def("getTransportTypeId",&ripley::RipleyDomain::getTransportTypeId,
330 caltinay 3670 args("solver", "preconditioner", "package", "symmetry"),
331 caltinay 3697 ":return: the identifier of the transport problem type to be used when a particular solver, preconditioner, package and symmetric matrix is used.\n"
332 caltinay 3670 ":rtype: ``int``\n"
333     ":param solver:\n:type solver: ``int``\n"
334     ":param preconditioner:\n:type preconditioner: ``int``\n"
335     ":param package:\n:type package: ``int``\n"
336 caltinay 3791 ":param symmetry:\n:type symmetry: ``int``"
337 caltinay 3670 )
338 caltinay 3691 .def("getX",&ripley::RipleyDomain::getX, ":return: locations in the FEM nodes\n\n"
339 caltinay 3670 ":rtype: `Data`")
340 caltinay 3691 .def("getNormal",&ripley::RipleyDomain::getNormal,
341 caltinay 3670 ":return: boundary normals at the quadrature point on the face elements\n"
342     ":rtype: `Data`")
343 caltinay 3691 .def("getSize",&ripley::RipleyDomain::getSize,":return: the element size\n"
344 caltinay 3670 ":rtype: `Data`")
345 caltinay 3691 .def("setTagMap",&ripley::RipleyDomain::setTagMap,args("name","tag"),
346 caltinay 3670 "Give a tag number a name.\n\n:param name: Name for the tag\n:type name: ``string``\n"
347     ":param tag: numeric id\n:type tag: ``int``\n:note: Tag names must be unique within a domain")
348 caltinay 3691 .def("getTag",&ripley::RipleyDomain::getTag,args("name"),":return: tag id for "
349 caltinay 3670 "``name``\n:rtype: ``string``")
350 caltinay 3691 .def("isValidTagName",&ripley::RipleyDomain::isValidTagName,args("name"),
351     ":return: True if ``name`` corresponds to a tag, otherwise False\n:rtype: ``bool``")
352     .def("showTagNames",&ripley::RipleyDomain::showTagNames,":return: A space separated list of tag names\n:rtype: ``string``")
353     .def("getMPISize",&ripley::RipleyDomain::getMPISize,":return: the number of processes used for this `Domain`\n:rtype: ``int``")
354     .def("getMPIRank",&ripley::RipleyDomain::getMPIRank,":return: the rank of this process\n:rtype: ``int``")
355     .def("MPIBarrier",&ripley::RipleyDomain::MPIBarrier,"Wait until all processes have reached this point")
356     .def("onMasterProcessor",&ripley::RipleyDomain::onMasterProcessor,":return: True if this code is executing on the master process\n:rtype: `bool`");
357 caltinay 3670
358 caltinay 3691 class_<ripley::Brick, bases<ripley::RipleyDomain> >("RipleyBrick", "", no_init);
359 jfenwick 4526 class_<ripley::Rectangle, bases<ripley::RipleyDomain> >("RipleyRectangle", "", no_init)
360 jfenwick 4586 .def("randomFill", &ripley::Rectangle::randomFill,":return: random data\n:rtype: `Data`\n:param seed: pass zero to use system generated seed\n:type seed: `int`\n"
361     ":param details: more info about the type of randomness\nCurrently, the only acceptable value for this tuple is ('gaussian', r, s) where r is the radius of the"
362     "guassian blur and s is the sigma parameter."
363     );
364 caltinay 3670 }
365    

  ViewVC Help
Powered by ViewVC 1.1.26