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 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.10354 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)8486 # 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())116106 # 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)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 d131133 # 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_rules179 180 # FIXME: Access subclasses through members instead, this is a partial refactoring state.181 -class FormRepresentation(FormRepresentationBase, 182 GeometryRepresentation, 183 QuadratureRuleRepresentation, 184 FormArgumentRepresentation):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 #=============================================================================== 235186 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)198200 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
Home | Trees | Indices | Help |
|
---|
Generated by Epydoc 3.0.1 on Thu Jan 31 03:52:26 2013 | http://epydoc.sourceforge.net |