Package sfc :: Package codegeneration :: Module finiteelementcg
[hide private]
[frames] | no frames]

Source Code for Module sfc.codegeneration.finiteelementcg

  1  #!/usr/bin/env python 
  2  # -*- coding: utf-8 -*- 
  3  """ 
  4  This module contains code generation tools for the ufc::finite_element class. 
  5  """ 
  6   
  7  # Copyright (C) 2008-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  8  # 
  9  # This file is part of SyFi. 
 10  # 
 11  # SyFi is free software: you can redistribute it and/or modify 
 12  # it under the terms of the GNU General Public License as published by 
 13  # the Free Software Foundation, either version 2 of the License, or 
 14  # (at your option) any later version. 
 15  # 
 16  # SyFi is distributed in the hope that it will be useful, 
 17  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 18  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 19  # GNU General Public License for more details. 
 20  # 
 21  # You should have received a copy of the GNU General Public License 
 22  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 23  # 
 24  # First added:  2008-08-13 
 25  # Last changed: 2009-03-19 
 26   
 27  from itertools import izip 
 28  import ufl 
 29  import swiginac 
 30   
 31  from sfc.codegeneration.codeformatting import indent, CodeFormatter, gen_switch, gen_token_assignments, gen_const_token_definitions 
 32  from sfc.symbolic_utils import symbol, symbols, TempSymbolContext 
 33  from sfc.geometry import gen_geometry_code 
 34  from sfc.common import sfc_warning, sfc_assert, sfc_error, sfc_info 
 35  from sfc.common.utilities import unique 
 36  import ufc_utils 
 37   
