Home | Trees | Indices | Help |
|
---|
|
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 4722450 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 79 80 # --- Basis function and coefficient data for arguments 8183 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 v8991 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 s109111 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 w121 126 127 # --- Basis function and coefficient data for derivatives of arguments 128 132134 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))138140 """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 dv150152 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)167169 """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 Dv187189 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)202204 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 w218220 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)
Home | Trees | Indices | Help |
|
---|
Generated by Epydoc 3.0.1 on Thu Jan 31 03:52:25 2013 | http://epydoc.sourceforge.net |