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

Source Code for Module sfc.codegeneration.integralcg_regular

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains code generation functionality for the ufc::*_integral classes. 
  4  """ 
  5   
  6  # Copyright (C) 2008-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  7  # 
  8  # This file is part of SyFi. 
  9  # 
 10  # SyFi is free software: you can redistribute it and/or modify 
 11  # it under the terms of the GNU General Public License as published by 
 12  # the Free Software Foundation, either version 2 of the License, or 
 13  # (at your option) any later version. 
 14  # 
 15  # SyFi is distributed in the hope that it will be useful, 
 16  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 17  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 18  # GNU General Public License for more details. 
 19  # 
 20  # You should have received a copy of the GNU General Public License 
 21  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 22  # 
 23  # First added:  2008-08-13 
 24  # Last changed: 2009-03-10 
 25   
 26  from itertools import tee 
 27  import swiginac 
 28  import ufc_utils 
 29   
 30  from sfc.common.output import sfc_debug, sfc_error, sfc_assert, sfc_warning 
 31  from sfc.common.utilities import indices_subset, unique 
 32   
 33  from sfc.codegeneration.codeformatting import (indent, CodeFormatter, 
 34                                                 gen_symbol_declarations, 
 35                                                 gen_const_token_definitions, 
 36                                                 gen_token_declarations, 
 37                                                 gen_token_assignments, 
 38                                                 gen_token_additions, 
 39                                                 gen_token_prints) 
 40   
 41  from sfc.codegeneration.support_code import (generate_sign_correction_support_code, 
 42                                               gen_A_reset_block) 
 43   
 44  from sfc.codegeneration.integralcgbase import IntegralCGBase 
 45   
 46   
 47  from sfc.quadrature import (gen_quadrature_rule_definition, 
 48                              gen_quadrature_rule_definitions) 
49 -class RegularIntegralCG(IntegralCGBase):
50 "Shared RegularIntegralCG code for several base classes."
51 - def __init__(self, itgrep):
52 IntegralCGBase.__init__(self, itgrep)
53
54 - def generate_support_code(self):
55 code = [] 56 57 # Write quadrature rule if we have one: 58 fr = self.itgrep.formrep 59 # TODO: Place quadrature rules in a separate shared file? 60 if self.itgrep._on_facet: 61 if fr.facet_quad_rules is not None: 62 code += [gen_quadrature_rule_definitions(fr.facet_quad_rules)] 63 else: 64 if fr.quad_rule is not None: 65 code += [gen_quadrature_rule_definition(fr.quad_rule)] 66 67 # Write sign correction code if needed 68 code += [generate_sign_correction_support_code(fr)] 69 70 return '\n'.join(code)
71
72 - def gen_pre_debug_code(self):
73 if not self.options.enable_debug_code: 74 return "" 75 code = CodeFormatter() 76 code.begin_debug() 77 code.begin_block() 78 code += 'std::cout << std::endl;' 79 code += 'std::cout << "SFC DEBUGGING OUTPUT: " << std::endl;' 80 code += 'std::cout << "void %s::tabulate_tensor(...)" << std::endl;' % self.itgrep.classname 81 code += 'std::cout << "{" << std::endl;' 82 fr = self.itgrep.formrep 83 fd = fr.formdata 84 for k in range(fr.num_coefficients): 85 element = fd.elements[k] 86 rep = fr.element_reps[element] 87 code += 'std::cout << " w[%d][:] = [ ";' % k 88 code += "for(int j=0; j<%d; j++)" % rep.local_dimension 89 code.begin_block() 90 code += 'std::cout << w[%d][j] << ", ";' % k 91 code.end_block() 92 code += 'std::cout << " ]" << std::endl;' 93 code += 'std::cout << " now computing element tensor..." << std::endl;' 94 code.end_block() 95 code.end_debug() 96 return str(code)
97
98 - def gen_post_debug_code(self):
99 if not self.options.enable_debug_code: 100 return "" 101 code = CodeFormatter() 102 code.begin_debug() 103 code.begin_block() 104 code += 'std::cout << " ... done computing element tensor." << std::endl;' 105 code += "for(int k=0; k<%d; k++)" % len(self.itgrep.indices) 106 code.begin_block() 107 # TODO: enhance this output by printing indices 108 code += 'std::cout << " A[" << k << "]" << " = " << A[k] << std::endl;' 109 code.end_block() 110 code += 'std::cout << "}" << std::endl;' 111 code.end_block() 112 code.end_debug() 113 return str(code)
114
115 - def gen_partition_block(self, data, deps, basis_functions):
116 "Generate code for a partition of the linearized computational graph." 117 code = gen_token_assignments( self.itgrep.iter_partition(data, deps, 118 basis_functions) ) 119 if code: 120 deps = ", ".join(deps) 121 if basis_functions: 122 deps += "; using basis functions %s" % str(basis_functions) 123 code = "// Partition depending on %s\n{\n%s\n}" % (deps, indent(code)) 124 return code
125
126 - def gen_geometry_block(self):
127 sfc_debug("Entering RegularIntegralCG.gen_geometry_block") 128 129 code = CodeFormatter() 130 131 # --- Sign correction code: # TODO 132 #sign_sym = itgrep.sign_sym 133 #for fe_name in sign_sym.keys(): 134 # nbf = len(sign_sym[fe_name]) 135 # code += use_sign_code2(fe_name, nbf, nsd) 136 137 # --- Generate code for for initial geometry tokens 138 # (vertex coordinates, affine map, normal vector, ...) 139 geometry_tokens = self.itgrep.iter_geometry_tokens() 140 141 # if we want to print debug tokens, duplicate the iterator 142 if self.options.enable_debug_code: 143 geometry_tokens, debug_geometry_tokens = tee(geometry_tokens) 144 145 code += "// Geometric quantities" 146 code += gen_const_token_definitions(geometry_tokens) 147 148 if self.options.enable_debug_code: 149 code += "" 150 code.begin_debug() 151 code += gen_token_prints(debug_geometry_tokens) 152 code.end_debug() 153 code += "" 154 155 # --- Facet tokens code: 156 if self.itgrep._on_facet: 157 fr = self.itgrep.formrep 158 code += "// Geometric quantities on each facet" 159 for facet in range(fr.cell.num_facets): 160 facet_tokens = self.itgrep.iter_facet_tokens(facet) 161 162 if facet == 0: 163 # if we want to print debug tokens, duplicate the iterator 164 if self.options.enable_debug_code: 165 facet_tokens, debug_facet_tokens = tee(facet_tokens) 166 167 # duplicate the iterator to generate declarations outside switch 168 facet_tokens, facet_tokens2 = tee(facet_tokens) 169 code += gen_token_declarations(facet_tokens2) 170 code.begin_switch("facet") 171 172 code.begin_case(facet, braces=True) 173 code += "// Geometric quantities on facet %d" % facet 174 code += gen_token_assignments(facet_tokens) 175 code.end_case() 176 code += "default:" 177 code.indent() 178 code += 'throw std::runtime_error("Invalid facet number.");' 179 code.dedent() 180 code.end_switch() 181 182 if self.options.enable_debug_code: 183 code += "" 184 code.begin_debug() 185 code += 'std::cout << "facet = " << facet << std::endl;' 186 code += gen_token_prints(debug_facet_tokens) 187 code.end_debug() 188 code += "" 189 190 sfc_debug("Leaving RegularIntegralCG.gen_geometry_block") 191 return str(code)
192
193 - def gen_quadrature_runtime_block(self, data):
194 code = gen_const_token_definitions(self.itgrep.iter_runtime_quad_tokens(data)) 195 if code: 196 code = "// Geometric quantities and coefficients\n" + code 197 return code
198
199 - def gen_A_assignment_block(self, data):
200 "Assigning values to element tensor (for pre-integrated tensor components)" 201 # Only one symbolic integral supported: 202 assert data.integral is self.itgrep.symbolic_integral 203 204 if self.itgrep._on_facet: 205 # TODO: Each facet block can be extracted as a separate function 206 fr = self.itgrep.formrep 207 code = CodeFormatter() 208 code.begin_switch("facet") 209 for facet in range(fr.cell.num_facets): 210 A_tokens = self.itgrep.iter_A_tokens(data, facet) 211 code.begin_case(facet, braces=True) 212 code += "// Integrated element tensor entries" 213 code += gen_token_assignments(A_tokens) 214 code.end_case() 215 code.end_switch() 216 code = str(code) 217 218 else: 219 A_tokens = self.itgrep.iter_A_tokens(data) 220 code = gen_token_assignments(A_tokens) 221 code = "// Integrated element tensor entries\n{\n%s\n}" % code 222 223 return code
224
225 - def gen_quadrature_begin_block(self, data):
226 "Beginning of quadrature loop." 227 228 assert data.integral in self.itgrep.quadrature_integrals 229 230 code = CodeFormatter() 231 232 if self.itgrep._on_facet: 233 num_points = data.facet_quad_rules[0].num_points 234 code += "for(int iq=0; iq<%d; iq++)" % num_points 235 code += "{" 236 code.indent() 237 code += "// Fetch quadrature rule on facet" 238 code += "const double quad_weight = quadrature_weights[facet][iq];" 239 code += "const double *qp = quadrature_points[facet][iq];" 240 else: 241 num_points = data.quad_rule.num_points 242 code += "for(int iq=0; iq<%d; iq++)" % num_points 243 code += "{" 244 code.indent() 245 code += "// Fetch quadrature rule" 246 code += "const double quad_weight = quadrature_weights[iq];" 247 code += "const double *qp = quadrature_points[iq];" 248 249 nsd = self.itgrep.formrep.cell.nsd 250 gr = self.itgrep.formrep.geomrep 251 if 0: # TODO: Want to use xi[%d] for local coords 252 coords = ', '.join('qp[%d]' % i for i in range(nsd)) 253 code += "const double xi[%d] = { %s };" % (nsd, coords) 254 assert gr.xi_sym[0].printc() == "xi[0]" 255 else: 256 coords = ', '.join('%s=qp[%d]' % ("xyz"[i], i) for i in range(nsd)) 257 code += "const double %s;" % coords 258 assert gr.xi_sym[0].printc() == "x" 259 260 if 0: # TODO: Want to use x[%d] for global coords 261 coords = ', '.join(gr.x_expr[i] for i in range(nsd)) 262 code += "const double x[%d] = { %s };" % (nsd, coords) 263 assert gr.x_sym[0].printc() == "x[0]" 264 elif 0: # Yielding in iter_runtime_quad_tokens instead 265 coords = ', '.join("%s=%s" % (gr.x_sym[i].printc(), gr.x_expr[i].printc()) for i in range(nsd)) 266 code += "const double %s;" % (coords,) 267 assert gr.x_sym[0].printc() == "x0" 268 269 code.dedent() 270 271 return str(code)
272
273 - def gen_A_accumulation_block(self, data):
274 "Accumulation of nonzero element tensor entries." 275 assert data.integral in self.itgrep.quadrature_integrals 276 A_tokens = self.itgrep.iter_A_tokens(data) 277 code = gen_token_additions((A_sym, A_expr) 278 for (A_sym, A_expr) in A_tokens 279 if A_expr != 0) 280 assert code 281 code = "// Accumulating element tensor entries\n{\n%s\n}" % indent(code) 282 return code
283
284 - def gen_quadrature_end_block(self):
285 return "}"
286
288 if self.itgrep._symbol_counter > 0: 289 code = "double s[%d];" % self.itgrep._symbol_counter 290 else: 291 code = "" 292 return code
293
294 - def gen_tabulate_tensor(self):
295 "Generic function to generate tabulate_tensor composed by reusable code blocks." 296 # TODO: Improve code generation here for "robustness", several possibilities: 297 # - Use a loop over both trial and test functions. 298 # - Use a loop over trial functions, keep inlined application of test functions. 299 # - Apply "outlining" to reduce function size, that is, generate helper functions 300 # instead of inlining all expressions explicitly. 301 # For example one for each integral and one for each block. 302 303 sfc_debug("Entering RegularIntegralCG.tabulate_tensor") 304 305 fr = self.itgrep.formrep 306 fd = fr.formdata 307 r = fr.rank 308 309 def generate_partition_blocks(data, bf_dep_list, known): 310 pblocks = [] 311 known = frozenset(known) 312 313 # Define sequence of blocks to handle 314 todo = [] 315 for (iota, keep) in bf_dep_list: 316 deps = known 317 if keep is not None: 318 deps |= frozenset("v%d" % j for (j,k) in enumerate(keep) if k) 319 if deps in data.partitions: 320 todo.append((deps, iota)) 321 322 # Handle each block in sequence 323 for t in todo: 324 deps, iota = t 325 code = self.gen_partition_block(data, deps, iota) 326 if code: 327 pblocks += [code] 328 return pblocks
329 330 # Build a list of code blocks 331 blocks = [] 332 333 # Print debugging info about input arguments in code 334 blocks += [self.gen_pre_debug_code()] 335 336 # Geometry variables like G, Ginv, detG, n 337 blocks += [self.gen_geometry_block()] 338 339 # Symbolic integration part or element tensor reset 340 integral = self.itgrep.symbolic_integral 341 if integral is not None: 342 data = self.itgrep.integral_data[integral.measure()] 343 # TODO: Do we need any intermediate variable blocks? 344 # Maybe if we apply some optimizations. 345 # We could f.ex. compute the integrands to be 346 # symbolically integrated from the graph, 347 # with the same tokens precomputed as with quadrature. 348 blocks += [self.gen_A_assignment_block(data)] 349 else: 350 blocks += [gen_A_reset_block(len(self.itgrep.indices))] 351 352 # List of ways that expressions may depend on basis functions 353 bf_dep_list = [((), None)] 354 355 if r == 1: 356 m = fr.element_reps[fd.elements[0]].local_dimension 357 bf_dep_list += [((i,), (True,)) for i in range(m)] 358 elif r == 2: 359 m = fr.element_reps[fd.elements[0]].local_dimension 360 n = fr.element_reps[fd.elements[1]].local_dimension 361 bf_dep_list += [((None, i), (False, True)) for i in range(m)] 362 bf_dep_list += [((i, None), (True, False)) for i in range(m)] 363 bf_dep_list += [((i, j), (True, True)) for i in range(m) 364 for j in range(n)] 365 elif r > 2: 366 sfc_error("Support for higher order tensors not implemented.") 367 368 # Code for each quadrature integral 369 for integral in self.itgrep.quadrature_integrals: 370 # Data about this particular integral, including integration mode 371 # and computational graph 372 data = self.itgrep.integral_data[integral.measure()] 373 374 # FIXME: These should be precomputed, no dependencies 375 blocks += generate_partition_blocks(data, bf_dep_list, ()) 376 377 # Compute blocks of intermediate variables, s[...] = ...; 378 # that are independent of coordinates 379 blocks += generate_partition_blocks(data, bf_dep_list, ("c",)) 380 381 # Begin quadrature loop 382 blocks += [self.gen_quadrature_begin_block(data)] 383 384 # Build quadrature loop body from blocks 385 qblocks = [] 386 # Just for indentation similar to generated code :) 387 if True: 388 qblocks += [indent(self.gen_quadrature_runtime_block(data))] 389 390 # FIXME: These should be precomputed for each quadrature point, 391 # no runtime dependencies 392 qblocks += generate_partition_blocks(data, bf_dep_list, ("x",)) 393 394 # Compute blocks of intermediate variables, s[...] = ...; 395 qblocks += generate_partition_blocks(data, bf_dep_list, ("c", "x")) 396 397 # Accumulate element tensor values, A[...] += ...; 398 qblocks += [self.gen_A_accumulation_block(data)] 399 400 blocks += [indent(code) for code in qblocks] 401 402 # End quadrature loop 403 blocks += [self.gen_quadrature_end_block()] 404 405 # Print debugging info about output values in code 406 blocks += [self.gen_post_debug_code()] 407 408 # Insert declaration for allocation of symbols at beginning 409 code = self.gen_symbol_allocation_block() 410 if code: 411 blocks.insert(0, code) 412 413 # Compose the final tabulate_tensor code! 414 final_code = "\n\n".join(blocks) 415 416 sfc_debug("Leaving RegularIntegralCG.tabulate_tensor") 417 return final_code
418
419 - def gen_tabulate_tensor_quadrature(self):
420 return 'throw std::runtime_error("Not implemented.");' # FIXME
421 422
423 -class RegularCellIntegralCG(RegularIntegralCG):
424 - def __init__(self, itgrep):
425 RegularIntegralCG.__init__(self, itgrep) 426 # Generator trick: 427 # Iterators over members are used in gen_members and 428 # gen_constructor, avoid double iteration with tee: 429 430 # FIXME: Update this for multiple integral data, and 431 # for precomputation without constructor code 432 for integral in self.itgrep.quadrature_integrals: 433 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!") 434 data = self.itgrep.integral_data[integral.measure()] 435 436 a, b = tee(self.itgrep.iter_member_quad_tokens(data)) 437 self._iter_member_quad_tokens1 = a 438 self._iter_member_quad_tokens2 = b
439
440 - def header_template(self):
441 return ufc_utils.cell_integral_header
442
443 - def impl_template(self):
444 return ufc_utils.cell_integral_implementation
445
446 - def gen_members(self):
447 sfc_debug("entering CellIntegral.gen_members") 448 code = "" 449 450 # FIXME: Update this for multiple integral data, and 451 # for precomputation without constructor code 452 for integral in self.itgrep.quadrature_integrals: 453 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!") 454 data = self.itgrep.integral_data[integral.measure()] 455 456 #member_quad_tokens = self.itgrep.iter_member_quad_tokens() 457 member_quad_tokens = self._iter_member_quad_tokens1 458 member_quad_token_names = (str(s[0]).split(".")[-1] for s in member_quad_tokens) 459 decl = gen_symbol_declarations(member_quad_token_names) 460 if decl: 461 code += "\n\n" 462 code += "struct member_quad_tokens\n{\n" 463 code += indent(decl) 464 code += "\n};\n" 465 code += "member_quad_tokens qt[%d];\n" % data.quad_rule.num_points 466 467 sfc_debug("leaving CellIntegral.gen_members") 468 return code
469
470 - def gen_constructor(self):
471 sfc_debug("entering CellIntegral.gen_constructor") 472 code = "" 473 474 # FIXME: Update this for multiple integral data, and 475 # for precomputation without constructor code 476 for integral in self.itgrep.quadrature_integrals: 477 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!") 478 data = self.itgrep.integral_data[integral.measure()] 479 480 #member_quad_tokens = self.itgrep.iter_member_quad_tokens() 481 member_quad_tokens = self._iter_member_quad_tokens2 482 member_quad_tokens_assignments = gen_token_assignments(member_quad_tokens) 483 if member_quad_tokens_assignments: 484 # Construct quadrature loop 485 quad_loop_start_code = self.gen_quadrature_begin_block(data) 486 487 # Stitch code pieces together 488 code += "// Quad loop start code:\n" 489 code += quad_loop_start_code + "\n" 490 code += "// Member quad tokens:\n" 491 code += indent(member_quad_tokens_assignments) 492 code += "\n}\n" 493 494 sfc_debug("leaving CellIntegral.gen_constructor") 495 return code
496
497 - def gen_destructor(self):
498 return ""
499 500
501 -class RegularExteriorFacetIntegralCG(RegularIntegralCG):
502 - def __init__(self, itgrep):
503 RegularIntegralCG.__init__(self, itgrep) 504 # Generator trick: 505 # Iterators over members are used in gen_members and 506 # gen_constructor, avoid double iteration with tee: 507 508 # FIXME: Update this for multiple integral data, and for precomputation without constructor code 509 for integral in self.itgrep.quadrature_integrals: 510 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!") 511 data = self.itgrep.integral_data[integral.measure()] 512 513 a, b = tee(self.itgrep.iter_member_quad_tokens(data)) 514 self._iter_member_quad_tokens1 = a 515 self._iter_member_quad_tokens2 = b
516
517 - def header_template(self):
518 return ufc_utils.exterior_facet_integral_header
519
520 - def impl_template(self):
521 return ufc_utils.exterior_facet_integral_implementation
522
523 - def gen_members(self):
524 sfc_debug("entering ExteriorFacetIntegral.gen_members") 525 code = "" 526 527 # FIXME: Update this for multiple integral data, and for precomputation without constructor code 528 for integral in self.itgrep.quadrature_integrals: 529 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!") 530 data = self.itgrep.integral_data[integral.measure()] 531 532 member_quad_tokens = self._iter_member_quad_tokens1 #self.itgrep.iter_member_quad_tokens() 533 member_quad_token_names = (str(s[0]).split(".")[-1] for s in member_quad_tokens) 534 decl = gen_symbol_declarations(member_quad_token_names) 535 if decl: 536 code += "\n\n" 537 code += "struct member_quad_tokens\n{\n" 538 code += indent(decl) 539 code += "\n};\n" 540 code += "member_quad_tokens qt[%d][%d];\n" % (self.itgrep.formrep.cell.num_facets, data.facet_quad_rules[0].num_points) 541 542 sfc_debug("leaving ExteriorFacetIntegral.gen_members") 543 return code
544
545 - def gen_constructor(self):
546 sfc_debug("entering ExteriorFacetIntegral.gen_constructor") 547 code = "" 548 549 # FIXME: Update this for multiple integral data, and for precomputation without constructor code 550 for integral in self.itgrep.quadrature_integrals: 551 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!") 552 data = self.itgrep.integral_data[integral.measure()] 553 554 member_quad_tokens = self._iter_member_quad_tokens2 # itgrep.iter_member_quad_tokens() 555 member_quad_tokens_assignments = gen_token_assignments(member_quad_tokens) 556 if member_quad_tokens_assignments: 557 # Construct loop over facets with quad loop inside 558 facet_loop_begin = "for(int facet=0; facet<%d; facet++)\n{\n" % self.itgrep.formrep.cell.num_facets 559 quad_loop_start_code = self.gen_quadrature_begin_block(data) 560 561 # Stitch code pieces together 562 code += facet_loop_begin 563 code += quad_loop_start_code 564 code += indent(member_quad_tokens_assignments) 565 code += "\n}\n" 566 code += "\n}\n" 567 568 sfc_debug("leaving ExteriorFacetIntegral.gen_constructor") 569 return code
570
571 - def gen_destructor(self):
572 return ""
573 574
575 -class RegularInteriorFacetIntegralCG(RegularIntegralCG):
576 - def __init__(self, itgrep):
577 RegularIntegralCG.__init__(self, itgrep)
578
579 - def header_template(self):
580 return ufc_utils.exterior_facet_integral_header
581
582 - def impl_template(self):
583 return ufc_utils.exterior_facet_integral_implementation
584
585 - def gen_constructor(self):
586 return ""
587
588 - def gen_destructor(self):
589 return ""
590
591 - def gen_members(self):
592 return ""
593
595 return 'throw std::runtime_error("tabulate_tensor for interior facet integrals not implemented");'
596
598 return 'throw std::runtime_error("tabulate_tensor for interior facet integrals not implemented");'
599 600 601 # For dynamic importing of code generation strategy: 602 CellIntegralCG = RegularCellIntegralCG 603 ExteriorFacetIntegralCG = RegularExteriorFacetIntegralCG 604 InteriorFacetIntegralCG = RegularInteriorFacetIntegralCG 605