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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4616 - (show 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
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 #include <ripley/Brick.h>
17 #include <ripley/Rectangle.h>
18 #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 namespace ripley {
29
30 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 escript::Data readBinaryGrid(std::string filename, escript::FunctionSpace fs,
57 const object& pyShape, double fill, int byteOrder, int dataType,
58 const object& pyFirst, const object& pyNum, const object& pyMultiplier,
59 const object& pyReverse)
60 {
61 int dim=fs.getDim();
62 GridParameters params;
63
64 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
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 dom->readBinaryGrid(res, filename, params);
78 return res;
79 }
80
81 escript::Data readNcGrid(std::string filename, std::string varname,
82 escript::FunctionSpace fs, const object& pyShape, double fill,
83 const object& pyFirst, const object& pyNum, const object& pyMultiplier,
84 const object& pyReverse)
85 {
86 int dim=fs.getDim();
87 GridParameters params;
88
89 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
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 dom->readNcGrid(res, filename, varname, params);
101 return res;
102 }
103
104 // These wrappers are required to make the shared pointers work through the
105 // Python wrapper
106
107 // The double for n? is just to keep python happy when people need to deal with
108 // truediv
109 escript::Domain_ptr _brick(double _n0, double _n1, double _n2, const object& l0,
110 const object& l1, const object& l2, int d0, int d1, int d2)
111 {
112 int n0=static_cast<int>(_n0), n1=static_cast<int>(_n1), n2=static_cast<int>(_n2);
113 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 return escript::Domain_ptr(new Brick(n0,n1,n2, x0,y0,z0, x1,y1,z1, d0,d1,d2));
151 }
152
153 const int _q[]={0x61686969,0x746c4144,0x79616e43};
154 escript::Domain_ptr _rectangle(double _n0, double _n1, const object& l0,
155 const object& l1, int d0, int d1)
156 {
157 int n0=static_cast<int>(_n0), n1=static_cast<int>(_n1);
158 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 return escript::Domain_ptr(new Rectangle(n0,n1, x0,y0, x1,y1, d0,d1));
184 }
185 std::string _who(){int a[]={_q[0]^42,_q[1]^42,_q[2]^42,0};return (char*)&a[0];}
186
187 } // end of namespace ripley
188
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 // params are: bool show_user_defined, bool show_py_signatures, bool show_cpp_signatures
195 docstring_options docopt(true, true, false);
196 #endif
197
198 register_exception_translator<ripley::RipleyException>(&(esysUtils::esysExceptionTranslator));
199
200 scope().attr("__doc__") = "To use this module, please import esys.ripley";
201 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
208 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 "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 ":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 ":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
220 def("Rectangle", ripley::_rectangle, (arg("n0"),arg("n1"),arg("l0")=1.0,arg("l1")=1.0,arg("d0")=-1,arg("d1")=-1),
221 "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 ":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 ":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 def("_theculprit_", ripley::_who);
229
230 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
236 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
241 class_<ripley::RipleyDomain, bases<escript::AbstractContinuousDomain>, boost::noncopyable >
242 ("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 .def("writeBinaryGrid", &ripley::RipleyDomain::writeBinaryGrid)
247
248 .def("dump", &ripley::RipleyDomain::dump, args("filename"),
249 "Dumps the mesh to a file with the given name.")
250 .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 .def("getDescription", &ripley::RipleyDomain::getDescription,
257 ":return: a description for this domain\n:rtype: ``string``")
258 .def("getDim", &ripley::RipleyDomain::getDim, ":rtype: ``int``")
259 .def("getDataShape", &ripley::RipleyDomain::getDataShape, args("functionSpaceCode"),
260 ":return: a pair (dps, ns) where dps=the number of data points per sample, and ns=the number of samples\n:rtype: ``tuple``")
261 .def("getNumDataPointsGlobal", &ripley::RipleyDomain::getNumDataPointsGlobal,
262 ":return: the number of data points summed across all MPI processes\n"
263 ":rtype: ``int``")
264 .def("addPDEToSystem",&ripley::RipleyDomain::addPDEToSystem,
265 args("mat", "rhs", "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact"),
266 "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 ":param y_contact:\n:type y_contact: `Data`"
277 )
278 .def("addPDEToRHS",&ripley::RipleyDomain::addPDEToRHS,
279 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 .def("addPDEToTransportProblem",&ripley::RipleyDomain::addPDEToTransportProblem,
288 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 ":param y_contact:\n:type y_contact: `Data`"
302 )
303 .def("newOperator",&ripley::RipleyDomain::newSystemMatrix,
304 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 ":param type:\n:type type: ``int``"
311 )
312 .def("newTransportProblem",&ripley::RipleyDomain::newTransportProblem,
313 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 ":param type:\n:type type: ``int``"
319 )
320 .def("getSystemMatrixTypeId",&ripley::RipleyDomain::getSystemMatrixTypeId,
321 args("solver", "preconditioner", "package", "symmetry"),
322 ":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 ":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 ":param symmetry:\n:type symmetry: ``int``"
328 )
329 .def("getTransportTypeId",&ripley::RipleyDomain::getTransportTypeId,
330 args("solver", "preconditioner", "package", "symmetry"),
331 ":return: the identifier of the transport problem type to be used when a particular solver, preconditioner, package and symmetric matrix is used.\n"
332 ":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 ":param symmetry:\n:type symmetry: ``int``"
337 )
338 .def("getX",&ripley::RipleyDomain::getX, ":return: locations in the FEM nodes\n\n"
339 ":rtype: `Data`")
340 .def("getNormal",&ripley::RipleyDomain::getNormal,
341 ":return: boundary normals at the quadrature point on the face elements\n"
342 ":rtype: `Data`")
343 .def("getSize",&ripley::RipleyDomain::getSize,":return: the element size\n"
344 ":rtype: `Data`")
345 .def("setTagMap",&ripley::RipleyDomain::setTagMap,args("name","tag"),
346 "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 .def("getTag",&ripley::RipleyDomain::getTag,args("name"),":return: tag id for "
349 "``name``\n:rtype: ``string``")
350 .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
358 class_<ripley::Brick, bases<ripley::RipleyDomain> >("RipleyBrick", "", no_init);
359 class_<ripley::Rectangle, bases<ripley::RipleyDomain> >("RipleyRectangle", "", no_init)
360 .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 }
365

  ViewVC Help
Powered by ViewVC 1.1.26