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

Source Code for Module sfc.representation.elementrepresentation

  1  #!/usr/bin/env python 
  2  # -*- coding: utf-8 -*- 
  3  """ 
  4  This module contains functions and classes converting  
  5  from UFL to internal SyFi representations of elements. 
  6  """ 
  7   
  8  # Copyright (C) 2008 Martin Sandve Alnes and Simula Resarch Laboratory 
  9  # 
 10  # This file is part of SyFi. 
 11  # 
 12  # SyFi is free software: you can redistribute it and/or modify 
 13  # it under the terms of the GNU General Public License as published by 
 14  # the Free Software Foundation, either version 2 of the License, or 
 15  # (at your option) any later version. 
 16  # 
 17  # SyFi is distributed in the hope that it will be useful, 
 18  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 19  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 20  # GNU General Public License for more details. 
 21  # 
 22  # You should have received a copy of the GNU General Public License 
 23  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 24  # 
 25  # Modified by Kent-Andre Mardal, 2010. 
 26  # 
 27  # First added:  2008-08-13 
 28  # Last changed: 2008-12-16 
 29   
 30  import operator 
 31  import ufl 
 32  import swiginac 
 33  import SyFi 
 34   
 35  from ufl.common import component_to_index, index_to_component 
 36   
 37  from sfc.symbolic_utils import grad, ddx, symbols, symbolic_matrix 
 38  from sfc.geometry import UFCCell 
 39  from sfc.common.names import finite_element_classname, dof_map_classname 
 40  from sfc.common import sfc_assert, sfc_info, sfc_warning, sfc_error, sfc_debug 
 41  from sfc.common.options import default_parameters 
 42   
