SyFi  0.3
cell_assembly.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 import ufl
00003 import dolfin
00004 from dolfin import Mesh, MeshEditor, assemble
00005 
00006 cell2dim = { "interval": 1,
00007              "triangle": 2,
00008              "tetrahedron": 3,
00009              "quadrilateral": 2,
00010              "hexahedron": 3,
00011              ufl.interval: 1,
00012              ufl.triangle: 2,
00013              ufl.tetrahedron: 3,
00014              ufl.quadrilateral: 2,
00015              ufl.hexahedron: 3,
00016              }
00017 
00018 cell2volume = { "interval": 1.0,
00019                 "triangle": 0.5,
00020                 "tetrahedron": 1.0/6.0,
00021                 "quadrilateral": 1.0,
00022                 "hexahedron": 1.0,
00023                 ufl.interval: 1.0,
00024                 ufl.triangle: 0.5,
00025                 ufl.tetrahedron: 1.0/6.0,
00026                 ufl.quadrilateral: 1.0,
00027                 ufl.hexahedron: 1.0,
00028                 }
00029 
00030 def UnitCell(celltype):
00031     tdim = cell2dim[celltype]
00032     gdim = tdim
00033     mesh = Mesh()
00034     editor = MeshEditor()
00035     editor.open(mesh, celltype, tdim, gdim)
00036     if celltype == "interval":
00037         vertices = [(0.0,),
00038                     (1.0,)]
00039     if celltype == "triangle":
00040         vertices = [(0.0, 0.0),
00041                     (1.0, 0.0),
00042                     (0.0, 1.0)]
00043     if celltype == "tetrahedron":
00044         vertices = [(0.0, 0.0, 0.0),
00045                     (1.0, 0.0, 0.0),
00046                     (0.0, 1.0, 0.0),
00047                     (0.0, 0.0, 1.0)]
00048     if celltype == "quadrilateral":
00049         vertices = [(0.0, 0.0),
00050                     (1.0, 0.0),
00051                     (1.0, 1.0),
00052                     (0.0, 1.0)]
00053     if celltype == "hexahedron":
00054         vertices = [(0.0, 0.0, 0.0),
00055                     (1.0, 0.0, 0.0),
00056                     (1.0, 1.0, 0.0),
00057                     (0.0, 1.0, 0.0),
00058                     (0.0, 0.0, 1.0),
00059                     (1.0, 0.0, 1.0),
00060                     (1.0, 1.0, 1.0),
00061                     (0.0, 1.0, 1.0)]
00062     editor.init_vertices(len(vertices))
00063     editor.init_cells(1)
00064     for i, p in enumerate(vertices):
00065         editor.add_vertex(i, *p)
00066     editor.add_cell(0, *range(len(vertices)))
00067     editor.close()
00068     return mesh
00069 
00070 def element_to_function_space(element, mesh):
00071     if isinstance(element, ufl.FiniteElement):
00072         return dolfin.FunctionSpace(mesh, element.family(), element.degree())
00073     if isinstance(element, ufl.VectorElement):
00074         dim = len(element._sub_elements) # Hack!
00075         return dolfin.VectorFunctionSpace(mesh, element.family(), element.degree(), dim)
00076     if isinstance(element, ufl.TensorElement):
00077         shape = element._shape # Hack!
00078         return dolfin.TensorFunctionSpace(mesh, element.family(), element.degree(),
00079                                           shape, element.symmetry())
00080     raise TypeError("Invalid element type, cannot construct function space.")
00081 
00082 def as_coeff(c):
00083     if isinstance(c, dolfin.cpp.GenericFunction):
00084         return c
00085     elif (isinstance(c, (float,int)) or
00086           (isinstance(c, tuple) and isinstance(c[0],(float,int))) or
00087           (isinstance(c, tuple) and isinstance(c[0],tuple) and isinstance(c[0][0],(float,int)))):
00088         return dolfin.Constant(c)
00089     elif isinstance(c, str) or (isinstance(c, tuple) and isinstance(c[0],str)):
00090         return dolfin.Expression(c)
00091     else:
00092         raise TypeError("Invalid coefficient type %s." % str(type(c)))
00093 
00094 def assemble_on_cell(form, celltype, coeffs, elements=()):
00095     "Assemble UFC form on a unit cell mesh and return the result as a float or numpy array."
00096 
00097     if not isinstance(celltype, str):
00098         celltype = celltype.domain()
00099 
00100     # Convert data to dolfin formats
00101     mesh = UnitCell(celltype)
00102     coefficients = [as_coeff(c) for c in coeffs]
00103     function_spaces = [element_to_function_space(e, mesh) for e in elements]
00104 
00105     # Assemble ufc form with dolfin
00106     A = assemble(form, mesh=mesh, function_spaces=function_spaces, coefficients=coefficients)
00107 
00108     # Return in python/numpy format
00109     if isinstance(A, float):
00110         return A
00111     return A.array()
00112 
00113 def num_integrals(form):
00114     return (form.num_cell_domains(),
00115             form.num_exterior_facet_domains(),
00116             form.num_interior_facet_domains())
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines