Actually Paraview supports quadratic elements pretty well. You can select
the number of "nonlinear subdivisions" for each object (on the properties
panel) and improve the resolution.
I have a draft script that exports a P1/P2 solution to xdmf from dolfin.
Still work in progress. Let me know if you plan to improve it or merge it
into dolfin.
You'll find it as attachment (run it and then try to wrap the domain with
the "displacement" field).
Simone
2014-02-07 21:31 GMT+01:00 Jed Brown <[email protected]>:
> David Bernstein <[email protected]> writes:
> > If you have it lying around could you send me the xdmf file with
> > quadratic elements? If I have time I might see if someone at kitware
> > knows what the problem is.
>
> It's not a "problem" per se. There have been research papers on
> visualizing high-order basis functions, but all the production software
> decimates to linear triangles. As far as VTK is concerned,
> VTK_QUADRATIC_TRIANGLE is just a compact way to specify these smaller
> triangles; it doesn't really visualize the higher order.
>
> _______________________________________________
> fenics mailing list
> [email protected]
> http://fenicsproject.org/mailman/listinfo/fenics
>
>
from dolfin import *
import sys
import numpy as np
import cStringIO
from textwrap import dedent
def xdmf_dataitem_dump(data, fmt) :
if isinstance(data, str) :
# assuming HDF5 reference
return "HDF", data
elif isinstance(data, np.ndarray) :
# dumping data into str
data_str = cStringIO.StringIO()
np.savetxt(data_str, data, fmt = fmt)
return "XML", data_str.getvalue()
else :
raise ValueError
def xdmf_topology(connectivity, num_cell, degree, geodim) :
vtk_types = [ [ "Polyline", "Triangle", "Tetrahedron" ], [ "Edge_3", "Tri_6", "Tet_10" ] ]
cell_dofs = [ [ 1, 3, 4 ], [ 3, 6, 10 ] ]
dataitem_type, connectivity_data = xdmf_dataitem_dump(connectivity, '%d')
topology_data = '''\
<Topology NumberOfElements="%(num_cell)d" TopologyType="%(cell_type)s">
<DataItem Format="%(dataitem_type)s" DataType="Int" Dimensions="%(num_cell)d %(cell_dof)d">
%(connectivity)s
</DataItem>
</Topology>
''' % { "num_cell" : num_cell,
"cell_type" : vtk_types[degree-1][geodim-1],
"cell_dof" : cell_dofs[degree-1][geodim-1],
"dataitem_type" : dataitem_type,
"connectivity" : connectivity_data }
return dedent(topology_data)
def xdmf_geometry(coordinates, num_points) :
dataitem_type, geometry_data = xdmf_dataitem_dump(coordinates, '%f')
geometry_data = '''\
<Geometry GeometryType="XYZ">
<DataItem Format="%(dataitem_type)s" DataType="Float" Dimensions="%(num_points)d 3">
%(geometry_data)s
</DataItem>
</Geometry>
''' % { "dataitem_type" : dataitem_type,
"num_points" : num_points,
"geometry_data" : geometry_data }
return dedent(geometry_data)
def xdmf_attribute(name, values, rank, num_vals, center = "Node", symmetry = False) :
vtk_shapes = [ "Scalar", "Vector", "Tensor" ]
comp = pow(3, rank)
if symmetry and rank == 2:
comp /= 2
vtk_shapes[2] += "6"
if center not in [ "Node", "Cell" ] :
raise ValueError
dataitem_type, attribute_data = xdmf_dataitem_dump(values, '%f')
attribute_data = '''\
<Attribute Name="%(name)s" AttributeType="%(shape)s" Center="%(center)s">
<DataItem Format="%(dataitem_type)s" DataType="Float" Dimensions="%(num_vals)d %(comp)d">
%(attribute_data)s
</DataItem>
</Attribute>
''' % { "name" : name,
"shape" : vtk_shapes[rank],
"center" : center,
"dataitem_type" : dataitem_type,
"num_vals" : num_vals,
"comp" : comp,
"attribute_data" : attribute_data }
return dedent(attribute_data)
def write_xdmf(filename, gridname, functions) :
# various checks
meshes = set([ f.function_space().mesh().id() for f in functions.values() ])
if len(meshes) > 1 :
error("Functions don't share the same mesh!")
spacemap = lambda elm : "%s_%d" % (elm.family(), elm.degree())
basespaces = set([ spacemap(f.function_space().ufl_element()) for f in functions.values() ])
if "Lagrange_1" in basespaces and "Lagrange_2" in basespaces :
error("CG1 and CG2 cannot be exported on the same grid!")
degrees = [ f.function_space().ufl_element().degree() for f in functions.values() ]
ranks = [ f.rank() for f in functions.values() ]
if max(degrees) > 2 :
error("Cannot export more than quadratic functions!")
if max(ranks) > 2 :
error("Cannot export tensor fields of order greater than 2!")
# generating grids and exporting
with open(filename, "w") as ofile :
ofile.write('<?xml version="1.0"?>\n')
ofile.write('<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>\n')
ofile.write('<Xdmf Version="2.0" xmlns:xi="http://www.w3.org/2001/XInclude";>\n')
ofile.write('<Domain>\n')
ofile.write('<Grid Name="%s" GridType="Uniform">\n' % gridname)
ofile.write(generate_xdmf_grid(functions))
ofile.write('</Grid>\n')
ofile.write('</Domain>\n')
ofile.write('</Xdmf>\n')
def generate_xdmf_grid(functions) :
first_pass = True
for uname, u in functions.iteritems() :
fspace = u.function_space()
family = fspace.ufl_element().family()
degree = fspace.ufl_element().degree()
mesh = fspace.mesh()
dim = mesh.geometry().dim()
# rank indicates the tensorial shape of the function
rank = u.rank()
dofmap = fspace.dofmap()
# not necessary the same as the geometric dimension
num_dofs = dofmap.global_dimension()
num_block = dofmap.block_size
# FIXME not valid for DG0
if rank == 0 :
# scalar function
dofs = dofmap.dofs()
coords = fspace.dofmap().tabulate_all_coordinates(mesh).reshape((-1,dim))
elif rank in [ 1, 2 ] :
# vector/tensor function
dofs = np.empty(shape = (num_dofs / num_block, num_block), dtype = np.uintp)
for block in xrange(0, num_block) :
dofs[:,block] = fspace.sub(block).dofmap().dofs()
coords = fspace.dofmap().tabulate_all_coordinates(mesh).reshape((-1,dim))[dofs[:,0]]
else :
warning("'%s' of rank %d not supported, skipping ..." % (family, rank))
return None
# FIXME check match of coordinates between functions
coords_map = np.argsort(dofs.reshape((1, -1), order = 'F')[0,:])
# topology
ncell = mesh.num_cells()
dof_per_cell = [ [ 1, 3 ], [ 3, 6 ], [ 4, 10 ] ][dim-1][degree-1]
cell_dofs = np.zeros(shape=(ncell, dof_per_cell), dtype=np.uint)
if degree == 2 :
perm = np.array([ [ 0, 1, 2 ],
[ 0, 1, 2, 5, 3, 4 ],
[ 0, 1, 2, 3, 9, 6, 8, 7, 5, 4 ] ][dim-1])
else :
perm = np.arange(0, dim+1)
for c in xrange(0, ncell) :
cell_dofs[c,:] = coords_map[dofmap.cell_dofs(c)[perm]]
# reorder blockwise the dofs
u_vals = u.vector().array()[dofs]
# generating XML file
num_points = coords.shape[0]
zerovec = np.zeros((num_points, 3-dim))
if first_pass :
# always exports 3d coordinates
coords = np.concatenate((coords, zerovec), axis=1)
topology_sec = xdmf_topology(cell_dofs, cell_dofs.shape[0], degree, dim)
geometry_sec = xdmf_geometry(coords, coords.shape[0])
attributes_sec = ""
first_pass = False
# attributes always exported in 3d
if rank == 1 :
u_vals = np.concatenate((u_vals, zerovec), axis=1)
elif rank == 2 :
# FIXME symmetry
if dim == 2 :
u_vals = np.concatenate((u_vals[:, 0:2], zerovec,
u_vals[:, 2:4], zerovec,
zerovec, zerovec, zerovec), axis=1)
# FIXME dim = 1 and rank != 0
attributes_sec += xdmf_attribute(uname, u_vals, rank, coords.shape[0])
return topology_sec + geometry_sec + attributes_sec
def test_laplacian(mesh, order = 1) :
V = FunctionSpace(mesh, "CG", order)
bc = DirichletBC(V, Constant(0.0), lambda x, onbc : onbc)
u, v = TrialFunction(V), TestFunction(V)
f = Constant(1.0)
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx
u = Function(V)
solve(a == L, u, bc)
return u
def test_linearelasticity(mesh, order = 1) :
V = VectorFunctionSpace(mesh, "CG", order)
if V.cell().d == 3 :
u_rb = Constant((+0.3, 0.0, 0.0))
u_lb = Constant((-0.3, 0.0, 0.0))
f = Constant((0.0, 0.0, 0.0))
elif V.cell().d == 2 :
u_rb = Constant((+0.3, 0.0))
u_lb = Constant((-0.3, 0.0))
f = Constant((0.0, 0.0))
bc = [ DirichletBC(V, u_rb, lambda x, onbc : near(x[0], 1.0) and onbc),
DirichletBC(V, u_lb, lambda x, onbc : near(x[0], 0.0) and onbc) ]
E, nu = 1.0, 0.3
G, lmbda = E/(2*(1+nu)), E*nu/((1+nu)*(1-2*nu))
u, v = TrialFunction(V), TestFunction(V)
e = lambda u : sym(grad(u))
a = Constant(2*G) * inner(e(u), e(v)) * dx \
+ Constant(lmbda) * inner(div(u), div(v)) * dx
L = inner(f, v) * dx
u = Function(V)
solve(a == L, u, bc)
# "recover" stress
I = Identity(V.cell().d)
sigma_expr = Constant(2*G)*e(u) + Constant(lmbda)*div(u)*I
sigma = project(sigma_expr, TensorFunctionSpace(mesh, "CG", order))
return u, sigma
if __name__ == "__main__" :
ndiv = 8
mesh = UnitSquareMesh(ndiv, ndiv, "crossed")
u_el, sigma_el = test_linearelasticity(mesh, 2)
u_lp = test_laplacian(mesh, 2)
write_xdmf("mixed.xdmf", "Cube", { "Displacement" : u_el, "Stress" : sigma_el, "Potential" : u_lp })
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
