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

Source Code for Module sfc.representation.formrepresentation

  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  from sfc.representation.formargumentrepresentation import FormArgumentRepresentation 
 48   
 49  # FIXME: Need some refactoring of this class, 
 50  # for example splitting out geometry and function code in separate classes, 
 51  # and splitting out naming from symbol creation. 
52 -class FormRepresentationBase(object):
53 - def __init__(self, formdata, element_reps, geomrep, options):
54 55 # UFL data about the form 56 self.formdata = formdata 57 58 # Mapping : ufl.FiniteElement -> sfc.ElementRepresentation 59 self.element_reps = element_reps 60 61 # Shared geometry representation data 62 self.geomrep = geomrep 63 64 # SFC options 65 self.options = options 66 67 # Fetch some common variables from formdata for shorter names in rest of code 68 self.form = self.formdata.preprocessed_form 69 self.rank = self.formdata.rank 70 self.num_coefficients = self.formdata.num_coefficients 71 72 self.polygon_name = self.formdata.cell.domain() 73 self.domain = self.polygon_name # TODO: Rename everywhere and remove 74 75 self.signature = self.formdata.signature 76 77 # Some strings 78 self._define_classnames() 79 80 # Build running count of all _unique_ ufl elements. 81 # Not the same as running count of form arguments! 82 # mapping : ufl_element -> running count of unique elements 83 self.element_number = build_count_dict(self.formdata.sub_elements)
84
85 - def _define_classnames(self):
86 # Generate names for all classes 87 f = self.form 88 fd = self.formdata 89 90 self.classname = form_classname(fd.signature, self.options) 91 self.fe_names = [finite_element_classname(e) for e in fd.elements] 92 self.dm_names = [dof_map_classname(e) for e in fd.elements] 93 94 self.itg_names = dict((itg, integral_classname(itg, self.classname)) \ 95 for itg in self.form.integrals()) 96 97 self.citg_names = dict((itg.measure().domain_id(), self.itg_names[itg]) 98 for itg in f.cell_integrals()) 99 self.eitg_names = dict((itg.measure().domain_id(), self.itg_names[itg]) 100 for itg in f.exterior_facet_integrals()) 101 self.iitg_names = dict((itg.measure().domain_id(), self.itg_names[itg]) 102 for itg in f.interior_facet_integrals())
103
104 -class QuadratureRuleRepresentation(object):
105 - def __init__(self, form, sfc_cell, options):
106 # Determine quadrature rule from options or elements 107 self.quad_order = decide_quad_order(options, form) 108 109 # Pick a rule for cell 110 self.quad_rule = find_quadrature_rule(sfc_cell.shape, self.quad_order) 111 112 # Pick a rule for reference facet 113 reference_quad_rule = find_quadrature_rule(sfc_cell.facet_shape, 114 self.quad_order, family="gauss") 115 self.facet_quad_rules = map_facet_quadrature_rules(reference_quad_rule, sfc_cell)
116
117 -def decide_quad_order(options, form):
118 # Get integration order from options or elements 119 d = options.code.integral.integration_order 120 121 if d is None: 122 d = extract_max_quadrature_element_degree(form) 123 if d is None: 124 # TODO: Option for which estimation method to choose! 125 d = estimate_quadrature_degree(form) 126 if d is None: 127 d = estimate_total_polynomial_degree(form) 128 129 sfc_assert(d is not None, "All quadrature degree estimation rules failed!") 130 return d
131
132 -def map_facet_quadrature_rules(reference_quad_rule, sfc_cell):
133 # Map rule to each facet polygon on the reference cell 134 facet_quad_rules = [] 135 for facet, facet_polygon in enumerate(sfc_cell.facet_polygons): 136 137 # TODO: That dimensions get right here is rather magical, 138 # the global dimension set by initSyFi plays a major role here... 139 140 # affine mapping between reference 141 facet_G, facet_x0 = geometry_mapping(facet_polygon) 142 143 # scaling factor for diagonal facets, 144 # for scaling quadrature rule to a 145 # different reference domain, this is NOT the 146 # same as the global/local scaling factor facet_D! 147 D = 1.0 148 if facet == 0: 149 if sfc_cell.shape == "triangle": 150 D = math.sqrt(2) 151 elif sfc_cell.shape == "tetrahedron": 152 D = math.sqrt(3) 153 154 # scale weights by scaling factor 155 weights = [float(w)*D for w in reference_quad_rule.weights] 156 157 # apply affine mapping to quadrature points 158 points = [] 159 for p in reference_quad_rule.points: 160 # extend (n-1)D point p into nD point x: 161 x = swiginac.matrix(len(facet_x0), 1, list(p)+[0.0]) 162 # apply mapping in nD space: 163 newx = facet_G * x + facet_x0 164 # store as a list: 165 xlist = [float(newx[i].evalf()) for i in range(len(newx))] 166 points.append(xlist) 167 168 # make a copy of reference quadrature rule and modify it: 169 facet_quad_rule = copy.copy(reference_quad_rule) 170 171 if sfc_cell.shape != "interval": 172 facet_quad_rule.nsd += 1 173 assert facet_quad_rule.nsd == sfc_cell.nsd 174 facet_quad_rule.comment += "\nMapped to polygon in higher dimension.\n" 175 facet_quad_rule.weights = weights 176 facet_quad_rule.points = points 177 facet_quad_rules.append(facet_quad_rule) 178 return facet_quad_rules
179 180 # FIXME: Access subclasses through members instead, this is a partial refactoring state.
181 -class FormRepresentation(FormRepresentationBase, 182 GeometryRepresentation, 183 QuadratureRuleRepresentation, 184 FormArgumentRepresentation):
185 - def __init__(self, formdata, element_reps, geomrep, options):
186 if 0: 187 quadrep = QuadratureRuleRepresentation(formdata.preprocessed_form, 188 geomrep.sfc_cell, options) 189 farep = FormArgumentRepresentation(self) 190 191 FormRepresentationBase.__init__(self, formdata, element_reps, geomrep, options) 192 193 if 1: 194 GeometryRepresentation.__init__(self, formdata.cell.domain()) 195 QuadratureRuleRepresentation.__init__(self, formdata.preprocessed_form, 196 geomrep.sfc_cell, options) 197 FormArgumentRepresentation.__init__(self, self)
198
199 - def __str__(self):
200 s = "" 201 s += "FormRepresentation:\n" 202 s += "\n" # TODO: Add more here 203 s += GeometryRepresentation.__str__(self) 204 s += FormArgumentRepresentation.__str__(self) 205 return s
206 207 #def element_representation(self, iarg): 208 # return self.element_reps[self.formdata.elements[iarg]] 209 210 #def unique_element_number(self, iarg): 211 # return self.element_number[self.formdata.elements[iarg]] 212 213 #=============================================================================== 214 # def _build_argument_cache(self): 215 # use_symbols = self.integration_method == "quadrature" 216 # 217 # self.v_cache = [] 218 # self.Dv_cache = [] 219 # for iarg in range(self.rank + self.num_coefficients): 220 # elm = self.formdata.elements[iarg] 221 # rep = self.element_reps[elm] 222 # self.v_cache.append([]) 223 # self.Dv_cache.append([]) 224 # for i in range(rep.local_dimension): 225 # for component in rep.value_components: 226 # ve = self.v_expr(iarg, i, component) 227 # vs = self.v_sym(iarg, i, component) 228 # self.v_cache[iarg][i][component] = (vs, ve) 229 # 230 # for derivatives in [(d,) for d in range(self.sfc_cell.nsd)]: 231 # Dve = self.Dv_expr(iarg, i, component, derivatives, use_symbols) 232 # Dvs = self.Dv_sym(iarg, i, component, derivatives) 233 # self.Dv_cache[iarg][i][component][derivatives] = (Dvs, Dve) 234 #=============================================================================== 235