~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to site-packages/dolfin/compile_functions.py

  • Committer: Kent-Andre Mardal
  • Date: 2008-06-07 09:10:29 UTC
  • mfrom: (2668.14.7 trunk)
  • mto: (2668.1.50 trunk)
  • mto: This revision was merged to the branch mainline in revision 2670.
  • Revision ID: kent-and@simula.no-20080607091029-ulyrrbk8kqd9xjw0
merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
"This module provides functionality for compilation of strings as dolfin Functions."
 
2
 
 
3
__author__ = "Martin Sandve Alnes (martinal@simula.no)"
 
4
__date__ = "2008-06-04 -- 2008-06-04"
 
5
__copyright__ = "Copyright (C) 2008-2008 Martin Sandve Alnes"
 
6
__license__  = "GNU LGPL Version 2.1"
 
7
 
 
8
import re
 
9
import numpy
 
10
import instant
 
11
 
 
12
__all__ = ["compile_functions",]
 
13
 
 
14
# FIXME: Extend this list, needed to autodetect variable names that are not builtins
 
15
_builtins = [
 
16
             # local symbols:
 
17
             "v", "x", "_dim", "pi",
 
18
             # cmath funcions:
 
19
             "cos", "sin", "tan", "acos", "asin", "atan", "atan2",
 
20
             "cosh", "sinh", "tanh",
 
21
             "exp", "frexp", "ldexp", "log", "log10", "modf",
 
22
             "pow", "sqrt", "ceil", "fabs", "floor", "fmod",
 
23
            ]
 
24
 
 
25
# Add utility code here
 
26
_header = """
 
27
// cmath functions
 
28
using std::cos;
 
29
using std::sin;
 
30
using std::tan;
 
31
using std::acos;
 
32
using std::asin;
 
33
using std::atan;
 
34
using std::atan2;
 
35
using std::cosh;
 
36
using std::sinh;
 
37
using std::tanh;
 
38
using std::exp;
 
39
using std::frexp;
 
40
using std::ldexp;
 
41
using std::log;
 
42
using std::log10;
 
43
using std::modf;
 
44
using std::pow;
 
45
using std::sqrt;
 
46
using std::ceil;
 
47
using std::fabs;
 
48
using std::floor;
 
49
using std::fmod;
 
50
 
 
51
const double pi = acos(-1.0);
 
52
"""
 
53
 
 
54
 
 
55
_function_template = """
 
56
class %(classname)s: public dolfin::Function
 
57
{
 
58
public:
 
59
  unsigned _dim;
 
60
%(members)s
 
61
 
 
62
  %(classname)s(dolfin::Mesh & mesh):
 
63
    dolfin::Function(mesh)
 
64
  {
 
65
    _dim = mesh.geometry().dim();
 
66
%(constructor)s
 
67
  }
 
68
 
 
69
  void eval(real* v, const real* x) const
 
70
  {
 
71
%(eval)s
 
72
  }
 
73
%(rank)s%(dim)s
 
74
};
 
75
"""
 
76
 
 
77
_ranktemplate = """\
 
78
  dolfin::uint rank() const
 
79
  {
 
80
    return %d;
 
81
  }
 
82
"""
 
83
 
 
84
_dimtemplate = """\
 
85
  dolfin::uint dim(dolfin::uint i) const
 
86
  {
 
87
    switch(i) {
 
88
%s
 
89
    }
 
90
    throw std::runtime_error("Invalid dimension i in dim(i).");
 
91
  }
 
92
"""
 
93
 
 
94
def expression_to_function(e, classname):
 
95
    "Generates code for a dolfin::Function subclass for a single expression."
 
96
    # Get shape from e and make e a flat tuple of strings
 
97
    if isinstance(e, str):
 
98
        e = (e,)
 
99
        shape = ()
 
100
    elif isinstance(e, tuple):
 
101
        if isinstance(e[0], str):
 
102
            shape = (len(e),)
 
