Package sfc :: Package representation :: Module formargumentrepresentation
[hide private]
[frames] | no frames]

Source Code for Module sfc.representation.formargumentrepresentation

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains a representation class for forms, storing information about all form aspects. 
  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-04-19 
 25   
 26  import copy 
 27  import math 
 28   
 29  import swiginac 
 30  import SyFi 
 31   
 32  import ufl 
 33  from ufl.algorithms import extract_max_quadrature_element_degree,\ 
 34      estimate_quadrature_degree, estimate_total_polynomial_degree 
 35   
 36  from sfc.common import sfc_assert 
 37  from sfc.common.names import finite_element_classname, dof_map_classname, \ 
 38                             integral_classname, form_classname 
 39  from sfc.common.utilities import (matrix_to_list, list_to_vector, 
 40                                    index_string, dot_product, 
 41                                    build_count_dict) 
 42  from sfc.geometry import geometry_mapping 
 43  from sfc.quadrature import find_quadrature_rule 
 44  from sfc.symbolic_utils import symbolic_matrix, symbolic_vector, symbol, symbols, cross, sub, inner 
 45  from sfc.representation.elementrepresentation import create_syfi_polygon 
 46  from sfc.representation.geometryrepresentation import GeometryRepresentation 
 47   
48 -class FormArgumentRepresentation(object):
49 - def __init__(self, formrep):
50 self.formrep = formrep 51 52 # Storing some references here for convenience. 53 # (Note that this class stems from a partial refactoring of old FormRepresentation class) 54 self.formdata = formrep.formdata 55 self.element_reps = formrep.element_reps 56 self.geomrep = formrep.geomrep 57 self.rank = formrep.rank 58 self.num_coefficients = formrep.num_coefficients 59 60 # ... Build symbols for all coefficient dofs as UFC input variables w[][] 61 # Get dimensions from element representations in coefficient ordering 62 dims = [self.element_reps[self.formdata.elements[self.rank+i]].local_dimension 63 for i in range(self.num_coefficients)] 64 # Build names 65 self.w_dof_names = [] 66 for i, dim in enumerate(dims): 67 wi_dof_names = ["w[%d][%d]" % (i, j) for j in range(dim)] 68 self.w_dof_names.append(wi_dof_names) 69 # Build symbols 70 self.w_dofs = [symbols(wdn) for wdn in self.w_dof_names]
71 72 # ... TODO: Build dofindex-templated names for all references to precomputed basis functions? 73
74 - def __str__(self):
75 s = "" 76 s += " Argument tokens:\n" 77 s += " self.w_dofs = %s\n" % self.w_dofs 78 return s
79 80 # --- Basis function and coefficient data for arguments 81
82 - def v_expr(self, iarg, i, component):
83 elm = self.formdata.elements[iarg] 84 #component, elm = elm.extract_component(component) 85 86 rep = self.element_reps[elm] 87 v = rep.basis_function(i, component) 88 return v
89
90 - def v_sym(self, iarg, i, component, on_facet):
91 elm = self.formdata.elements[iarg] 92 93 e = self.v_expr(iarg, i, component) 94 e2 = e.evalf() 95 if e2.nops() == 0: 96 # FIXME: Remove this line after fixing code generation, 97 # this leads to "double 0;" statements: 98 if e2 == 0: 99 return e2 100 101 scomponent, selm = elm.extract_component(component) 102 num = self.formrep.element_number[selm] 103 104 prefix = "qt[facet][iq]." if on_facet else "qt[iq]." 105 suffix = index_string((num, i) + scomponent) 106 name = "%sv_%s" % (prefix, suffix) 107 s = symbol(name) 108 return s
109
110 - def w_expr(self, iarg, component, use_symbols, on_facet):
111 elm = self.formdata.elements[self.rank+iarg] 112 rep = self.element_reps[elm] 113 dim = rep.local_dimension 114 if use_symbols: 115 vs = [self.v_sym(iarg+self.rank, i, component, on_facet=on_facet) 116 for i in range(dim)] 117 else: 118 vs = [self.v_expr(iarg+self.rank, i, component) for i in range(dim)] 119 w = dot_product(self.w_dofs[iarg], vs) 120 return w
121
122 - def w_sym(self, iarg, component):
123 name = "w_%s" % index_string((iarg,) + component) 124 s = symbol(name) 125 return s
126 127 # --- Basis function and coefficient data for derivatives of arguments 128
129 - def ddxi(self, f, i):
130 gr = self.formrep.geomrep 131 return swiginac.diff(f, gr.xi_sym[i])
132
133 - def ddx(self, f, i):
134 gr = self.formrep.geomrep 135 # TODO: Verify this line (in particular transposed or not!) 136 return sum(gr.Ginv_sym[j,i]*self.ddxi(f, j) 137 for j in range(gr.sfc_cell.nsd))
138
139 - def dv_expr(self, iarg, i, component, d):
140 """Return expression for dv/dxi, with v being a particular 141 component of basis function i in argument space iarg, 142 and d is a tuple (i.e. multiindex) of xi directions (local coordinates).""" 143 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 144 d = tuple(sorted(d)) 145 v = self.v_expr(iarg, i, component) 146 dv = v 147 for k in d: 148 dv = self.ddxi(dv, k) 149 return dv
150
151 - def dv_sym(self, iarg, i, component, d, on_facet):
152 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 153 d = tuple(sorted(d)) 154 155 e = self.dv_expr(iarg, i, component, d) 156 e2 = e.evalf() 157 if e2.nops() == 0: 158 # FIXME: Remove this line after fixing code generation, 159 # this leads to "double 0;" statements: 160 if e2 == 0: 161 return e2 162 163 prefix = "qt[facet][iq]." if on_facet else "qt[iq]." 164 suffix = index_string((iarg, i) + component + d) 165 name = "%sdv_dxi_%s" % (prefix, suffix) 166 return symbol(name)
167
168 - def Dv_expr(self, iarg, i, component, d, use_symbols, on_facet):
169 """Return expression for dv/dx, with v being a particular 170 component of basis function i in argument space iarg, 171 and d is a tuple (i.e. multiindex) of x directions (global coordinates).""" 172 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 173 d = tuple(sorted(d)) 174 # Get local derivatives in all directions 175 gr = self.formrep.geomrep 176 if use_symbols: 177 dv = [self.dv_sym(iarg, i, component, (j,), on_facet=on_facet) 178 for j in range(gr.sfc_cell.nsd)] 179 else: 180 dv = [self.dv_expr(iarg, i, component, (j,)) 181 for j in range(gr.sfc_cell.nsd)] 182 # Apply mapping to dv 183 j, = d 184 # TODO: Verify this line (in particular transposed or not!) 185 Dv = sum(gr.Ginv_sym[k, j]*dv[k] for k in range(gr.sfc_cell.nsd)) 186 return Dv
187
188 - def Dv_sym(self, iarg, i, component, d, on_facet):
189 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 190 d = tuple(sorted(d)) 191 192 e = self.Dv_expr(iarg, i, component, d, False, on_facet=on_facet) 193 e2 = e.evalf() 194 if e2.nops() == 0: 195 # FIXME: Remove this line after fixing code generation, 196 # this leads to "double 0;" statements: 197 if e2 == 0: 198 return e2 199 200 name = "Dv_%s" % index_string((iarg, i) + component + d) 201 return symbol(name)
202
203 - def Dw_expr(self, iarg, component, d, use_symbols, on_facet):
204 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 205 d = tuple(sorted(d)) 206 elm = self.formdata.elements[self.rank+iarg] 207 rep = self.formrep.element_reps[elm] 208 dim = rep.local_dimension 209 if use_symbols: 210 vs = [self.Dv_sym(self.rank+iarg, i, component, d, on_facet=on_facet) 211 for i in range(dim)] 212 else: 213 vs = [self.Dv_expr(self.rank+iarg, i, component, d, 214 use_symbols=False, on_facet=on_facet) 215 for i in range(dim)] 216 w = dot_product(self.w_dofs[iarg], vs) 217 return w
218
219 - def Dw_sym(self, iarg, component, d):
220 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 221 d = tuple(sorted(d)) 222 name = "Dw_%s" % index_string((iarg,) + component + d) 223 return symbol(name)
224