/[escript]/trunk/speckley/src/speckleycpp.cpp
ViewVC logotype

Contents of /trunk/speckley/src/speckleycpp.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6141 - (show annotations)
Wed Apr 6 03:51:30 2016 UTC (2 years, 8 months ago) by caltinay
File size: 22200 byte(s)
more namespacing of defines.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include <speckley/AbstractAssembler.h>
18 #include <speckley/Brick.h>
19 #include <speckley/Rectangle.h>
20
21 #include <escript/ExceptionTranslators.h>
22 #include <escript/SubWorld.h>
23
24 #include <boost/python.hpp>
25 #include <boost/python/module.hpp>
26 #include <boost/python/def.hpp>
27 #include <boost/python/detail/defaults_gen.hpp>
28 #include <boost/version.hpp>
29
30 using namespace boost::python;
31
32 namespace speckley {
33
34 template<typename T>
35 std::vector<T> extractPyArray(const object& obj, const std::string& name,
36 int expectedLength=0)
37 {
38 std::vector<T> result;
39 if (extract<tuple>(obj).check() || extract<list>(obj).check()) {
40 if (expectedLength==0 || len(obj)==expectedLength) {
41 for (int i=0; i<len(obj); i++) {
42 result.push_back(extract<T>(obj[i]));
43 }
44 } else {
45 std::stringstream ssError;
46 ssError << "argument '" << name << "' has wrong length";
47 std::string error(ssError.str());
48 throw SpeckleyException(error.c_str());
49 }
50 } else {
51 std::stringstream ssError;
52 ssError << "argument '" << name << "' must be a tuple or list";
53 std::string error(ssError.str());
54 throw SpeckleyException(error.c_str());
55 }
56
57 return result;
58 }
59
60 escript::Data readBinaryGrid(std::string filename, escript::FunctionSpace fs,
61 const object& pyShape, double fill, int byteOrder, int dataType,
62 const object& pyFirst, const object& pyNum, const object& pyMultiplier,
63 const object& pyReverse)
64 {
65 int dim=fs.getDim();
66 ReaderParameters params;
67
68 params.first = extractPyArray<index_t>(pyFirst, "first", dim);
69 params.numValues = extractPyArray<dim_t>(pyNum, "numValues", dim);
70 params.multiplier = extractPyArray<int>(pyMultiplier, "multiplier", dim);
71 params.reverse = extractPyArray<int>(pyReverse, "reverse", dim);
72 params.byteOrder = byteOrder;
73 params.dataType = dataType;
74 std::vector<int> shape(extractPyArray<int>(pyShape, "shape"));
75
76 const SpeckleyDomain* dom=dynamic_cast<const SpeckleyDomain*>(fs.getDomain().get());
77 if (!dom)
78 throw SpeckleyException("Function space must be on a speckley domain");
79
80 escript::Data res(fill, shape, fs, true);
81 dom->readBinaryGrid(res, filename, params);
82 return res;
83 }
84
85 escript::Data readBinaryGridFromZipped(std::string filename, escript::FunctionSpace fs,
86 const object& pyShape, double fill, int byteOrder, int dataType,
87 const object& pyFirst, const object& pyNum, const object& pyMultiplier,
88 const object& pyReverse)
89 {
90 #ifdef ESYS_HAVE_BOOST_IO
91 int dim=fs.getDim();
92 ReaderParameters params;
93
94 params.first = extractPyArray<index_t>(pyFirst, "first", dim);
95 params.numValues = extractPyArray<dim_t>(pyNum, "numValues", dim);
96 params.multiplier = extractPyArray<int>(pyMultiplier, "multiplier", dim);
97 params.reverse = extractPyArray<int>(pyReverse, "reverse", dim);
98 params.byteOrder = byteOrder;
99 params.dataType = dataType;
100 std::vector<int> shape(extractPyArray<int>(pyShape, "shape"));
101
102 const SpeckleyDomain* dom=dynamic_cast<const SpeckleyDomain*>(fs.getDomain().get());
103 if (!dom)
104 throw SpeckleyException("Function space must be on a speckley domain");
105
106 escript::Data res(fill, shape, fs, true);
107 dom->readBinaryGridFromZipped(res, filename, params);
108 return res;
109 #else
110 throw SpeckleyException("Speckley was not built with zip support!");
111 #endif
112 }
113
114 escript::Data readNcGrid(std::string filename, std::string varname,
115 escript::FunctionSpace fs, const object& pyShape, double fill,
116 const object& pyFirst, const object& pyNum, const object& pyMultiplier,
117 const object& pyReverse)
118 {
119 int dim=fs.getDim();
120 ReaderParameters params;
121
122 params.first = extractPyArray<index_t>(pyFirst, "first", dim);
123 params.numValues = extractPyArray<dim_t>(pyNum, "numValues", dim);
124 params.multiplier = extractPyArray<int>(pyMultiplier, "multiplier", dim);
125 params.reverse = extractPyArray<int>(pyReverse, "reverse", dim);
126 std::vector<int> shape(extractPyArray<int>(pyShape, "shape"));
127
128 const SpeckleyDomain* dom=dynamic_cast<const SpeckleyDomain*>(fs.getDomain().get());
129 if (!dom)
130 throw SpeckleyException("Function space must be on a speckley domain");
131
132 escript::Data res(fill, shape, fs, true);
133 dom->readNcGrid(res, filename, varname, params);
134 return res;
135 }
136
137 // These wrappers are required to make the shared pointers work through the
138 // Python wrapper
139
140 // The double for n? is just to keep python happy when people need to deal with
141 // truediv
142 escript::Domain_ptr _brick(int order, double _n0, double _n1, double _n2,
143 const object& l0, const object& l1, const object& l2, int d0,
144 int d1, int d2, const object& objpoints,
145 const object& objtags, escript::SubWorld_ptr world)
146 {
147 dim_t n0=static_cast<dim_t>(_n0), n1=static_cast<dim_t>(_n1), n2=static_cast<dim_t>(_n2);
148 double x0=0., x1=1., y0=0., y1=1., z0=0., z1=1.;
149 if (order < 2 || order > 10)
150 throw SpeckleyException("Order must be in the range 2 to 10");
151 if (extract<tuple>(l0).check()) {
152 tuple x=extract<tuple>(l0);
153 if (len(x)==2) {
154 x0=extract<double>(x[0]);
155 x1=extract<double>(x[1]);
156 } else
157 throw SpeckleyException("Argument l0 must be a float or 2-tuple");
158 } else if (extract<double>(l0).check()) {
159 x1=extract<double>(l0);
160 } else
161 throw SpeckleyException("Argument l0 must be a float or 2-tuple");
162
163 if (extract<tuple>(l1).check()) {
164 tuple y=extract<tuple>(l1);
165 if (len(y)==2) {
166 y0=extract<double>(y[0]);
167 y1=extract<double>(y[1]);
168 } else
169 throw SpeckleyException("Argument l1 must be a float or 2-tuple");
170 } else if (extract<double>(l1).check()) {
171 y1=extract<double>(l1);
172 } else
173 throw SpeckleyException("Argument l1 must be a float or 2-tuple");
174
175 if (extract<tuple>(l2).check()) {
176 tuple z=extract<tuple>(l2);
177 if (len(z)==2) {
178 z0=extract<double>(z[0]);
179 z1=extract<double>(z[1]);
180 } else
181 throw SpeckleyException("Argument l2 must be a float or 2-tuple");
182 } else if (extract<double>(l2).check()) {
183 z1=extract<double>(l2);
184 } else
185 throw SpeckleyException("Argument l2 must be a float or 2-tuple");
186 boost::python::list pypoints=extract<boost::python::list>(objpoints);
187 boost::python::list pytags=extract<boost::python::list>(objtags);
188 int numpts=extract<int>(pypoints.attr("__len__")());
189 int numtags=extract<int>(pytags.attr("__len__")());
190 std::vector<double> points;
191 std::vector<int> tags;
192 tags.resize(numtags, -1);
193 for (int i=0;i<numpts;++i) {
194 tuple temp = extract<tuple>(pypoints[i]);
195 int l=extract<int>(temp.attr("__len__")());
196 if (l != 3)
197 throw SpeckleyException("Number of coordinates for each dirac point must match dimensions.");
198 for (int k=0; k<l; ++k) {
199 points.push_back(extract<double>(temp[k]));
200 }
201 }
202 TagMap tagstonames;
203 int curmax=40;
204 // but which order to assign tags to names?????
205 for (int i=0;i<numtags;++i) {
206 extract<int> ex_int(pytags[i]);
207 extract<std::string> ex_str(pytags[i]);
208 if (ex_int.check()) {
209 tags[i]=ex_int();
210 if (tags[i]>= curmax) {
211 curmax=tags[i]+1;
212 }
213 } else if (ex_str.check()) {
214 std::string s=ex_str();
215 TagMap::iterator it=tagstonames.find(s);
216 if (it!=tagstonames.end()) {
217 // we have the tag already so look it up
218 tags[i]=it->second;
219 } else {
220 tagstonames[s]=curmax;
221 tags[i]=curmax;
222 curmax++;
223 }
224 } else {
225 throw SpeckleyException("Error - Unable to extract tag value.");
226 }
227 }
228 if (numtags != numpts)
229 throw SpeckleyException("Number of tags does not match number of points.");
230 return escript::Domain_ptr(new Brick(order, n0,n1,n2, x0,y0,z0, x1,y1,z1, d0,d1,d2,
231 points, tags, tagstonames, world));
232 }
233
234 escript::Domain_ptr _rectangle(int order, double _n0, double _n1,
235 const object& l0, const object& l1, int d0, int d1,
236 const object& objpoints, const object& objtags,
237 escript::SubWorld_ptr world
238 )
239 {
240 if (order < 2 || order > 10)
241 throw SpeckleyException("Order must be in the range 2 to 10");
242 dim_t n0=static_cast<dim_t>(_n0), n1=static_cast<dim_t>(_n1);
243 double x0=0., x1=1., y0=0., y1=1.;
244 if (extract<tuple>(l0).check()) {
245 tuple x=extract<tuple>(l0);
246 if (len(x)==2) {
247 x0=extract<double>(x[0]);
248 x1=extract<double>(x[1]);
249 } else
250 throw SpeckleyException("Argument l0 must be a float or 2-tuple");
251 } else if (extract<double>(l0).check()) {
252 x1=extract<double>(l0);
253 } else
254 throw SpeckleyException("Argument l0 must be a float or 2-tuple");
255
256 if (extract<tuple>(l1).check()) {
257 tuple y=extract<tuple>(l1);
258 if (len(y)==2) {
259 y0=extract<double>(y[0]);
260 y1=extract<double>(y[1]);
261 } else
262 throw SpeckleyException("Argument l1 must be a float or 2-tuple");
263 } else if (extract<double>(l1).check()) {
264 y1=extract<double>(l1);
265 } else
266 throw SpeckleyException("Argument l1 must be a float or 2-tuple");
267 boost::python::list pypoints=extract<boost::python::list>(objpoints);
268 boost::python::list pytags=extract<boost::python::list>(objtags);
269 int numpts=extract<int>(pypoints.attr("__len__")());
270 int numtags=extract<int>(pytags.attr("__len__")());
271 std::vector<double> points;
272 std::vector<int> tags;
273 tags.resize(numtags, -1);
274 for (int i=0;i<numpts;++i) {
275 tuple temp = extract<tuple>(pypoints[i]);
276 int l=extract<int>(temp.attr("__len__")());
277 if (l != 2)
278 throw SpeckleyException("Number of coordinates for each dirac point must match dimensions.");
279 for (int k=0;k<l;++k) {
280 points.push_back(extract<double>(temp[k]));
281 }
282 }
283 TagMap tagstonames;
284 int curmax=40;
285 // but which order to assign tags to names?????
286 for (int i=0;i<numtags;++i) {
287 extract<int> ex_int(pytags[i]);
288 extract<std::string> ex_str(pytags[i]);
289 if (ex_int.check()) {
290 tags[i]=ex_int();
291 if (tags[i]>= curmax) {
292 curmax=tags[i]+1;
293 }
294 } else if (ex_str.check()) {
295 std::string s=ex_str();
296 TagMap::iterator it=tagstonames.find(s);
297 if (it!=tagstonames.end()) {
298 // we have the tag already so look it up
299 tags[i]=it->second;
300 } else {
301 tagstonames[s]=curmax;
302 tags[i]=curmax;
303 curmax++;
304 }
305 } else {
306 throw SpeckleyException("Error - Unable to extract tag value.");
307 }
308 }
309 if (numtags != numpts)
310 throw SpeckleyException("Number of tags does not match number of points.");
311 return escript::Domain_ptr(new Rectangle(order, n0,n1, x0,y0, x1,y1, d0,d1,
312 points, tags, tagstonames, world));
313 }
314
315 } // end of namespace speckley
316
317
318 BOOST_PYTHON_MODULE(speckleycpp)
319 {
320 // This feature was added in boost v1.34
321 #if ((BOOST_VERSION/100)%1000 > 34) || (BOOST_VERSION/100000 >1)
322 // params are: bool show_user_defined, bool show_py_signatures, bool show_cpp_signatures
323 docstring_options docopt(true, true, false);
324 #endif
325
326 // register escript's default translators
327 REGISTER_ESCRIPT_EXCEPTION_TRANSLATORS;
328 register_exception_translator<speckley::SpeckleyException>(&escript::RuntimeErrorTranslator);
329
330 scope().attr("__doc__") = "To use this module, please import esys.speckley";
331 scope().attr("BYTEORDER_NATIVE") = (int)speckley::BYTEORDER_NATIVE;
332 scope().attr("BYTEORDER_LITTLE_ENDIAN") = (int)speckley::BYTEORDER_LITTLE_ENDIAN;
333 scope().attr("BYTEORDER_BIG_ENDIAN") = (int)speckley::BYTEORDER_BIG_ENDIAN;
334 scope().attr("DATATYPE_INT32") = (int)speckley::DATATYPE_INT32;
335 scope().attr("DATATYPE_FLOAT32") = (int)speckley::DATATYPE_FLOAT32;
336 scope().attr("DATATYPE_FLOAT64") = (int)speckley::DATATYPE_FLOAT64;
337
338 def("Brick", speckley::_brick, (arg("order"),arg("n0"),arg("n1"),arg("n2"),arg("l0")=1.0,arg("l1")=1.0,arg("l2")=1.0,
339 arg("d0")=-1,arg("d1")=-1,arg("d2")=-1,arg("diracPoints")=list(),arg("diracTags")=list(), arg("escriptworld")=escript::SubWorld_ptr()),
340 "Creates a hexagonal mesh with n0 x n1 x n2 elements over the brick [0,l0] x [0,l1] x [0,l2].\n\n"
341 ":param order: the number of quadrature points to have in each dimension\n:type n0: ``int``\n"
342 ":param n0: number of elements in direction 0\n:type n0: ``int``\n"
343 ":param n1: number of elements in direction 1\n:type n1: ``int``\n"
344 ":param n2: number of elements in direction 2\n:type n2: ``int``\n"
345 ":param l0: length of side 0 or coordinate range of side 0\n:type l0: ``float`` or ``tuple``\n"
346 ":param l1: length of side 1 or coordinate range of side 1\n:type l1: ``float`` or ``tuple``\n"
347 ":param l2: length of side 2 or coordinate range of side 2\n:type l2: ``float`` or ``tuple``\n"
348 ":param d0: number of subdivisions in direction 0\n:type d0: ``int``\n"
349 ":param d1: number of subdivisions in direction 1\n:type d1: ``int``\n"
350 ":param d2: number of subdivisions in direction 2\n:type d2: ``int``");
351
352 def("Rectangle", speckley::_rectangle, (arg("order"),arg("n0"),arg("n1"),arg("l0")=1.0,arg("l1")=1.0,arg("d0")=-1,arg("d1")=-1,arg("diracPoints")=list(),arg("diracTags")=list(), arg("escriptworld")=escript::SubWorld_ptr()),
353 "Creates a rectangular mesh with n0 x n1 elements over the rectangle [0,l0] x [0,l1].\n\n"
354 ":param n0: number of elements in direction 0\n:type n0: ``int``\n"
355 ":param n1: number of elements in direction 1\n:type n1: ``int``\n"
356 ":param l0: length of side 0 or coordinate range of side 0\n:type l0: ``float`` or ``tuple``\n"
357 ":param l1: length of side 1 or coordinate range of side 1\n:type l1: ``float`` or ``tuple``\n"
358 ":param d0: number of subdivisions in direction 0\n:type d0: ``int``\n"
359 ":param d1: number of subdivisions in direction 1\n:type d1: ``int``");
360
361 def("readBinaryGrid", &speckley::readBinaryGrid, (arg("filename"),
362 arg("functionspace"), arg("shape"), arg("fill")=0.,
363 arg("byteOrder"), arg("dataType"), arg("first"),
364 arg("numValues"), arg("multiplier"), arg("reverse")),
365 "Reads a binary Grid");
366 def("_readBinaryGridFromZipped", &speckley::readBinaryGridFromZipped, (arg("filename"),
367 arg("functionspace"), arg("shape"), arg("fill")=0.,
368 arg("byteOrder"), arg("dataType"), arg("first"),
369 arg("numValues"), arg("multiplier"), arg("reverse")),
370 "Reads a binary Grid");
371 def("_readNcGrid", &speckley::readNcGrid, (arg("filename"), arg("varname"),
372 arg("functionspace"), arg("shape"), arg("fill"), arg("first"),
373 arg("numValues"), arg("multiplier"), arg("reverse")),
374 "Reads a grid from a netCDF file");
375
376 class_<speckley::SpeckleyDomain, bases<escript::AbstractContinuousDomain>, boost::noncopyable >
377 ("SpeckleyDomain", "", no_init)
378 .def("print_mesh_info", &speckley::SpeckleyDomain::Print_Mesh_Info, (arg("full")=false),
379 "Prints out a summary about the mesh.\n"
380 ":param full: whether to output additional data\n:type full: ``bool``")
381 .def("writeBinaryGrid", &speckley::SpeckleyDomain::writeBinaryGrid)
382
383 .def("dump", &speckley::SpeckleyDomain::dump, args("filename"),
384 "Dumps the mesh to a file with the given name.")
385 .def("getGridParameters", &speckley::SpeckleyDomain::getGridParameters,
386 "Returns the tuple (origin, spacing, elements) where the entries are tuples:\n"
387 "``origin`` the coordinates of the domain's global origin,\n"
388 "``spacing`` the element size (=node spacing) of the domain,\n"
389 "``elements`` the global number of elements in all dimensions\n\n"
390 ":rtype: ``tuple``")
391 .def("getDescription", &speckley::SpeckleyDomain::getDescription,
392 ":return: a description for this domain\n:rtype: ``string``")
393 .def("getDim", &speckley::SpeckleyDomain::getDim, ":rtype: ``int``")
394 .def("getDataShape", &speckley::SpeckleyDomain::getDataShape, args("functionSpaceCode"),
395 ":return: a pair (dps, ns) where dps is the number of data points"
396 " per sample, and ns isthe number of samples\n"
397 ":rtype: ``tuple``")
398 .def("getNumDataPointsGlobal", &speckley::SpeckleyDomain::getNumDataPointsGlobal,
399 ":return: the number of data points summed across all MPI processes\n"
400 ":rtype: ``int``")
401 .def("addToSystem",&speckley::SpeckleyDomain::addToSystemFromPython,
402 args("mat", "rhs", "data"),
403 "adds a PDE to the system, results depend on domain\n\n"
404 ":param mat:\n:type mat: `OperatorAdapter`\n"
405 ":param rhs:\n:type rhs: `Data`\n"
406 ":param data:\n:type data: `list`")
407 .def("addToRHS",&speckley::SpeckleyDomain::addToRHSFromPython,
408 args("rhs", "data"),
409 "adds a PDE onto the stiffness matrix mat and a rhs, "
410 "results depends on domain\n\n"
411 ":param rhs:\n:type rhs: `Data`\n"
412 ":param data:\n:type data: `list`")
413 .def("getOrder",&speckley::SpeckleyDomain::getOrder,
414 ":return: the order of the domain\n"
415 ":rtype: ``int``")
416 .def("createAssembler", &speckley::SpeckleyDomain::createAssemblerFromPython,
417 args("typename", "options"),
418 "request from the domain an assembler of the specified type, if "
419 "supported, using the supplied options (if provided)"
420 ":param typename:\n:type typename: `string`\n"
421 ":param options:\n:type options: `list`\n")
422 .def("newOperator",&speckley::SpeckleyDomain::newSystemMatrix,
423 args("row_blocksize", "row_functionspace", "column_blocksize", "column_functionspace", "type"),
424 "creates a SystemMatrixAdapter stiffness matrix and initializes it with zeros\n\n"
425 ":param row_blocksize:\n:type row_blocksize: ``int``\n"
426 ":param row_functionspace:\n:type row_functionspace: `FunctionSpace`\n"
427 ":param column_blocksize:\n:type column_blocksize: ``int``\n"
428 ":param column_functionspace:\n:type column_functionspace: `FunctionSpace`\n"
429 ":param type:\n:type type: ``int``"
430 )
431 .def("getSystemMatrixTypeId",&speckley::SpeckleyDomain::getSystemMatrixTypeId,
432 args("options"),
433 ":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"
434 ":rtype: ``int``\n"
435 ":param solver:\n:type solver: ``int``\n"
436 ":param preconditioner:\n:type preconditioner: ``int``\n"
437 ":param package:\n:type package: ``int``\n"
438 ":param symmetry:\n:type symmetry: ``int``"
439 )
440 .def("writeBinaryGrid", &speckley::SpeckleyDomain::writeBinaryGrid)
441 .def("getX",&speckley::SpeckleyDomain::getX, ":return: locations in the FEM nodes\n\n"
442 ":rtype: `Data`")
443 .def("getNormal",&speckley::SpeckleyDomain::getNormal,
444 ":return: boundary normals at the quadrature point on the face elements\n"
445 ":rtype: `Data`")
446 .def("getSize",&speckley::SpeckleyDomain::getSize,":return: the element size\n"
447 ":rtype: `Data`")
448 .def("setTagMap",&speckley::SpeckleyDomain::setTagMap,args("name","tag"),
449 "Give a tag number a name.\n\n:param name: Name for the tag\n:type name: ``string``\n"
450 ":param tag: numeric id\n:type tag: ``int``\n:note: Tag names must be unique within a domain")
451 .def("getTag",&speckley::SpeckleyDomain::getTag,args("name"),":return: tag id for "
452 "``name``\n:rtype: ``string``")
453 .def("isValidTagName",&speckley::SpeckleyDomain::isValidTagName,args("name"),
454 ":return: True if ``name`` corresponds to a tag, otherwise False\n:rtype: ``bool``")
455 .def("showTagNames",&speckley::SpeckleyDomain::showTagNames,":return: A space separated list of tag names\n:rtype: ``string``")
456 .def("getMPISize",&speckley::SpeckleyDomain::getMPISize,":return: the number of processes used for this `Domain`\n:rtype: ``int``")
457 .def("getMPIRank",&speckley::SpeckleyDomain::getMPIRank,":return: the rank of this process\n:rtype: ``int``")
458 .def("MPIBarrier",&speckley::SpeckleyDomain::MPIBarrier,"Wait until all processes have reached this point")
459 .def("onMasterProcessor",&speckley::SpeckleyDomain::onMasterProcessor,":return: True if this code is executing on the master process\n:rtype: `bool`");
460 /* These two class exports are necessary to ensure that the extra methods added by speckley make it to python.
461 * This change became necessary when the Brick and Rectangle constructors turned into factories instead of classes */
462 class_<speckley::Brick, bases<speckley::SpeckleyDomain> >("SpeckleyBrick", "", no_init);
463 class_<speckley::Rectangle, bases<speckley::SpeckleyDomain> >("SpeckleyRectangle", "", no_init);
464 class_<speckley::AbstractAssembler, speckley::Assembler_ptr, boost::noncopyable >
465 ("AbstractAssembler", "", no_init);
466 }
467

  ViewVC Help
Powered by ViewVC 1.1.26