Home | Trees | Indices | Help |
|
---|
|
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 1618 """C++ formatter class for SFC code generation.""" 1922121 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"2931 CppDefaultFormatter.update_terminals(self, terminals) 32 self.facg.precompute(terminals, self.required)33 34 # --- Functions for generating definitions and other statements 35 3941 # 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)7375 # Coefficients in DG0 and R do not need to be computed, 76 # we can just refer directly to w[count][component]. 77 return ""7880 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)8587 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 code9496 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 code107109 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 code118120 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 code129131 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 code139141 # 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 code156 159 160 # --- Explicitly selecting rules for geometry here just for clarity: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 code174 175 # n[0] etc. 176 #facet_normal = CppDefaultFormatter.facet_normal178 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 code186 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: 200202 #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 code214216 #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
Home | Trees | Indices | Help |
|
---|
Generated by Epydoc 3.0.1 on Thu Jan 31 03:52:26 2013 | http://epydoc.sourceforge.net |