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

Source Code for Module sfc.codegeneration.geometrycg

  1   
  2  from sfc.representation.geometryrepresentation import GeometryRepresentation 
  3   
4 -class GeometryCG(object):
5 - def __init__(self, geomrep):
6 self.geomrep = geomrep
7
8 - def generate_declarations(self, 9 cellvolume, 10 circumradius, 11 cellsurfacearea, 12 facetarea, 13 facetnormal):
14 # FIXME: Skip what's not needed 15 code = [] 16 code += ['', '// Vertex coordinates in more compact notation'] 17 code += self.gen_vx_code() 18 19 code += ['', '// Offset for affine geometry mapping'] 20 code += self.gen_x0_code() 21 22 code += ['', '// Jacobi of geometry mapping'] 23 code += self.gen_G_code() 24 25 code += ['', '// Reference cell volume ratio'] 26 code += self.gen_detGtmp_code() 27 code += self.gen_detG_code() 28 code += self.gen_detG_sign_code() 29 30 code += ['', '// Inverted Jacobi of geometry mapping'] 31 code += self.gen_Ginv_code() 32 33 code += ['', '// Cell derived UFL quantities'] 34 if cellvolume: 35 code += self.gen_cell_volume_code() 36 if circumradius: 37 code += self.gen_circumradius_code() 38 if cellsurfacearea: 39 code += self.gen_cell_surface_code() 40 if facetarea: 41 code += self.gen_facet_area_code() 42 if facetnormal: 43 code += self.gen_facet_normal_code() 44 45 return code
46
47 - def generate_x_dependent_declarations(self, x, D, facet_D):
48 gr = self.geomrep 49 nsd = gr.sfc_cell.nsd 50 51 code = [] 52 code += ['', '// Local coordinates'] 53 code += self.gen_xi_code() 54 55 if x: 56 code += ['', '// Global coordinates'] 57 code += self.gen_x_code() 58 59 assert not (D and facet_D) 60 if D or facet_D: 61 code += ['', '// Integration point scaling factor'] 62 if D: 63 code += self.gen_D_code() 64 if facet_D: 65 code += self.gen_facet_D_code() 66 67 return code
68
69 - def gen_xi_code(self):
70 code = [] 71 gr = self.geomrep 72 nsd = gr.sfc_cell.nsd 73 if 0: 74 coords = ', '.join('qp[%d]' % i for i in range(nsd)) 75 code += ["const double xi[%d] = { %s };" % (nsd, coords)] 76 assert gr.xi_sym[0].printc() == "xi[0]" 77 else: 78 code += ["const double %s = qp[%d];" % (gr.xi_sym[i].printc(), i) 79 for i in range(nsd)] 80 assert gr.xi_sym[0].printc() == "x" 81 return code
82
83 - def gen_x_names(self):
84 gr = self.geomrep 85 nsd = gr.sfc_cell.nsd 86 code = [gr.x_sym[i].printc() for i in range(nsd)] 87 return code
88 - def gen_x_code(self):
89 gr = self.geomrep 90 nsd = gr.sfc_cell.nsd 91 code = [] 92 93 if 0: 94 coords = ', '.join('qp[%d]' % i for i in range(nsd)) 95 code += ["const double x[%d] = { %s };" % (nsd, coords)] 96 assert gr.x_sym[0].printc() == "x[0]" 97 else: 98 code += ["const double %s = %s;" % (gr.x_sym[i].printc(), gr.x_expr[i].printc()) 99 for i in range(nsd)] 100 assert gr.x_sym[0].printc() == "x0" 101 return code
102
103 - def gen_vx_names(self):
104 "Return names of vertex coordinate " 105 sc = self.geomrep.sfc_cell 106 d = sc.nsd 107 num_vertices = sc.num_vertices 108 return [gr.vx_sym[i][j].printc() 109 for i in range(nv) 110 for j in range(d)]
111 - def gen_vx_code(self):
112 gr = self.geomrep 113 sc = gr.sfc_cell 114 d = sc.nsd 115 nv = sc.num_vertices 116 code = ["const double %s = %s;" % (gr.vx_sym[i][j].printc(), 117 gr.vx_expr[i][j].printc()) 118 for i in range(nv) 119 for j in range(d)] 120 return code
121
122 - def gen_x0_names(self):
123 d = self.geomrep.sfc_cell.nsd 124 return [self.geomrep.x0_sym[i].printc() % i for i in range(d)]
125 - def gen_x0_code(self):
126 # x0 just just set equal to vertex 0 127 return []
128 - def xgen_x0_code(self):
129 d = self.geomrep.sfc_cell.nsd 130 expr = ", ".join(self.geomrep.x0_expr[i].printc() for i in range(d)) 131 code = "const double x0[%d] = { %s };" % (d, expr) 132 return [code]
133
134 - def gen_G_names(self):
135 d = self.geomrep.sfc_cell.nsd 136 return [self.geomrep.G_sym[i,j].printc() 137 for i in range(d) 138 for j in range(d)]
139 - def gen_G_code(self):
140 gr = self.geomrep 141 d = gr.sfc_cell.nsd 142 if 0: 143 expr = ", ".join(gr.G_expr[i,j].printc() 144 for i in range(d) 145 for j in range(d)) 146 code = ["const double G[%d] = { %s };" % (d*d, expr)] 147 else: 148 code = ["const double %s = %s;" % (gr.G_sym[i,j].printc(), 149 gr.G_expr[i,j].printc()) 150 for i in range(d) for j in range(d)] 151 return code
152
153 - def gen_Ginv_names(self):
154 d = self.geomrep.sfc_cell.nsd 155 return [self.geomrep.Ginv_sym[i,j].printc() 156 for i in range(d) 157 for j in range(d)]
158 - def gen_Ginv_code(self):
159 gr = self.geomrep 160 d = gr.sfc_cell.nsd 161 if 0: 162 expr = ", ".join(gr.Ginv_expr[i,j].printc() 163 for i in range(d) 164 for j in range(d)) 165 code = ["const double Ginv[%d] = { %s };" % (d*d, expr)] 166 else: 167 code = ["const double %s = %s;" % (gr.Ginv_sym[i,j].printc(), 168 gr.Ginv_expr[i,j].printc()) 169 for i in range(d) for j in range(d)] 170 return code
171
172 - def gen_detG_names(self):
173 return [self.geomrep.detG_sym.printc()]
174 - def gen_detG_code(self):
175 expr = self.geomrep.detG_expr.printc() 176 code = "const double detG = %s;" % expr 177 return [code]
178
179 - def gen_detGtmp_names(self):
180 return [self.geomrep.detGtmp_sym.printc()]
181 - def gen_detGtmp_code(self):
182 expr = self.geomrep.detGtmp_expr.printc() 183 code = "const double detGtmp = %s;" % expr 184 return [code]
185
186 - def gen_detG_sign_names(self):
187 return [self.geomrep.detG_sign_sym.printc()]
188 - def gen_detG_sign_code(self):
189 expr = self.geomrep.detG_sign_expr.printc() 190 code = "const double detG_sign = %s;" % expr 191 return [code]
192
193 - def gen_facet_D_names(self):
194 return ["facet_D"]
195 - def gen_facet_D_code(self):
196 gr = self.geomrep 197 code = [] 198 code += ["switch (facet)", "{"] 199 for f in range(gr.sfc_cell.num_facets): 200 code += ["case %d:" % f] 201 code += ["%s = %s;" % (gr.facet_D_sym.printc(), 202 gr.facet_D_expr[f].printc())] 203 code += ["}"] 204 # D = facet scaling times quadrature weight 205 code += ["const double %s = %s * %s;" % (gr.D_sym.printc(), 206 gr.facet_D_sym.printc(), 207 gr.quad_weight_sym.printc())] 208 return code
209
210 - def gen_D_names(self):
211 return ["D"]
212 - def gen_D_code(self):
213 gr = self.geomrep 214 code = [] 215 # FIXME: Consistent quad_weight scaling! * gr.quad_weight_sym? 216 code += ["const double %s = %s;" % (gr.D_sym.printc(), gr.D_expr.printc())] 217 return code
218
219 - def gen_facet_normal_names(self):
220 d = self.geomrep.sfc_cell.nsd 221 return [self.geomrep.n_sym[0].printc() % i for i in range(d)]
222 - def gen_facet_normal_code(self):
223 gr = self.geomrep 224 d = gr.sfc_cell.nsd 225 code = [] 226 assert self.geomrep.n_sym[0].printc() == "n0" 227 code = [] 228 code += ["double %s;" % (', '.join('n%d'%i for i in range(d)),)] 229 code += ["switch (facet)", "{"] 230 for f in range(gr.sfc_cell.num_facets): 231 code += ["case %d:" % f] 232 # Mapping local facet normal to global cell orientation: 233 for i in range(d): 234 name = gr.n_tmp_sym[i].printc() 235 expr = gr.n_tmp_expr[f][i].printc() 236 code += [" %s = %s;" % (name, expr)] 237 code += ["}"] 238 # Normalize global facet normal: 239 #normscale = self.detG_sign_sym.printc(), # FIXME: scale by sign? 240 normscale = "1.0" 241 norm2 = ' + '.join((gr.n_sym[i].printc() + '*' + gr.n_sym[i].printc()) for i in range(d)) 242 code += ['const double nnorm = %s * sqrt(%s);' % (normscale, norm2)] 243 for i in range(d): 244 code += ['%s /= nnorm;' % (gr.n_sym[i].printc(),)] 245 return [code]
246
247 - def gen_cell_volume_names(self):
248 return ["K_vol"]
249 - def gen_cell_volume_code(self):
250 gr = self.geomrep 251 sc = gr.sfc_cell 252 expr = (gr.detG_sym * sc.reference_volume).printc() # FIXME: Validate this 253 code = "const double K_vol = %s;" % expr 254 return [code]
255
256 - def gen_circumradius_names(self):
257 return ["K_rad"]
258 - def gen_circumradius_code(self):
259 expr = "FIXME" 260 code = "const double K_rad = %s;" % expr 261 return [code]
262
263 - def gen_cell_surface_names(self):
264 return ["K_surf"]
265 - def gen_cell_surface_code(self):
266 expr = "FIXME" 267 code = "const double K_surf = %s;" % expr 268 return [code]
269
270 - def gen_facet_area_names(self):
271 return ["F_area"]
272 - def gen_facet_area_code(self):
273 expr = "FIXME" 274 code = "const double F_area = %s;" % expr 275 return [code]
276