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

Source Code for Module sfc.codegeneration.integralcg

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains shared code generation tools 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  from sfc.common.output import sfc_debug, sfc_error, sfc_assert, sfc_warning 
 28  from sfc.common.utilities import indices_subset, unique 
 29  from sfc.quadrature import gen_quadrature_rule_definition, gen_quadrature_rule_definitions 
 30  from sfc.codegeneration.codeformatting import indent, CodeFormatter, \ 
 31      gen_const_token_definitions, gen_token_prints, \ 
 32      gen_token_declarations, gen_token_assignments, gen_token_additions 
 33   
34 -class IntegralCG(object):
35 - def __init__(self, itgrep):
36 self.itgrep = itgrep 37 self.options = self.itgrep.formrep.options.code.integral
38
39 - def hincludes(self):
40 l = [] 41 return l
42
43 - def cincludes(self):
44 l = [] 45 return l
46
47 - def generate_code_dict(self):
48 vars = { 49 "classname" : self.itgrep.classname, 50 "constructor" : indent(self.gen_constructor()), 51 "constructor_arguments" : indent(self.gen_constructor_arguments()), 52 "initializer_list" : indent(self.gen_initializer_list()), 53 "destructor" : indent(self.gen_destructor()), 54 "members" : indent(self.gen_members()), 55 "tabulate_tensor" : indent(self.gen_tabulate_tensor()), 56 "tabulate_tensor_quadrature" : indent(self.gen_tabulate_tensor_quadrature()), 57 } 58 return vars
59
60 - def generate_support_code(self):
61 ccode = "" 62 63 # Write quadrature rule if we have one: # TODO: Place quadrature rules in a separate shared file? 64 fr = self.itgrep.formrep 65 if fr.quad_rule is not None: 66 ccode += gen_quadrature_rule_definition(fr.quad_rule) 67 68 if fr.facet_quad_rules is not None: 69 ccode += gen_quadrature_rule_definitions(fr.facet_quad_rules) 70 71 # Write sign correction code if needed: TODO: Update sign code 72 #nsd = fr.cell.nsd 73 #handled_fe_names = set() 74 #fe_list = fr.fe_list 75 #for fe in fe_list: 76 # if fe_is_signed(fe): 77 # fe_name = strings.finite_element_classname(fe) 78 # if not fe_name in handled_fe_names: 79 # handled_fe_names.add(fe_name) 80 # ccode += tabulate_sign_code(fe, nsd) 81 82 return ccode
83
84 - def gen_constructor(self):
85 raise NotImplementedError
86
88 return ""
89
90 - def gen_initializer_list(self):
91 return ""
92
93 - def gen_destructor(self):
94 raise NotImplementedError
95
96 - def gen_members(self):
97 raise NotImplementedError
98
99 - def gen_pre_debug_code(self):
100 if not self.options.enable_debug_code: 101 return "" 102 code = CodeFormatter() 103 code.begin_debug() 104 code.begin_block() 105 code += 'std::cout << std::endl;' 106 code += 'std::cout << "SFC DEBUGGING OUTPUT: " << std::endl;' 107 code += 'std::cout << "void %s::tabulate_tensor(...)" << std::endl;' % self.itgrep.classname 108 code += 'std::cout << "{" << std::endl;' 109 fr = self.itgrep.formrep 110 fd = fr.formdata 111 for k in range(fr.num_coefficients): 112 element = fd.elements[k] 113 rep = fr.element_reps[element] 114 code += 'std::cout << " w[%d][:] = [ ";' % k 115 code += "for(int j=0; j<%d; j++)" % rep.local_dimension 116 code.begin_block() 117 code += 'std::cout << w[%d][j] << ", ";' % k 118 code.end_block() 119 code += 'std::cout << " ]" << std::endl;' 120 code += 'std::cout << " now computing element tensor..." << std::endl;' 121 code.end_block() 122 code.end_debug() 123 return str(code)
124
125 - def gen_post_debug_code(self):
126 if not self.options.enable_debug_code: 127 return "" 128 code = CodeFormatter() 129 code.begin_debug() 130 code.begin_block() 131 code += 'std::cout << " ... done computing element tensor." << std::endl;' 132 code += "for(int k=0; k<%d; k++)" % len(self.itgrep.indices) 133 code.begin_block() 134 # TODO: enhance this output by printing indices 135 code += 'std::cout << " A[" << k << "]" << " = " << A[k] << std::endl;' 136 code.end_block() 137 code += 'std::cout << "}" << std::endl;' 138 code.end_block() 139 code.end_debug() 140 return str(code)
141
142 - def gen_partition_block(self, data, deps, basis_functions):
143 "Generate code for a partition of the linearized computational graph." 144 code = gen_token_assignments( self.itgrep.iter_partition(data, deps, basis_functions) ) 145 if code: 146 deps = ", ".join(deps) 147 if basis_functions: 148 deps += "; using basis functions %s" % str(basis_functions) 149 code = "// Partition depending on %s\n{\n%s\n}" % (deps, indent(code)) 150 return code
151
152 - def gen_geometry_block(self):
153 sfc_debug("Entering IntegralCG.gen_geometry_block") 154 155 code = CodeFormatter() 156 157 # --- Sign correction code: # TODO 158 #sign_sym = itgrep.sign_sym 159 #for fe_name in sign_sym.keys(): 160 # nbf = len(sign_sym[fe_name]) 161 # code += use_sign_code2(fe_name, nbf, nsd) 162 163 # --- Generate code for for initial geometry tokens 164 # (vertex coordinates, affine map, normal vector, ...) 165 geometry_tokens = self.itgrep.iter_geometry_tokens() 166 167 # if we want to print debug tokens, duplicate the iterator 168 if self.options.enable_debug_code: 169 geometry_tokens, debug_geometry_tokens = tee(geometry_tokens) 170 171 code += "// Geometric quantities" 172 code += gen_const_token_definitions(geometry_tokens) 173 174 if self.options.enable_debug_code: 175 code += "" 176 code.begin_debug() 177 code += gen_token_prints(debug_geometry_tokens) 178 code.end_debug() 179 code += "" 180 181 # --- Facet tokens code: 182 if self.itgrep._on_facet: 183 fr = self.itgrep.formrep 184 code += "// Geometric quantities on each facet" 185 for facet in range(fr.cell.num_facets): 186 facet_tokens = self.itgrep.iter_facet_tokens(facet) 187 188 if facet == 0: 189 # if we want to print debug tokens, duplicate the iterator 190 if self.options.enable_debug_code: 191 facet_tokens, debug_facet_tokens = tee(facet_tokens) 192 193 # duplicate the iterator to generate declarations outside switch 194 facet_tokens, facet_tokens2 = tee(facet_tokens) 195 code += gen_token_declarations(facet_tokens2) 196 code.begin_switch("facet") 197 198 code.begin_case(facet, braces=True) 199 code += "// Geometric quantities on facet %d" % facet 200 code += gen_token_assignments(facet_tokens) 201 code.end_case() 202 code += "default:" 203 code.indent() 204 code += 'throw std::runtime_error("Invalid facet number.");' 205 code.dedent() 206 code.end_switch() 207 208 if self.options.enable_debug_code: 209 code += "" 210 code.begin_debug() 211 code += 'std::cout << "facet = " << facet << std::endl;' 212 code += gen_token_prints(debug_facet_tokens) 213 code.end_debug() 214 code += "" 215 216 sfc_debug("Leaving IntegralCG.gen_geometry_block") 217 return str(code)
218
219 - def gen_quadrature_runtime_block(self, data):
220 code = gen_const_token_definitions(self.itgrep.iter_runtime_quad_tokens(data)) 221 if code: 222 code = "// Geometric quantities and coefficients\n" + code 223 return code
224
225 - def gen_A_assignment_block(self, data):
226 "Assigning values to element tensor (for pre-integrated tensor components)" 227 # Only one symbolic integral supported: 228 assert data.integral is self.itgrep.symbolic_integral 229 230 if self.itgrep._on_facet: 231 # TODO: Each facet block can be extracted as a separate function 232 fr = self.itgrep.formrep 233 code = CodeFormatter() 234 code.begin_switch("facet") 235 for facet in range(fr.cell.num_facets): 236 A_tokens = self.itgrep.iter_A_tokens(data, facet) 237 code.begin_case(facet, braces=True) 238 code += "// Integrated element tensor entries" 239 code += gen_token_assignments(A_tokens) 240 code.end_case() 241 code.end_switch() 242 code = str(code) 243 244 else: 245 A_tokens = self.itgrep.iter_A_tokens(data) 246 code = gen_token_assignments(A_tokens) 247 code = "// Integrated element tensor entries\n{\n%s\n}" % code 248 249 return code
250
251 - def gen_A_reset_block(self):
252 # Zero element tensor 253 code = "// Reset element tensor\n" 254 code += "memset(A, 0, sizeof(double)*%d);" % len(self.itgrep.indices) 255 return code
256
257 - def gen_quadrature_begin_block(self, data):
258 "Beginning of quadrature loop." 259 260 assert data.integral in self.itgrep.quadrature_integrals 261 262 code = CodeFormatter() 263 264 if self.itgrep._on_facet: 265 num_points = data.facet_quad_rules[0].num_points 266 code += "for(int iq=0; iq<%d; iq++)" % num_points 267 code += "{" 268 code.indent() 269 code += "// Fetch quadrature rule" 270 code += "const double quad_weight = facet_quad_weights[facet][iq];" 271 code += "const double *_p = facet_quad_points[facet][iq];" 272 else: 273 num_points = data.quad_rule.num_points 274 code += "for(int iq=0; iq<%d; iq++)" % num_points 275 code += "{" 276 code.indent() 277 code += "// Fetch quadrature rule" 278 code += "const double quad_weight = quad_weights[iq];" 279 code += "const double *_p = quad_points[iq];" 280 281 nsd = self.itgrep.formrep.cell.nsd 282 if nsd > 0: code += "const double x = _p[0];" 283 if nsd > 1: code += "const double y = _p[1];" 284 if nsd > 2: code += "const double z = _p[2];" 285 286 code.dedent() 287 288 return str(code)
289
290 - def gen_A_accumulation_block(self, data):
291 "Accumulation of nonzero element tensor entries." 292 assert data.integral in self.itgrep.quadrature_integrals 293 A_tokens = self.itgrep.iter_A_tokens(data) 294 code = gen_token_additions((A_sym, A_expr) for (A_sym, A_expr) in A_tokens if A_expr != 0) 295 assert code 296 code = "// Accumulating element tensor entries\n{\n%s\n}" % indent(code) 297 return code
298
299 - def gen_quadrature_end_block(self):
300 return "}"
301
303 if self.itgrep._symbol_counter > 0: 304 code = "double s[%d];" % self.itgrep._symbol_counter 305 else: 306 code = "" 307 return code
308
309 - def gen_tabulate_tensor(self):
310 "Generic function to generate tabulate_tensor composed by reusable code blocks." 311 # TODO: Improve code generation here for "robustness", several possibilities: 312 # - Use a loop over both trial and test functions. 313 # - Use a loop over trial functions, keep inlined application of test functions. 314 # - Apply "outlining" to reduce function size, that is, generate helper functions 315 # instead of inlining all expressions explicitly. 316 # For example one for each integral and one for each block. 317 318 sfc_debug("Entering IntegralCG.tabulate_tensor") 319 320 fr = self.itgrep.formrep 321 fd = fr.formdata 322 r = fr.rank 323 324 def generate_partition_blocks(data, bf_dep_list, known): 325 pblocks = [] 326 known = frozenset(known) 327 328 # Define sequence of blocks to handle 329 todo = [] 330 for (iota, keep) in bf_dep_list: 331 deps = known 332 if keep is not None: 333 deps |= frozenset("v%d" % j for (j,k) in enumerate(keep) if k) 334 if deps in data.partitions: 335 todo.append((deps, iota)) 336 337 # Handle each block in sequence 338 for t in todo: 339 deps, iota = t 340 code = self.gen_partition_block(data, deps, iota) 341 if code: 342 pblocks += [code] 343 return pblocks
344 345 # Build a list of code blocks 346 blocks = [] 347 348 # Print debugging info about input arguments in code 349 blocks += [self.gen_pre_debug_code()] 350 351 # Geometry variables like G, Ginv, detG, n 352 blocks += [self.gen_geometry_block()] 353 354 # Symbolic integration part or element tensor reset 355 integral = self.itgrep.symbolic_integral 356 if integral is not None: 357 data = self.itgrep.integral_data[integral.measure()] 358 # TODO: Do we need any intermediate variable blocks? 359 # Maybe if we apply some optimizations. 360 # We could f.ex. compute the integrands to be 361 # symbolically integrated from the graph, 362 # with the same tokens precomputed as with quadrature. 363 blocks += [self.gen_A_assignment_block(data)] 364 else: 365 blocks += [self.gen_A_reset_block()] 366 367 # List of ways that expressions may depend on basis functions 368 bf_dep_list = [((), None)] 369 370 if r == 1: 371 m = fr.element_reps[fd.elements[0]].local_dimension 372 bf_dep_list += [((i,), (True,)) for i in range(m)] 373 elif r == 2: 374 m = fr.element_reps[fd.elements[0]].local_dimension 375 n = fr.element_reps[fd.elements[1]].local_dimension 376 bf_dep_list += [((None, i), (False, True)) for i in range(m)] 377 bf_dep_list += [((i, None), (True, False)) for i in range(m)] 378 bf_dep_list += [((i, j), (True, True)) for i in range(m) for j in range(n)] 379 elif r > 2: 380 sfc_error("Support for higher order tensors not implemented.") 381 382 # Code for each quadrature integral 383 for integral in self.itgrep.quadrature_integrals: 384 # Data about this particular integral, including integration mode and computational graph 385 data = self.itgrep.integral_data[integral.measure()] 386 387 # FIXME: These should be precomputed, no dependencies 388 blocks += generate_partition_blocks(data, bf_dep_list, ()) 389 390 # Compute blocks of intermediate variables, s[...] = ...; that are independent of coordinates 391 blocks += generate_partition_blocks(data, bf_dep_list, ("c",)) 392 393 # Begin quadrature loop 394 blocks += [self.gen_quadrature_begin_block(data)] 395 396 # Build quadrature loop body from blocks 397 qblocks = [] 398 # Just for indentation similar to generated code :) 399 if True: 400 qblocks += [indent(self.gen_quadrature_runtime_block(data))] 401 402 # FIXME: These should be precomputed for each quadrature point, no runtime dependencies 403 qblocks += generate_partition_blocks(data, bf_dep_list, ("x",)) 404 405 # Compute blocks of intermediate variables, s[...] = ...; 406 qblocks += generate_partition_blocks(data, bf_dep_list, ("c", "x")) 407 408 # Accumulate element tensor values, A[...] += ...; 409 qblocks += [self.gen_A_accumulation_block(data)] 410 411 blocks += [indent(code) for code in qblocks] 412 413 # End quadrature loop 414 blocks += [self.gen_quadrature_end_block()] 415 416 # Print debugging info about output values in code 417 blocks += [self.gen_post_debug_code()] 418 419 # Insert declaration for allocation of symbols at beginning 420 code = self.gen_symbol_allocation_block() 421 if code: 422 blocks.insert(0, code) 423 424 # Compose the final tabulate_tensor code! 425 final_code = "\n\n".join(blocks) 426 427 sfc_debug("Leaving IntegralCG.tabulate_tensor") 428 return final_code
429
430 - def gen_tabulate_tensor_quadrature(self):
431 return 'throw std::runtime_error("Not implemented.");' # FIXME
432