SyFi
0.3
|
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())