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, product 
 36   
 37  from sfc.symbolic_utils import grad, ddx, symbols, symbolic_matrix 
 38  from sfc.geometry import SFCCell 
 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  from sfc.representation.geometryrepresentation import create_syfi_polygon 
 43   
44 -def create_syfi_element(ufl_element, syfi_polygon, default_order=1):
45 "Create a basic element with SyFi." 46 sfc_debug("Entering create_syfi_element") 47 48 sfc_assert(not isinstance(ufl_element, ufl.MixedElement), 49 "Only creating SyFi elements for basic elements.") 50 f = ufl_element.family() 51 52 # Ensure that the element always have some degree, 53 # this should probably be handled somewhere else 54 d = ufl_element.degree() 55 if not isinstance(d, int): 56 d = default_order 57 58 if f in ("Lagrange", "CG"): 59 fe = SyFi.Lagrange(syfi_polygon, d) 60 elif f in ("Discontinuous Lagrange", "DG"): 61 if d == 0: 62 fe = SyFi.P0(syfi_polygon, 0) 63 else: 64 fe = SyFi.DiscontinuousLagrange(syfi_polygon, d) 65 elif f in ("Crouzeix-Raviart", "CR"): 66 fe = SyFi.CrouzeixRaviart(syfi_polygon, d) 67 elif f in ("Bubble", "B"): 68 fe = SyFi.Bubble(syfi_polygon, d) 69 70 elif f in ("Brezzi-Douglas-Marini", "BDM"): 71 raise NotImplementedError("Not implemented element family '%s'." % f) 72 elif f in ("Brezzi-Douglas-Fortin-Marini", "BDFM"): 73 raise NotImplementedError("Not implemented element family '%s'." % f) 74 elif f in ("Raviart-Thomas", "RT"): 75 raise NotImplementedError("Not implemented element family '%s'." % f) 76 elif f in ("Nedelec 1st kind H(div)", "N1div"): 77 raise NotImplementedError("Not implemented element family '%s'." % f) 78 elif f in ("Nedelec 2nd kind H(div)", "N2div"): 79 raise NotImplementedError("Not implemented element family '%s'." % f) 80 elif f in ("Nedelec 1st kind H(curl)", "N1curl"): 81 raise NotImplementedError("Not implemented element family '%s'." % f) 82 elif f in ("Nedelec 2nd kind H(curl)", "N2curl"): 83 raise NotImplementedError("Not implemented element family '%s'." % f) 84 elif f in ("Quadrature", "Q"): 85 raise NotImplementedError("Not implemented element family '%s'." % f) 86 elif f in ("Boundary Quadrature", "BQ"): 87 raise NotImplementedError("Not implemented element family '%s'." % f) 88 89 else: 90 raise NotImplementedError("Unknown element family '%s'." % f) 91 92 sfc_debug("Leaving create_syfi_element") 93 return fe
94 95 96 97 #=============================================================================== 98 # # swiginac.matrix(nsd, 1, dof_xi(fe,i)) 99 # def dof_xi(fe, i): 100 # """Return a swiginac column vector with the reference coordinates of dof i in fe, 101 # assuming elements with point evaluation dofs.""" 102 # dofi = fe.dof(i) 103 # # check if the element is a scalar or vector element 104 # if isinstance(dofi[0], swiginac.numeric): 105 # dofi0 = dofi 106 # elif isinstance(dofi[0], list): 107 # dofi0 = dofi[0] 108 # return dofi0 109 #=============================================================================== 110 111
112 -class ElementRepresentation(object):
113 __slots__ = (#Administrative data: 114 "options", "signature", 115 "dof_map_classname", "finite_element_classname", 116 # Element in other representations: 117 "ufl_element", "syfi_element", 118 # Cell data: 119 "geomrep", 120 "quad_rule", "ufl_cell", "polygon", "cell", 121 # Dimensions: 122 "local_dimension", "geometric_dimension", "topological_dimension", 123 # Value shape info: 124 "value_shape", "value_rank", "value_size", "value_components", 125 # Subelement data 126 "sub_elements", "sub_element_dof_offsets", "sub_element_value_offsets", 127 # Caches for basis functions 128 "_basis_function_cache", "_basis_function_derivative_cache", 129 # Coordinate information: 130 "dof_xi", "dof_x", 131 # Topology information: 132 "entity_dofs", "dof_entity", "num_entity_dofs", "facet_dofs", "num_facet_dofs", 133 # Geometry symbols: 134 "p0", "p", "G", "GinvT", 135 ) 136
137 - def __init__(self, ufl_element, geomrep, quad_rule=None, options=None, cache=None):
138 sfc_debug("Entering ElementRepresentation.__init__") 139 140 # Handle input and default values 141 assert isinstance(ufl_element, ufl.FiniteElementBase) 142 self.ufl_element = ufl_element 143 self.quad_rule = quad_rule 144 if options is None: 145 self.options = default_parameters() 146 else: 147 self.options = options 148 if cache is None: 149 cache = {} 150 151 # Some derived strings 152 self.signature = repr(self.ufl_element) 153 self.dof_map_classname = dof_map_classname(self.ufl_element) 154 self.finite_element_classname = finite_element_classname(self.ufl_element) 155 156 # Geometry information 157 # FIXME: Get stuff from geomrep 158 self.geomrep = geomrep 159 #self.ufl_cell = self.ufl_element.cell() 160 self.polygon = geomrep.syfi_polygon #create_syfi_polygon(self.ufl_cell.domain()) 161 self.cell = geomrep.sfc_cell #SFCCell(self.polygon) 162 163 self.geometric_dimension = self.cell.nsd 164 self.topological_dimension = self.cell.nsd 165 166 # Handy information about value shape 167 self.value_shape = self.ufl_element.value_shape() 168 self.value_rank = len(self.value_shape) 169 self.value_size = product(self.value_shape) # FIXME: Symmetries? 170 self.value_components = ufl.permutation.compute_indices(self.value_shape) 171 172 # Representations of subelements 173 self.sub_elements = [] 174 self.sub_element_dof_offsets = [] 175 self.sub_element_value_offsets = [] 176 177 # FIXME: It doesn't look like symmetries are handled properly? Write symmetry unit tests. 178 if isinstance(self.ufl_element, ufl.MixedElement): 179 dof_offset = 0 180 value_size_offset = 0 181 182 # Create ElementRepresentation objects for subelements, reuse if possible 183 for se in ufl_element.sub_elements(): 184 rep = cache.get(se, None) 185 if rep is None: 186 rep = ElementRepresentation(se, geomrep, self.quad_rule, self.options, cache) 187 cache[se] = rep 188 self.sub_elements.append(rep) 189 190 # Determine numbering offsets of subelements for dofs and values 191 self.sub_element_dof_offsets.append(dof_offset) 192 self.sub_element_value_offsets.append(value_size_offset) 193 194 dof_offset += rep.local_dimension 195 value_size_offset += rep.value_size 196 197 # Appending final sizes makes some algorithms more elegant 198 self.sub_element_dof_offsets.append(dof_offset) 199 self.sub_element_value_offsets.append(value_size_offset) 200 201 # No SyFi element, local dimension is the sum of subelement dimensions 202 self.syfi_element = None 203 self.local_dimension = dof_offset 204 205 elif self.ufl_element.family() == "Quadrature": 206 # No SyFi element, local dimension is the number of quadrature points 207 self.syfi_element = None 208 self.local_dimension = self.quad_rule.num_points 209 210 elif self.ufl_element.family() == "Boundary Quadrature": 211 # No SyFi element, local dimension is the number of quadrature points 212 # (TODO: On one facet or on all facets?) 213 sfc_error("Boundary Quadrature elements not implemented!") 214 self.syfi_element = None 215 self.local_dimension = self.facet_quad_rule.num_points # TODO: *num_facets? 216 217 elif self.ufl_element.family() == "Real": 218 # No SyFi element, local dimension equals value size 219 self.syfi_element = None 220 self.local_dimension = self.value_size 221 222 else: 223 # Make SyFi element 224 do = self.options.code.finite_element.default_order_of_element 225 self.syfi_element = create_syfi_element(self.ufl_element, self.polygon, default_order=do) 226 self.local_dimension = self.syfi_element.nbf() 227 228 # utility symbols 229 self._def_symbols() 230 231 # compute dof coordinates 232 self._precomp_coords() 233 234 # compute dof entity relations 235 self._build_entity_dofs() 236 self._count_entity_dofs() 237 self._build_dof_entity() 238 self._build_facet_dofs() 239 240 # initialize cache structures 241 self._basis_function_cache = {} 242 self._basis_function_derivative_cache = {} 243 244 sfc_debug("Leaving ElementRepresentation.__init__")
245
246 - def _def_symbols(self):
247 nsd = self.cell.nsd 248 249 # ... x,y,z symbols for convenience 250 #self.p0 = self.geomrep.x0_sym # FIXME: Use these 251 self.p = self.geomrep.xi_sym 252 self.p0 = swiginac.matrix(nsd, 1, symbols(["x0", "y0", "z0"][:nsd])) 253 #self.p = swiginac.matrix(nsd, 1, symbols(["x", "y", "z"][:nsd])) 254 255 # ... affine mapping symbols for convenience # TODO: Use J instead? 256 self.G = self.geomrep.G_sym 257 self.GinvT = symbolic_matrix(nsd, nsd, "GinvT")
258 #self.GinvT = self.geomrep.GinvT_sym # FIXME: Clean up transpose issues with Ginv 259
260 - def _dof_xi(self, i):
261 "Compute local dof coordinate of dof i." 262 nsd = self.cell.nsd 263 264 if self.sub_elements: 265 sub_element_index = self.dof_to_sub_element_index(i) 266 sub_dof = i - self.sub_element_dof_offsets[sub_element_index] 267 xi = self.sub_elements[sub_element_index]._dof_xi(sub_dof) 268 269 elif self.ufl_element.family() == "Quadrature": 270 xi = swiginac.matrix(nsd, 1, self.quad_rule.points[i]) 271 272 elif self.ufl_element.family() == "Real": 273 xi = swiginac.matrix(nsd, 1, [0.0]*nsd) 274 275 else: 276 # check if the element is a scalar or vector element 277 dofi = self.syfi_element.dof(i) 278 if isinstance(dofi[0], swiginac.numeric): 279 # scalar element 280 dof_xi_list = dofi 281 elif isinstance(dofi[0], list): 282 # vector element 283 if isinstance(dofi[0][0], list): 284 # compute midpoints 285 midpoint = [0 for i in range(len(dofi[0][0]))] 286 for d in dofi[0][0:]: 287 for p, dp in enumerate(d): 288 midpoint[p] += dp 289 for p in range(len(d)): 290 midpoint[p] /= len(dofi[0]) 291 dof_xi_list = midpoint 292 else: 293 # use coordinate directly 294 dof_xi_list = dofi[0] 295 xi = swiginac.matrix(nsd, 1, dof_xi_list) 296 297 return xi
298
299 - def _precomp_coords(self):
300 "Precompute dof coordinates." 301 self.dof_xi = [] 302 self.dof_x = [] 303 304 for i in range(self.local_dimension): 305 # point coordinates for this dof in reference coordinates 306 dof_xi = self._dof_xi(i) 307 self.dof_xi.append(dof_xi) 308 309 # apply geometry mapping to get global coordinates 310 dof_x = (self.G.mul(dof_xi).add(self.p0)).evalm() # TODO: Assumes affine mapping! 311 self.dof_x.append(dof_x)
312
313 - def _build_entity_dofs(self): # XXXReal: Does not make sense for Real
314 "Build dof vs mesh entity relations." 315 316 # TODO: This may be optimized if necessary for mixed elements, but maybe we don't care. 317 318 # Build empty basic structure with a list of lists of lists 319 # entity_dofs[top_dim][entity_num] = [] 320 entity_dofs = [ [ [] for j in range(self.cell.num_entities[i]) ] 321 for i in range(self.topological_dimension+1) ] 322 323 fam = self.ufl_element.family() 324 if fam in ("Discontinuous Lagrange", "Quadrature"): 325 # associate all dofs with the cell 326 entity_dofs[self.topological_dimension][0].extend(range(self.local_dimension)) 327 328 elif fam == "Boundary Quadrature": 329 sfc_error("Boundary Quadrature not handled.") 330 331 elif fam == "Real": 332 pass #sfc_error("Real elements not handled.") 333 334 else: 335 # NB! Assuming location dof_xi coordinates match topological entity! 336 # build dof list for each cell entity 337 for k in range(self.local_dimension): 338 (d, i) = self.cell.find_entity(self.dof_xi[k]) 339 entity_dofs[d][i].append(k) 340 341 self.entity_dofs = entity_dofs
342
343 - def _count_entity_dofs(self):
344 # count number of dofs per entity 345 self.num_entity_dofs = tuple(len(self.entity_dofs[d][0]) 346 for d in range(self.topological_dimension+1)) 347 348 # assert that all entities have the same number of associated dofs 349 # (there's a theoretical risk of floating point comparisons messing up find_entity) 350 for d in range(self.topological_dimension+1): 351 # each doflist is a list of local dofs associated 352 # with a particular mesh entity of dimension d 353 sfc_assert(all(len(doflist) == self.num_entity_dofs[d] 354 for doflist in self.entity_dofs[d]), 355 "Entities have different number of associated dofs.")
356
357 - def _build_dof_entity(self):
358 # Build the inverse mapping: idof -> ((d, i), j) (Not currently used for anything) 359 dof_entity = [None]*self.local_dimension 360 for d in range(self.topological_dimension+1): 361 for i in range(self.cell.num_entities[d]): 362 for j, k in enumerate(self.entity_dofs[d][i]): 363 sfc_assert(dof_entity[k] is None, "Expected to set each dof only once.") 364 dof_entity[k] = (d, i, j) 365 self.dof_entity = dof_entity
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 # cache hit? 415 key = (i, component) 416 N = self._basis_function_cache.get(key) 417 if N is not None: 418 return N 419 420 if self.sub_elements: 421 # get sub element representation corresponding to this dof 422 sub_element_index = self.dof_to_sub_element_index(i) 423 sub_element = self.sub_elements[sub_element_index] 424 425 # dof in sub element numbering 426 sub_dof = i - self.sub_element_dof_offsets[sub_element_index] 427 428 # component in flattened sub element value index space 429 comp_index = component_to_index(component, self.value_shape) 430 value_offset = self.sub_element_value_offsets[sub_element_index] 431 432 # check that the component is in the value range of the subelement 433 is_nonzero = ( (comp_index >= value_offset) and 434 (comp_index < (value_offset + sub_element.value_size)) ) 435 if is_nonzero: 436 # component in unflattened sub element value index space 437 sub_comp_index = comp_index - value_offset 438 sub_component = index_to_component(sub_comp_index, sub_element.value_shape) 439 440 # basis_function from computed component of subelement 441 N = sub_element.basis_function(sub_dof, sub_component) 442 443 else: 444 N = swiginac.numeric(0.0) 445 446 elif "Quadrature" in self.ufl_element.family(): 447 sfc_error("Cannot compute basis functions for quadrature element.") 448 N = swiginac.numeric(1.0) 449 450 elif "Real" in self.ufl_element.family(): 451 # XXXReal: Will we use this or not? 452 #sfc_error("Cannot compute basis functions for Real element.") 453 N = swiginac.numeric(1.0) 454 455 else: 456 # Basic element, get basis function from SyFi 457 N = self.syfi_element.N(i) 458 if isinstance(N, swiginac.matrix): 459 sfc_assert((int(N.rows()), int(N.cols())) == self.value_shape, "Shape mismatch") 460 return N[component] 461 sfc_assert(component == (), 462 ("Found scalar basic element, expecting no component, got %s." 463 % repr(component))) 464 465 # store in cache 466 self._basis_function_cache[key] = N 467 return N
468
469 - def basis_function_derivative(self, i, component, directions):
470 # d/dx and d/dy commute so we sort the derivative variables: 471 directions = tuple(sorted(directions)) 472 key = (i, component, directions) 473 474 # cache hit? 475 DN = self._basis_function_derivative_cache.get(key, None) 476 if DN is None: 477 # compute derivative 478 DN = self.basis_function(i, component) 479 for j in directions: 480 DN = ddx(DN, j, self.GinvT) 481 # store in cache 482 self._basis_function_derivative_cache[key] = DN 483 484 return DN
485 486 487 # def component_to_sub_element_index(self, component): 488 # sfc_assert(self.sub_elements, "Only mixed elements have sub elements.") 489 # 490 # # FIXME! 491 # sfc_error("FIXME: In component_to_sub_element_index: not sure where this will be used?") 492 # 493 # component_value = flatten_component(component, self.value_shape) 494 # 495 # n = len(self.sub_elements) 496 # for i in range(1,n+1): 497 # if component_value < self.sub_element_value_offsets[i]: 498 # sub_component_value = self.sub_element_value_offsets[i] - ccomponent_value 499 # sub_element_index = i-1 500 # sub_element = self.sub_elements[sub_element_index] 501 # sub_component = unflatten_component(sub_component_value, sub_element) 502 # return sub_component 503 # sfc_error("Invalid component value!") 504 # 505 # #sub_element_component, sub_element =\ 506 # # self.ufl_element.extract_component(component) 507 # #return sub_element_component, sub_element 508