103
        elif isinstance(e[0], tuple):
 
104
            shape = (len(e),len(e[0]))
 
105
            assert isinstance(e[0][0], str)
 
106
            e = sum(e, ())
 
107
    else:
 
108
        raise RuntimeError("Invalid expression %s" % e)
 
109
 
 
110
    # Autodetect variables from function strings
 
111
    variables = set()
 
112
    for c in e:
 
113
        # Find groups of connected alphanumeric letters
 
114
        variables.update(re.findall(r"([\w]+)", c))
 
115
    variables.difference_update(_builtins)
 
116
    numerals = [v for v in variables if v[0] in "0123456789"]
 
117
    variables.difference_update(numerals)
 
118
    
 
119
    # Generate code for member variables
 
120
    memberscode = "\n".join("  double %s;" % name for name in variables)
 
121
    
 
122
    # Generate constructor code for initialization of member variables
 
123
    constructorcode = "\n".join("    %s = 0.0;" % name for name in variables)
 
124
    
 
125
    # Generate code for rank and dim
 
126
    if len(shape) > 0:
 
127
        rankcode = _ranktemplate % len(shape)
 
128
        cases = "\n".join("      case %d: return %d;" % (d,n) for (d,n) in enumerate(shape))
 
129
        dimcode = _dimtemplate % cases
 
130
    else:
 
131
        dimcode  = ""
 
132
        rankcode = ""
 
133
    
 
134
    # Generate code for the actual function evaluation
 
135
    evalcode = "\n".join("    v[%d] = %s;" % (i, c) for (i,c) in enumerate(e))
 
136
    
 
137
    # Connect the code fragments using the function template code
 
138
    fragments = {}
 
139
    fragments["classname"] = classname
 
140
    fragments["members"]   = memberscode
 
141
    fragments["rank"]      = rankcode
 
142
    fragments["dim"]       = dimcode
 
143
    fragments["eval"]      = evalcode
 
144
    fragments["constructor"] = constructorcode
 
145
    code = _function_template % fragments
 
146
    return code
 
147
 
 
148
 
 
149
_function_count = 0
 
150
def generate_functions(expressions):
 
151
    "Generates code for dolfin::Function subclasses for a list of expressions."
 
152
    global _function_count
 
153
    assert isinstance(expressions, list)
 
154
    code = ""
 
155
    classnames = []
 
156
    for e in expressions:
 
157
        classname = "function_%d" % _function_count
 
158
        code += expression_to_function(e, classname)
 
159
        classnames.append(classname)
 
160
        _function_count += 1
 
161
    return code, classnames
 
162
 
 
163
# NB! This code is highly dependent on the dolfin swig setup!
 
164
_additional_declarations = r"""
 
165
%rename(cpp_Function) Function;
 
166
%include exception.i
 
167
%include cpointer.i
 
168
%pointer_class(int, intp);
 
169
%pointer_class(double, doublep);
 
170
%include std_vector.i
 
171
%template(STLVectorFunctionPtr) std::vector<dolfin::Function *>;
 
172
%import dolfin.i
 
173
"""
 
174
    
 
175
_additional_definitions  = """
 
176
#include <dolfin.h>
 
177
using namespace dolfin;
 
178
#include <numpy/arrayobject.h>
 
179
"""
 
180
 
 
181
_module_count = 0
 
182
def compile_function_code(code, mesh, classnames=None, module_name=None):
 
183
    # Create unique module name for this application run
 
184
    global _module_count, _header, _additional_definitions, _additional_declarations
 
185
    if module_name is None:
 
186
        module_name = "functions_%d" % _module_count
 
187
        _module_count += 1
 
188
 
 
189
    # Autodetect classnames:
 
190
    _classnames = re.findall(r"class[ ]+([\w]+).*", code)
 
191
    # Just a little assertion for safety:
 
192
    if classnames is None:
 
193
        classnames = _classnames
 
194
    else:
 
