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

Source Code for Module sfc.codegeneration.sfc_uflacs_formatter

  1   
  2  from sfc.common.output import sfc_assert 
  3   
  4  from ufl.classes import (Terminal, Indexed, Grad, Restricted, 
  5                           Coefficient, Argument, 
  6                           SpatialCoordinate, 
  7                           FacetNormal, FacetArea, 
  8                           CellVolume, CellSurfaceArea, FacetArea, Circumradius) 
  9  from ufl.common import product, component_to_index 
 10  import uflacs 
 11  from uflacs.codeutils.target_formatter import (CppDefaultFormatter, 
 12                                                 analyse_terminalish) 
 13  from sfc.geometry.geometric_quantities import * 
 14  from sfc.codegeneration.geometrycg import GeometryCG 
 15  from sfc.codegeneration.formargumentcg import FormArgumentCG 
 16   
17 -class SfcUflacsFormatter(CppDefaultFormatter):
18 """C++ formatter class for SFC code generation.""" 19
20 - def __init__(self, itgrep, object_names=None):
21 CppDefaultFormatter.__init__(self, object_names) 22 self.itgrep = itgrep 23 self.geomcg = GeometryCG(itgrep.formrep.geomrep) 24 self.facg = FormArgumentCG(itgrep.formrep, on_facet=self.itgrep._on_facet) 25 self.form_argument_mapping = itgrep.formrep.formdata.function_replace_map 26 27 # Configuration options for statement formatting 28 self._output_basename = "A"
29
30 - def update_terminals(self, terminals):
31 CppDefaultFormatter.update_terminals(self, terminals) 32 self.facg.precompute(terminals, self.required)
33 34 # --- Functions for generating definitions and other statements 35
36 - def define_registers(self, num_registers):
37 return ['// Intermediate variables for subexpressions', 38 'double s[%d];' % num_registers]
39
41 # Default to no quantities 42 cellvolume = False 43 circumradius = False 44 cellsurfacearea = False 45 facetarea = False 46 facetnormal = False 47 48 # Set the ones needed 49 for o in self.terminals: 50 (t, ngrads, r) = analyse_terminalish(o) 51 if isinstance(t, FacetNormal): 52 assert ngrads == 0 53 facetnormal = True 54 if isinstance(t, CellVolume): 55 assert ngrads == 0 56 cellvolume = True 57 if isinstance(t, Circumradius): 58 assert ngrads == 0 59 circumradius = True 60 if isinstance(t, CellSurfaceArea): 61 assert ngrads == 0 62 cellsurfacearea = True 63 if isinstance(t, FacetArea): 64 assert ngrads == 0 65 facetarea = True 66 67 # Outsource to the GeometryCG class 68 return self.geomcg.generate_declarations(cellvolume=cellvolume, 69 circumradius=circumradius, 70 cellsurfacearea=cellsurfacearea, 71 facetarea=facetarea, 72 facetnormal=facetnormal)
73
75 # Coefficients in DG0 and R do not need to be computed, 76 # we can just refer directly to w[count][component]. 77 return ""
78
80 d = self.itgrep.formrep.formdata.cell.d 81 points = ((0.0,)*d,) 82 weigths = (1.0,) 83 return CppDefaultFormatter.define_quadrature_for_loop(self, points=points, 84 weights=weigths)
85
87 code = [] 88 89 # TODO: define num_quadrature_points elsewhere so tabulate_tensor_quadrature can just work 90 code += ["const int num_quadrature_points = %d;" % self.itgrep.num_quadrature_points] 91 92 code += ["for (int iq = 0; iq < num_quadrature_points; ++iq)"] 93 return code
94
96 if self.itgrep._on_facet: 97 code = [] 98 code += ["// Fetch quadrature point and weight on facet"] 99 code += ["const double quad_weight = quadrature_weights[facet][iq];"] 100 code += ["const double *qp = quadrature_points[facet][iq];"] 101 else: 102 code = [] 103 code += ["// Fetch quadrature point and weight in cell"] 104 code += ["const double quad_weight = quadrature_weights[iq];"] 105 code += ["const double *qp = quadrature_points[iq];"] 106 return code
107
109 code = [] 110 111 # FIXME: finalize D/facet_D vs quad_weight scaling 112 f = self.itgrep._on_facet 113 code += self.geomcg.generate_x_dependent_declarations(x=True, D=not f, facet_D=f) 114 115 # Compute x dependent geometry, if any? 116 # A spatially varying geometry mapping could be computed here. 117 return code
118
120 code = [] 121 for o in self.terminals: 122 (t, ngrads, r) = analyse_terminalish(o) 123 124 c = self.form_argument_mapping.get(t, t) 125 126 if isinstance(c, Coefficient) and not c.is_cellwise_constant(): 127 code += self.facg.generate_coefficient_declarations(c, ngrads, bool(r)) 128 return code
129
130 - def define_argument_for_loop(self, argument_count):
131 iname = "i%d" % (argument_count,) 132 133 fr = self.itgrep.formrep 134 dim = fr.element_reps[fr.formdata.elements[argument_count]].local_dimension 135 isize = str(dim) 136 137 code = ["for (int %s = 0; %s < %s; ++%s)" % (iname, iname, isize, iname)] 138 return code
139
140 - def define_argument_loop_vars(self, argument_count):
141 # Generate code for all needed derivatives of arguments with this count 142 code = [] 143 handled = set() 144 for o in self.terminals: 145 (t, ngrads, r) = analyse_terminalish(o) 146 147 c = self.form_argument_mapping.get(t, t) 148 149 if isinstance(c, Argument) and c.count() == argument_count: 150 # Avoid duplicates (restrictions will occur with + and -) 151 key = (c, ngrads, bool(r)) 152 if key not in handled: 153 handled.add(key) 154 code += self.facg.generate_argument_declarations(c, ngrads, bool(r)) 155 return code
156
157 - def output_variable_names(self, num_variables):
158 return ['%s[%d]' % (self._output_basename, i,) for i in xrange(num_variables)]
159 160 # --- Explicitly selecting rules for geometry here just for clarity:
161 - def spatial_coordinate(self, o, component=(), derivatives=(), restriction=None):
162 # Restriction has no effect on x 163 #postfix = {None:"", "+": "_p", "-": "_m"}[restriction] 164 165 sfc_assert(not derivatives, "I don't think this should happen.") 166 if derivatives: 167 sfc_warning("Compiler should be able to simplify derivatives of x!") 168 code = "1" if (len(derivatives) == 1 and i == derivatives[0]) else "0" 169 else: 170 i = component[0] if component else 0 171 code = self.geomcg.gen_x_names()[i] 172 self.require(o, component, derivatives, restriction, code) 173 return code
174 175 # n[0] etc. 176 #facet_normal = CppDefaultFormatter.facet_normal
177 - def facet_normal(self, o, component=(), derivatives=(), restriction=None):
178 sfc_assert(not derivatives, "I don't think this should happen.") 179 i = component[0] if component else 0 180 181 #postfix = {None:"", "+": "_p", "-": "_m"}[restriction] 182 code = self.geomcg.gen_facet_normal_names()[i] 183 184 self.require(o, component, derivatives, restriction, code) 185 return code
186 187 # K_vol 188 cell_volume = CppDefaultFormatter.cell_volume 189 190 # K_rad 191 circumradius = CppDefaultFormatter.circumradius 192 193 # K_surf 194 cell_surface_area = CppDefaultFormatter.cell_surface_area 195 196 # F_area 197 facet_area = CppDefaultFormatter.facet_area 198 199 # --- Rules for coefficients and arguments: 200
201 - def coefficient(self, o, component=(), derivatives=(), restriction=None):
202 #basename = self.coefficient_names[o] 203 o = self.form_argument_mapping.get(o, o) 204 205 if o.is_cellwise_constant(): 206 # For piecewise constant functions, refer directly to the dof symbol w[][] 207 assert not derivatives 208 code = self.facg.gen_constant_coefficient_name(o, component, restriction) 209 else: 210 code = self.facg.gen_coefficient_name(o, derivatives, component, restriction) 211 212 self.require(o, component, derivatives, restriction, code) 213 return code
214
215 - def argument(self, o, component=(), derivatives=(), restriction=None):
216 #basename = self.argument_names[o] 217 o = self.form_argument_mapping.get(o, o) 218 code = self.facg.gen_argument_name(o, derivatives, component, restriction) 219 self.require(o, component, derivatives, restriction, code) 220 return code
221