38 -def cse(*args):
39 sfc_error("FIXME")
40
41 -class FiniteElementCG:
42 """Code generator for ufc::finite_element implementations."""
43 - def __init__(self, elementrep):
44 self.rep = elementrep 45 self.classname = elementrep.finite_element_classname 46 self.signature = repr(self.rep.ufl_element) 47 self.options = self.rep.options.code.finite_element 48 self.DN_order = self.options.evaluate_basis_derivatives_order
49
50 - def hincludes(self):
51 l = [] 52 return l
53
54 - def cincludes(self):
55 l = [] 56 return l
57
58 - def local_includes(self):
59 return unique("%s.h" % e.finite_element_classname 60 for e in self.rep.sub_elements)
61
62 - def header_template(self):
63 return ufc_utils.finite_element_header
64
65 - def impl_template(self):
66 return ufc_utils.finite_element_implementation
67
68 - def generate_code_dict(self):
69 # generate the code components: 70 vars = { 71 "classname" : self.classname, 72 "constructor" : indent(self.gen_constructor()), 73 "constructor_arguments" : indent(self.gen_constructor_arguments()), 74 "initializer_list" : indent(self.gen_initializer_list()), 75 "destructor" : indent(self.gen_destructor()), 76 "create" : indent(self.gen_create()), 77 "signature" : indent(self.gen_signature()), 78 "cell_shape" : indent(self.gen_cell_shape()), 79 "geometric_dimension" : indent(self.gen_geometric_dimension()), 80 "topological_dimension" : indent(self.gen_topological_dimension()), 81 "space_dimension" : indent(self.gen_space_dimension()), 82 "value_rank" : indent(self.gen_value_rank()), 83 "value_dimension" : indent(self.gen_value_dimension()), 84 "map_from_reference_cell" : indent(self.gen_map_from_reference_cell()), 85 "map_to_reference_cell" : indent(self.gen_map_to_reference_cell()), 86 "evaluate_basis" : indent(self.gen_evaluate_basis()), 87 "evaluate_basis_all" : indent(self.gen_evaluate_basis_all()), 88 "evaluate_basis_derivatives" : indent(self.gen_evaluate_basis_derivatives()), 89 "evaluate_basis_derivatives_all" : indent(self.gen_evaluate_basis_derivatives_all()), 90 "evaluate_dof" : indent(self.gen_evaluate_dof()), 91 "evaluate_dofs" : indent(self.gen_evaluate_dofs()), 92 "interpolate_vertex_values" : indent(self.gen_interpolate_vertex_values()), 93 "num_sub_elements" : indent(self.gen_num_sub_elements()), 94 "create_sub_element" : indent(self.gen_create_sub_element()), 95 "members" : indent(self.gen_members()), 96 } 97 return vars
98
99 - def generate_support_code(self):
100 return ""
101
102 - def gen_constructor(self):
103 return ""
104
106 return ""
107
108 - def gen_initializer_list(self):
109 return ""
110
111 - def gen_destructor(self):
112 return ""
113
114 - def gen_create(self):
115 code = "return new %s();" % self.classname 116 return code
117
118 - def gen_signature(self):
119 """const char* signature() const""" 120 return 'return "%s";' % self.signature
121
122 - def gen_cell_shape(self):
123 """shape cell_shape() const""" 124 return "return ufc::%s;" % self.rep.cell.shape
125
126 - def gen_geometric_dimension(self):
127 return "return %d;" % self.rep.cell.nsd
128
130 return "return %d;" % self.rep.cell.nsd
131
132 - def gen_space_dimension(self):
133 """unsigned int space_dimension() const""" 134 return "return %d;" % self.rep.local_dimension
135
136 - def gen_value_rank(self):
137 """unsigned int value_rank() const""" 138 return "return %d;" % self.rep.value_rank
139
140 - def gen_value_dimension(self):
141 """unsigned int value_dimension(unsigned int i) const""" 142 if self.rep.value_rank == 0: 143 code = 'throw std::runtime_error("Rank 0 value has no dimension.");\n' 144 #elif self.rep.value_rank == 1: 145 # code = "return %d;" % self.rep.value_shape[0] 146 else: 147 dims = self.rep.value_shape 148 cases = [(i, "return %d;" % d) for (i,d) in enumerate(dims)] 149 code = gen_switch("i", cases) 150 code += 'throw std::runtime_error("Invalid dimension for rank %d value.");\n' % self.rep.value_rank 151 152 return code
153
155 return 'throw std::runtime_error("Not implemented.");' # FIXME
156
158 return 'throw std::runtime_error("Not implemented.");' # FIXME
159
160 - def gen_evaluate_basis(self):
161 """void evaluate_basis(unsigned int i, 162 double* values, 163 const double* coordinates, 164 const cell& c) const 165 """ 166 if not self.options.enable_evaluate_basis: 167 return 'throw std::runtime_error("evaluate_basis not implemented.");' 168 169 nsd = self.rep.cell.nsd 170 nbf = self.rep.local_dimension 171 value_shape = self.rep.value_shape 172 value_size = self.rep.value_size 173 174 if self.rep.ufl_element.family() == "Real": 175 code = [] 176 if value_size > 1: 177 code += ["memset(values, 0, sizeof(double)*%d);" % value_size] 178 code += ['values[i] = 1.0;'] 179 return '\n'.join(code) 180 181 # symbols for output values 182 val_sym = symbols(["values[%d]" % d for d in range(value_size)]) 183 184 # generate code body 185 code = CodeFormatter() 186 187 code += gen_geometry_code(nsd, detG=False, GinvT=True) 188 189 coordinates = swiginac.matrix(nsd, 1, symbols(["coordinates[%d]" % i for i in range(nsd)])) 190 xi = self.rep.GinvT.transpose().mul(coordinates - self.rep.p0).evalm() 191 192 for i in range(nsd): 193 code += "const double %s = %s;" % (("x","y","z")[i], xi[i].printc()) 194 195 if value_size > 1: 196 code += "memset(values, 0, sizeof(double)*%d);" % value_size 197 code += "" 198 199 # begin switch 200 if nbf>1: 201 code.begin_switch("i") 202 203 # generate one case for each basis function 204 for i in range(nbf): 205 206 # make token list for basis function i 207 values = [] 208 for c in self.rep.value_components: 209 values.append( self.rep.basis_function(i, c) ) 210 values_tokens = [] 211 for d in range(value_size): 212 if not values[d].expand().is_zero(): 213 values_tokens.append( (val_sym[d], values[d]) ) 214 # now values_tokens is a token list with "values[i] = expression" for all output values 215 216 if self.options.optimize_basis: 217 # split tokens into temporary variables (1) and output variables (2) 218 temp_symbol = TempSymbolContext() 219 values_tokens1 = [] 220 values_tokens2 = [] 221 for s, e in values_tokens: # handle s = e: 222 ts = temp_symbol() # construct ts 223 values_tokens1.append( (ts, e) ) # set ts = e 224 values_tokens2.append( (s, ts) ) # set s = ts 225 values_tokens = None 226 227 # optimize temporary variable list 228 values_tokens1, repmap = cse(values_tokens1, temp_symbol) 229 230 # generate code 231 values_code = gen_const_token_definitions(values_tokens1) + "\n" 232 values_code += gen_token_assignments(values_tokens2) 233 #values_codelines= chain(const_token_definitions(values_tokens1), token_assignments(values_tokens2)) 234 else: 235 values_code = gen_token_assignments(values_tokens) 236 #values_codelines = token_assignments(values_tokens) 237 238 # generate case code 239 if nbf>1: 240 code.begin_case(i, braces=True) 241 code.new_line( indent(values_code) ) 242 code.end_case() 243 #with code.case(i, braces=True): 244 # code.add_lines(values_codelines) 245 else: 246 code += values_code 247 248 if nbf>1: 249 code.end_switch() 250 251 return str(code)
252
253 - def gen_evaluate_basis_all(self): # TODO: implement optimized version of this
254 """/// Evaluate all basis functions at given point in cell 255 virtual void evaluate_basis_all(double* values, 256 const double* coordinates, 257 const cell& c) const 258 """ 259 260 code = CodeFormatter() 261 code += "for(unsigned i = 0; i < %d; i++)" % self.rep.local_dimension 262 code.begin_block() 263 code += "evaluate_basis(i, values+i*%d, coordinates, c);" % self.rep.value_size 264 code.end_block() 265 return str(code)
266
267 - def gen_evaluate_basis_derivatives(self):
268 """/// Evaluate order n derivatives of basis function i at given point in cell 269 void evaluate_basis_derivatives(unsigned int i, 270 unsigned int n, 271 double* values, 272 const double* coordinates, 273 const ufc::cell& c) const 274 """ 275 if not self.options.enable_evaluate_basis_derivatives: 276 return 'throw std::runtime_error("evaluate_basis_derivatives not implemented.");' 277 278 nsd = self.rep.cell.nsd 279 nbf = self.rep.local_dimension 280 value_shape = self.rep.value_shape 281 value_size = self.rep.value_size 282 283 if self.rep.ufl_element.family() == "Real": 284 return '\n'.join('values[%d] = 0.0;' % d for d in range(value_size)) 285 286 sfc_assert(self.DN_order <= 2, "Don't support computing higher order derivatives (yet).") 287 288 code = CodeFormatter() 289 code.begin_if("n > 2") 290 code += 'throw std::runtime_error("evaluate_basis_derivatives not implemented for the wanted derivative order.");' 291 code.end_if() 292 293 # define GinvT 294 code += gen_geometry_code(nsd, detG=False, GinvT=True) 295 296 # define x,y,z (spatial symbols) 297 p = symbols(["x", "y", "z"][:nsd]) 298 coords = symbols(["coordinates[%d]" % i for i in range(nsd)]) 299 code += gen_const_token_definitions( izip(p, coords) ) 300 301 # switch on derivative order 302 code.begin_switch("n") 303 304 for order in range(1, self.DN_order+1): 305 code.begin_case(order, braces=True) 306 307 # zero output array before filling nonzeros 308 num_derivatives = nsd ** order 309 do_memset = (value_size*num_derivatives) > 6 310 if do_memset: 311 code += "memset(values, 0, sizeof(double) * %d * %d);" % (value_size, num_derivatives) 312 code += "" 313 314 # switch on basis function number 315 code.begin_switch("i") 316 for ibf in range(nbf): 317 code.begin_case(ibf, braces=True) 318 temp_symbol = TempSymbolContext() 319 320 DN_tokens = [] 321 all_directions = ufl.permutation.compute_permutations(self.DN_order, nsd) 322 for (j,directions) in enumerate(all_directions): 323 for (k,c) in enumerate(self.rep.value_components): 324 DN = self.rep.basis_function_derivative(ibf, c, directions) 325 # DN is now the expression for the derivative wrt directions of component c of basis function ibf 326 if not (do_memset and DN.expand().is_zero()): 327 values_sym = symbol("values[%d * %d + %d]" % (j, value_size, k)) 328 DN_tokens.append( (values_sym, DN) ) 329 330 if self.options.optimize_basis: 331 # split tokens into temporary variables and output variables 332 DN_tokens1 = [] 333 DN_tokens2 = [] 334 for s, e in DN_tokens: # s = e 335 ts = temp_symbol() # -> 336 DN_tokens1.append( (ts, e) ) # ts = e 337 DN_tokens2.append( (s, ts) ) # s = ts 338 # optimize temporary variable list 339 DN_tokens1, repmap = cse(DN_tokens1, temp_symbol) 340 # generate code 341 code += gen_const_token_definitions(DN_tokens1) 342 code += gen_token_assignments(DN_tokens2) 343 else: 344 code += gen_token_assignments(DN_tokens) 345 code.end_case() 346 code.end_switch() 347 code.end_case() 348 code += "default:" 349 code.indent() 350 code += 'throw std::runtime_error("Derivatives of this order are not supported in evaluate_basis_derivatives.");' 351 code.dedent() 352 code.end_switch() 353 return str(code)
354
355 - def gen_evaluate_basis_derivatives_all(self): # TODO: implement optimized version of this
356 """/// Evaluate order n derivatives of all basis functions at given point in cell 357 virtual void evaluate_basis_derivatives_all(unsigned int n, 358 double* values, 359 const double* coordinates, 360 const cell& c) const 361 """ 362 code = CodeFormatter() 363 code.begin_switch("n") 364 for order in range(1, self.DN_order+1): 365 code.begin_case(order, braces=True) 366 offset = self.rep.value_size * (self.rep.cell.nsd ** order) 367 code += "for(unsigned i = 0; i < %d; i++)" % self.rep.local_dimension 368 code.begin_block() 369 code += "evaluate_basis_derivatives(n, i, values+i*%d, coordinates, c);" % offset 370 code.end_block() 371 code.end_case() 372 code += "default:" 373 code.indent() 374 code += 'throw std::runtime_error("Derivatives of this order are not implemented in evaluate_basis_derivatives_all.");' 375 code.dedent() 376 code.end_switch() 377 return str(code) 378
379 - def gen_evaluate_dof(self):
380 """double evaluate_dof(unsigned int i, 381 const function& f, 382 const cell& c) const 383 This implementation is general for all elements with point evaluation dofs. 384 385 TODO: Implement support for normal and tangential component dofs. 386 TODO: Implement for elements without point evaluation dofs. 387 """ 388 # some useful variables 389 nsd = self.rep.cell.nsd 390 nbf = self.rep.local_dimension 391 392 # initial code 393 code = CodeFormatter() 394 code.new_text( gen_geometry_code(nsd, detG=False) ) 395 code += "double v[%d];" % self.rep.value_size 396 code += "double x[%d];" % nsd 397 398 # fill global coordinates of dof in x[] 399 if nbf == 1: 400 # skip the switch if only one basis function 401 for k in range(nsd): 402 code += "x[%d] = %s;" % (k, self.rep.dof_x[0][k].printc()) 403 else: 404 code.begin_switch("i") 405 for i in range(nbf): 406 code.begin_case(i) 407 # compute dof coordinate i 408 for k in range(nsd): 409 code += "x[%d] = %s;" % (k, self.rep.dof_x[i][k].printc()) 410 code.end_case() 411 code.end_switch() 412 413 # Evaluate the function (this evaluates all value components!) 414 code += "f.evaluate(v, x, c);" 415 416 # dofs for a single sub element are numbered contiguously: 417 # i = (nbf / valsize) * value_component 418 if nbf == 1: 419 code += "return v[0];" 420 else: 421 code += "return v[i / %d];" % (nbf // self.rep.value_size) 422 423 return str(code)
424
425 - def gen_evaluate_dofs(self): # TODO: implement optimized version of this
426 """/// Evaluate linear functionals for all dofs on the function f 427 virtual void evaluate_dofs(double* values, 428 const function& f, 429 const cell& c) const 430 """ 431 code = CodeFormatter() 432 code += "for(unsigned i=0; i<%d; i++)" % self.rep.local_dimension 433 code.begin_block() 434 code += "values[i] = evaluate_dof(i, f, c);" 435 code.end_block() 436 return str(code) 437
438 - def gen_interpolate_vertex_values(self):
439 """void interpolate_vertex_values(double* vertex_values, 440 const double* dof_values, 441 const cell& c) const 442 """ 443 # some helper variables 444 nbf = self.rep.local_dimension 445 nsd = self.rep.cell.nsd 446 nv = self.rep.cell.num_entities[0] 447 value_size = self.rep.value_size 448 value_shape = self.rep.value_shape 449 450 # symbols for input array entries 451 dof_values_sym = symbols("dof_values[%d]" % i for i in xrange(nbf)) 452 vertex_values_sym = symbols("vertex_values[%d]" % i for i in xrange(nv*value_size)) 453 454 # the spatial symbols u is expressed in 455 p = self.rep.p 456 457 # construct expressions for the linear combinations of basis functions 458 u = [] 459 for component in self.rep.value_components: 460 u.append( sum(self.rep.basis_function(j, component)*dof_values_sym[j] for j in range(nbf)) ) 461 462 # for each vertex 463 vertex_values = [] 464 repmap = swiginac.exmap() 465 for i in range(nv): 466 # replacement map for coordinates 467 vx = self.rep.polygon.vertex(i) 468 for k in range(nsd): 469 repmap[p[k]] = vx[k] 470 # evaluate functions for each component in coordinate 471 for uval in u: 472 vertex_values.append(uval.subs(repmap)) 473 474 code = gen_token_assignments( izip(vertex_values_sym, vertex_values) ) 475 return code
476
477 - def gen_num_sub_elements(self):
478 """unsigned int num_sub_elements() const""" 479 return "return %d;" % len(self.rep.sub_elements)
480
481 - def gen_create_sub_element(self):
482 """finite_element* create_sub_element(unsigned int i) const""" 483 if len(self.rep.sub_elements) > 1: 484 code = CodeFormatter() 485 code.begin_switch("i") 486 for i, fe in enumerate(self.rep.sub_elements): 487 code += "case %d: return new %s();" % (i, fe.finite_element_classname) 488 code.end_switch() 489 code += 'throw std::runtime_error("Invalid index in create_sub_element.");' 490 else: 491 code = "return new %s();" % self.classname # FIXME: Should we throw error here instead now? 492 return str(code)
493
494 - def gen_members(self):
495 return ""
496 #code = CodeFormatter() 497 #code += "public:"; 498 #code += "protected:"; 499 #code.indent() 500 #code += "unsigned int foo = 0;" 501 #code.dedent() 502 #return str(code) 503