43 -def product(sequence):
44 return reduce(operator.__mul__, sequence, 1)
45 46 # TODO: Replace nsd with topological and geometric dimensions 47 # everywhere for clarity and future flexibility 48 49
50 -def test_polygon(polygon):
51 print type(polygon) 52 print dir(polygon) 53 print polygon.no_space_dim() 54 print polygon.str()
55
56 -def create_syfi_polygon(cell):
57 sfc_debug("Entering create_syfi_polygon") 58 if cell == "interval": p = SyFi.ReferenceLine() 59 elif cell == "triangle": p = SyFi.ReferenceTriangle() 60 elif cell == "tetrahedron": p = SyFi.ReferenceTetrahedron() 61 elif cell == "quadrilateral": p = SyFi.ReferenceRectangle() 62 elif cell == "hexahedron": p = SyFi.ReferenceBox() 63 else: raise ValueError("Unknown element cell '%s'." % cell) 64 sfc_debug("Leaving create_syfi_polygon") 65 return p
66
67 -def create_syfi_element(e, polygon, default_order=1):
68 "Create a basic element with SyFi." 69 sfc_debug("Entering create_syfi_element") 70 71 sfc_assert(not isinstance(e, ufl.MixedElement), "Only creating SyFi elements for basic elements.") 72 f = e.family() 73 74 # ensure that the element always have some degree, this should probably be 75 d = e.degree() 76 if not isinstance(d, int): d = default_order 77 78 if f in ("Lagrange", "CG"): 79 fe = SyFi.Lagrange(polygon, d) 80 elif f in ("Discontinuous Lagrange", "DG"): 81 if d == 0: fe = SyFi.P0(polygon, 0) 82 else: fe = SyFi.DiscontinuousLagrange(polygon, d) 83 elif f in ("Crouzeix-Raviart", "CR"): 84 fe = SyFi.CrouzeixRaviart(polygon, d) 85 elif f in ("Bubble", "B"): 86 fe = SyFi.Bubble(polygon, d) 87 elif f in ("Brezzi-Douglas-Marini", "BDM"): 88 raise NotImplementedError("Not implemented element family '%s'." % f) 89 elif f in ("Brezzi-Douglas-Fortin-Marini", "BDFM"): 90 raise NotImplementedError("Not implemented element family '%s'." % f) 91 elif f in ("Raviart-Thomas", "RT"): 92 raise NotImplementedError("Not implemented element family '%s'." % f) 93 elif f in ("Nedelec 1st kind H(div)", "N1div"): 94 raise NotImplementedError("Not implemented element family '%s'." % f) 95 elif f in ("Nedelec 2nd kind H(div)", "N2div"): 96 raise NotImplementedError("Not implemented element family '%s'." % f) 97 elif f in ("Nedelec 1st kind H(curl)", "N1curl"): 98 raise NotImplementedError("Not implemented element family '%s'." % f) 99 elif f in ("Nedelec 2nd kind H(curl)", "N2curl"): 100 raise NotImplementedError("Not implemented element family '%s'." % f) 101 elif f in ("Quadrature", "Q"): 102 raise NotImplementedError("Not implemented element family '%s'." % f) 103 elif f in ("Boundary Quadrature", "BQ"): 104 raise NotImplementedError("Not implemented element family '%s'." % f) 105 else: 106 raise NotImplementedError("Unknown element family '%s'." % f) 107 108 sfc_debug("Leaving create_syfi_element") 109 return fe
110 111 112 113 #=============================================================================== 114 # # swiginac.matrix(nsd, 1, dof_xi(fe,i)) 115 # def dof_xi(fe, i): 116 # """Return a swiginac column vector with the reference coordinates of dof i in fe, 117 # assuming elements with point evaluation dofs.""" 118 # dofi = fe.dof(i) 119 # # check if the element is a scalar or vector element 120 # if isinstance(dofi[0], swiginac.numeric): 121 # dofi0 = dofi 122 # elif isinstance(dofi[0], list): 123 # dofi0 = dofi[0] 124 # return dofi0 125 #=============================================================================== 126 127
128 -class ElementRepresentation(object):
129 __slots__ = (#Administrative data: 130 "options", "signature", 131 "dof_map_classname", "finite_element_classname", 132 # Element in other representations: 133 "ufl_element", "syfi_element", 134 # Cell data: 135 "quad_rule", "ufl_cell", "polygon", "cell", 136 # Dimensions: 137 "local_dimension", "geometric_dimension", "topological_dimension", 138 # Value shape info: 139 "value_shape", "value_rank", "value_size", "value_components", 140 # Subelement data 141 "sub_elements", "sub_element_dof_offsets", "sub_element_value_offsets", 142 # Caches for basis functions 143 "_basis_function_cache", "_basis_function_derivative_cache", 144 # Coordinate information: 145 "dof_xi", "dof_x", 146 # Topology information: 147 "entity_dofs", "dof_entity", "num_entity_dofs", "facet_dofs", "num_facet_dofs", 148 # Geometry symbols: 149 "p0", "p", "G", "GinvT", 150 ) 151
152 - def __init__(self, ufl_element, quad_rule=None, options=None, cache=None):
153 sfc_debug("Entering ElementRepresentation.__init__") 154 155 # Handle input and default values 156 assert isinstance(ufl_element, ufl.FiniteElementBase) 157 self.ufl_element = ufl_element 158 self.quad_rule = quad_rule 159 if options is None: 160 self.options = default_parameters() 161 else: 162 self.options = options 163 if cache is None: 164 cache = {} 165 166 # Some derived strings 167 self.signature = repr(self.ufl_element) 168 self.dof_map_classname = dof_map_classname(self.ufl_element) 169 self.finite_element_classname = finite_element_classname(self.ufl_element) 170 171 # Geometry information 172 self.ufl_cell = self.ufl_element.cell() 173 self.polygon = create_syfi_polygon(self.ufl_cell.domain()) 174 self.cell = UFCCell(self.polygon) 175 176 self.geometric_dimension = self.cell.nsd 177 self.topological_dimension = self.cell.nsd 178 179 # Handy information about value shape 180 self.value_shape = self.ufl_element.value_shape() 181 self.value_rank = len(self.value_shape) 182 self.value_size = product(self.value_shape) 183 self.value_components = ufl.permutation.compute_indices(self.value_shape) 184 185 # Representations of subelements 186 self.sub_elements = [] 187 self.sub_element_dof_offsets = [] 188 self.sub_element_value_offsets = [] 189 190 if isinstance(self.ufl_element, ufl.MixedElement): 191 dof_offset = 0 192 value_size_offset = 0 193 194 # Create ElementRepresentation objects for subelements, reuse if possible 195 for se in ufl_element.sub_elements(): 196 rep = cache.get(se, None) 197 if rep is None: 198 rep = ElementRepresentation(se, self.quad_rule, self.options, cache) 199 cache[se] = rep 200 self.sub_elements.append(rep) 201 202 # Determine numbering offsets of subelements for dofs and values 203 self.sub_element_dof_offsets.append(dof_offset) 204 self.sub_element_value_offsets.append(value_size_offset) 205 206 dof_offset += rep.local_dimension 207 value_size_offset += rep.value_size 208 209 # Appending final sizes makes some algorithms more elegant 210 self.sub_element_dof_offsets.append(dof_offset) 211 self.sub_element_value_offsets.append(value_size_offset) 212 213 # No SyFi element, local dimension is the sum of subelement dimensions 214 self.syfi_element = None 215 self.local_dimension = dof_offset 216 217 elif self.ufl_element.family() == "Quadrature": 218 # No SyFi element, local dimension is the number of quadrature points 219 self.syfi_element = None 220 self.local_dimension = self.quad_rule.num_points 221 222 elif self.ufl_element.family() == "Boundary Quadrature": 223 # No SyFi element, local dimension is the number of quadrature points (TODO: On one facet or on all facets?) 224 sfc_error("Boundary Quadrature elements not implemented!") 225 self.syfi_element = None 226 self.local_dimension = self.facet_quad_rule.num_points # TODO: *num_facets? 227 228 elif self.ufl_element.family() == "Real": 229 # No SyFi element, local dimension equals value size 230 self.syfi_element = None 231 self.local_dimension = self.value_size 232 233 else: 234 # Make SyFi element 235 self.syfi_element = create_syfi_element(self.ufl_element, self.polygon, default_order=self.options.code.finite_element.default_order_of_element) 236 self.local_dimension = self.syfi_element.nbf() 237 238 # utility symbols 239 self._def_symbols() 240 241 # compute dof coordinates 242 self._precomp_coords() 243 244 # compute dof entity relations 245 self._build_entity_dofs() 246 self._build_facet_dofs() 247 248 # initialize cache structures 249 self._basis_function_cache = {} 250 self._basis_function_derivative_cache = {} 251 252 sfc_debug("Leaving ElementRepresentation.__init__")
253
254 - def _def_symbols(self):
255 nsd = self.cell.nsd 256 257 # ... x,y,z symbols for convenience 258 self.p0 = swiginac.matrix(nsd, 1, symbols(["x0", "y0", "z0"][:nsd])) 259 self.p = swiginac.matrix(nsd, 1, symbols(["x", "y", "z"][:nsd])) 260 #self.x = swiginac.matrix(nsd, 1, symbols(["x0", "x1", "x2"][:nsd])) # TODO: Use these everywhere! self.x are global coordinates, self.xi are locals. 261 #self.xi = swiginac.matrix(nsd, 1, symbols(["xi0", "xi1", "xi2"][:nsd])) 262 263 # ... affine mapping symbols for convenience # TODO: Use J instead? 264 self.G = symbolic_matrix(nsd, nsd, "G") 265 self.GinvT = symbolic_matrix(nsd, nsd, "GinvT")
266
267 - def _dof_xi(self, i):
268 "Compute local dof coordinate of dof i." 269 nsd = self.cell.nsd 270 271 if self.sub_elements: 272 sub_element_index = self.dof_to_sub_element_index(i) 273 sub_dof = i - self.sub_element_dof_offsets[sub_element_index] 274 xi = self.sub_elements[sub_element_index]._dof_xi(sub_dof) 275 276 elif self.ufl_element.family() == "Quadrature": 277 xi = swiginac.matrix(nsd, 1, self.quad_rule.points[i]) 278 279 elif self.ufl_element.family() == "Real": 280 xi = swiginac.matrix(nsd, 1, [0.0]*nsd) 281 282 else: 283 # check if the element is a scalar or vector element 284 dofi = self.syfi_element.dof(i) 285 if isinstance(dofi[0], swiginac.numeric): 286 # scalar element 287 dof_xi_list = dofi 288 elif isinstance(dofi[0], list): 289 # vector element 290 if isinstance(dofi[0][0], list): 291 # compute midpoints 292 midpoint = [0 for i in range(len(dofi[0][0]))] 293 for d in dofi[0][0:]: 294 for p, dp in enumerate(d): 295 midpoint[p] += dp 296 for p in range(len(d)): 297 midpoint[p] /= len(dofi[0]) 298 dof_xi_list = midpoint 299 else: 300 # use coordinate directly 301 dof_xi_list = dofi[0] 302 xi = swiginac.matrix(nsd, 1, dof_xi_list) 303 304 return xi
305
306 - def _precomp_coords(self):
307 "Precompute dof coordinates." 308 self.dof_xi = [] 309 self.dof_x = [] 310 311 for i in range(self.local_dimension): 312 # point coordinates for this dof in reference coordinates 313 dof_xi = self._dof_xi(i) 314 self.dof_xi.append(dof_xi) 315 316 # apply geometry mapping to get global coordinates 317 dof_x = (self.G.mul(dof_xi).add(self.p0)).evalm() # TODO: Assumes affine mapping! 318 self.dof_x.append(dof_x)
319
320 - def _build_entity_dofs(self): # XXXReal: Does not make sense for Real
321 "Build dof vs mesh entity relations." 322 323 # TODO: This may be optimized if necessary for mixed elements, but maybe we don't care. 324 325 # The basic structure we're building here is a list of lists of lists 326 self.entity_dofs = [] 327 for i in range(self.topological_dimension+1): 328 lists = [[] for j in range(self.cell.num_entities[i])] 329 self.entity_dofs.append(lists) 330 331 if self.ufl_element.family() in ("Discontinuous Lagrange", "Quadrature"): 332 # associate all dofs with the cell 333 self.entity_dofs[self.topological_dimension][0] = list(range(self.local_dimension)) 334 335 elif self.ufl_element.family() == "Boundary Quadrature": 336 sfc_error("Boundary Quadrature not handled.") 337 338 elif self.ufl_element.family() == "Real": 339 pass #sfc_error("Real elements not handled.") 340 341 else: 342 # NB! Assuming location dof_xi coordinates match topological entity! 343 # build dof list for each cell entity 344 for k in range(self.local_dimension): 345 (d, i) = self.cell.find_entity(self.dof_xi[k]) 346 self.entity_dofs[d][i].append(k) 347 348 # Build the inverse mapping: idof -> ((d, i), j) (Not currently used for anything) 349 self.dof_entity = [None]*self.local_dimension 350 for d in range(self.topological_dimension+1): 351 for i in range(self.cell.num_entities[d]): 352 for j, k in enumerate(self.entity_dofs[d][i]): 353 sfc_assert(self.dof_entity[k] is None, "Expected to set each dof only once.") 354 self.dof_entity[k] = (d, i, j) 355 356 # count number of dofs per entity 357 self.num_entity_dofs = tuple(len(self.entity_dofs[d][0]) for d in range(self.topological_dimension+1)) 358 359 # assert that all entities have the same number of associated dofs 360 # (there's a theoretical risk of floating point comparisons messing up the above logic) 361 for d in range(self.topological_dimension+1): 362 for doflist in self.entity_dofs[d]: 363 # each doflist is a list of local dofs associated 364 # with a particular mesh entity of dimension d 365 assert len(doflist) == self.num_entity_dofs[d]
366
367 - def _build_facet_dofs(self):
368 "Build facet vs dof relations." 369 # ... build a list of dofs for each facet: 370 self.facet_dofs = [[] for j in range(self.cell.num_facets)] 371 372 if self.ufl_element.family() in ("Discontinuous Lagrange", "Quadrature", "Real"): 373 pass # no dofs on facets 374 375 elif self.ufl_element.family() == "Boundary Quadrature": 376 sfc_error("Boundary Quadrature not handled.") 377 378 else: 379 # for each facet j, loop over the reference coordinates 380 # for all dofs i and check if the dof is on the facet 381 for j in range(self.cell.num_facets): 382 for (i,p) in enumerate(self.dof_xi): 383 if self.cell.facet_check(j, p): 384 self.facet_dofs[j].append(i) 385 386 # ... count number of dofs for each facet (assuming this is constant!) 387 self.num_facet_dofs = len(self.facet_dofs[0]) 388 389 # verify that this number is constant for all facets 390 sfc_assert(all(len(fdofs) == self.num_facet_dofs for fdofs in self.facet_dofs), 391 "Not the same number of dofs on each facet. This breaks an assumption in UFC.")
392 393 # --- subelement access 394
395 - def dof_to_sub_element_index(self, dof):
396 "Return the index of the sub element the given dof is part of." 397 n = len(self.sub_elements) 398 sfc_assert(n, "Only mixed elements have sub elements.") 399 for i in range(n+1): 400 if dof < self.sub_element_dof_offsets[i]: 401 return i-1 402 sfc_error("Invalid dof value!")
403
404 - def sub_element_to_dofs(self, i):
405 "Return a list of all dof indices for sub element with index i." 406 sfc_assert(self.sub_elements, "Only mixed elements have sub elements.") 407 a = self.sub_element_dof_offsets[i] 408 b = self.sub_element_dof_offsets[i+1] 409 return range(a, b)
410 411 # --- function space 412
413 - def basis_function(self, i, component):
414 # hit cache? 415 N = self._basis_function_cache.get((i,component), None) 416 if N is not None: 417 return N 418 419 if self.sub_elements: 420 # get sub element representation corresponding to this dof 421 sub_element_index = self.dof_to_sub_element_index(i) 422 sub_element = self.sub_elements[sub_element_index] 423 424 # dof in sub element numbering 425 sub_dof = i - self.sub_element_dof_offsets[sub_element_index] 426 427 # component in flattened sub element value index space 428 comp_index = component_to_index(component, self.value_shape) 429 value_offset = self.sub_element_value_offsets[sub_element_index] 430 431 # check that the component is in the value range of the subelement 432 is_nonzero = (comp_index >= value_offset) and (comp_index < (value_offset + sub_element.value_size)) 433 if is_nonzero: 434 # component in unflattened sub element value index space 435 sub_comp_index = comp_index - value_offset 436 sub_component = index_to_component(sub_comp_index, sub_element.value_shape) 437 438 # basis_function from computed component of subelement 439 N = sub_element.basis_function(sub_dof, sub_component) 440 441 else: 442 N = swiginac.numeric(0.0) 443 444 elif "Quadrature" in self.ufl_element.family(): 445 sfc_error("Cannot compute basis functions for quadrature element.") 446 N = swiginac.numeric(1.0) 447 448 elif "Real" in self.ufl_element.family(): 449 #sfc_error("Cannot compute basis functions for Real element.") # XXXReal: Will we use this or not? 450 N = swiginac.numeric(1.0) 451 452 else: 453 # Basic element, get basis function from SyFi 454 N = self.syfi_element.N(i) 455 if isinstance(N, swiginac.matrix): 456 sfc_assert((int(N.rows()), int(N.cols())) == self.value_shape, "Shape mismatch") 457 return N[component] 458 sfc_assert(component == (), "Found scalar basic element, expecting no component, got %s." % repr(component)) 459 460 # put in cache 461 self._basis_function_cache[(i,component)] = N 462 return N
463
464 - def basis_function_derivative(self, i, component, directions):
465 # d/dx and d/dy commute so we sort the derivative variables: 466 directions = tuple(sorted(directions)) 467 468 # cache hit? 469 DN = self._basis_function_derivative_cache.get((i,component,directions), None) 470 if DN is not None: 471 return DN 472 473 # compute derivative 474 DN = self.basis_function(i, component) 475 for j in directions: 476 DN = ddx(DN, j, self.GinvT) 477 478 # put in cache 479 self._basis_function_derivative_cache[(i,component,directions)] = DN 480 return DN
481 482 483 # def component_to_sub_element_index(self, component): 484 # sfc_assert(self.sub_elements, "Only mixed elements have sub elements.") 485 # 486 # # FIXME! 487 # sfc_error("FIXME: In component_to_sub_element_index: not sure where this will be used?") 488 # 489 # component_value = flatten_component(component, self.value_shape) 490 # 491 # n = len(self.sub_elements) 492 # for i in range(1,n+1): 493 # if component_value < self.sub_element_value_offsets[i]: 494 # sub_component_value = self.sub_element_value_offsets[i] - ccomponent_value 495 # sub_element_index = i-1 496 # sub_element = self.sub_elements[sub_element_index] 497 # sub_component = unflatten_component(sub_component_value, sub_element) 498 # return sub_component 499 # sfc_error("Invalid component value!") 500 # 501 # #sub_element_component, sub_element =\ 502 # # self.ufl_element.extract_component(component) 503 # #return sub_element_component, sub_element 504