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