195
        assert all(a == b for (a,b) in zip(classnames, _classnames))
 
196
    
 
197
    # Get system configuration   
 
198
    (includes, flags, libraries, libdirs) = instant.header_and_libs_from_pkgconfig("dolfin")
 
199
    dolfin_include_dir = includes[0] # FIXME: is this safe?
 
200
    numpy_dir = numpy.get_include()
 
201
    includes.append(numpy_dir)
 
202
    
 
203
    sysheaders = ["cmath", "iostream", "stdexcept",
 
204
                  "dolfin.h", "dolfin/mesh/Mesh.h", "dolfin/function/Function.h"]
 
205
    
 
206
    # FIXME: use dolfin flags?
 
207
    cppargs = flags
 
208
    cppargs = " ".join(flags)
 
209
    cppargs = ""
 
210
    
 
211
    # Let swig see the installed dolfin swig files
 
212
    swigopts = "-c++ -I%s -I%s/dolfin/swig" % (dolfin_include_dir, dolfin_include_dir)
 
213
    
 
214
    # Compile extension module with instant
 
215
    compiled_module = instant.create_extension(\
 
216
             code           = _header + code,
 
217
             module         = module_name,
 
218
             swigopts       = swigopts,
 
219
             system_headers = sysheaders,
 
220
             include_dirs   = includes,
 
221
             cppargs        = cppargs,
 
222
             libraries      = libraries,
 
223
             library_dirs   = libdirs,
 
224
             additional_definitions  = _additional_definitions,
 
225
             additional_declarations = _additional_declarations
 
226
             )
 
227
    compiled_module = __import__(module_name)
 
228
    
 
229
    # Construct instances of the compiled functor classes
 
230
    functions = [eval("compiled_module.%s(mesh)" % name) for name in classnames]
 
231
    return functions
 
232
 
 
233
 
 
234
def compile_functions(expressions, mesh):
 
235
    """Compiles C++ string expressions into dolfin::Function instances.
 
236
    
 
237
    The variable 'expressions' can either be a str or a list.
 
238
    
 
239
    If 'expressions' is a str, it is interpreted as a C++ string
 
240
    with complete implementations of subclasses of dolfin::Function.
 
241
    The compiled functions returned will be in the same order 
 
242
    as they are defined in this code.
 
243
 
 
244
    If it is a list, each item of the list is interpreted as
 
245
    a function expression, and the compiled functions returned
 
246
    will be in the same order as they occur in this list.
 
247
    
 
248
    Each expression item in the list can either be a str,
 
249
    in which case it is interpreted as a scalar expression
 
250
    and a scalar Function is generated, or it can be a tuple.
 
251
    
 
252
    A tuple of str objects is interpreted as a vector expression,
 
253
    and a rank 1 Function is generated.
 
254
    
 
255
    A tuple of tuples of str objects is interpreted as a matrix
 
256
    expression, and a rank 2 Function is generated.
 
257
 
 
258
    If an expression string contains a name, it is assumed to
 
259
    be a scalar variable name, and is added as a public member
 
260
    of the generated function. The exceptions are set in the
 
261
    variable dolfin.compile_functions._builtins."""
 
262
    #, which contains:
 
263
    #    %s
 
264
    #""" % "\n".join("        " + b for b in _builtins)
 
265
    if isinstance(expressions, list):
 
266
        code, classnames = generate_functions(expressions)
 
267
        functions = compile_function_code(code, mesh, classnames)
 
268
    else:
 
269
        functions = compile_function_code(expressions, mesh)
 
270
    return functions
 
271
 
 
272
 
 
273
if __name__ == "__main__":
 
274
    code, cn = generate_functions([
 
275
                              "exp(alpha)",
 
276
                              ("sin(x[0])", "cos(x[1])", "0.0"),
 
277
                              (("sin(x[0])", "cos(x[1])"), ("0.0", "1.0"))
 
278
                             ])
 
279
    print code
 
280
    print cn
 
281
 
 
282