def create_mesh(self, extra_nodes=True):
"""
Create a mesh from the field region, optionally including the field
extra nodes.
"""
mesh = self.domain.mesh
if self.approx_order != 0:
conns, mat_ids, descs = [], [], []
for ig, ap in self.aps.iteritems():
group = self.domain.groups[ig]
if extra_nodes:
conn = ap.econn
else:
offset = group.shape.n_ep
conn = ap.econn[:,:offset]
conns.append(conn)
mat_ids.append(mesh.mat_ids[ig])
descs.append(mesh.descs[ig])
if extra_nodes:
coors = self.coors
else:
coors = self.coors[:self.n_vertex_dof]
mesh = Mesh.from_data(self.name, coors, None, conns,
mat_ids, descs)
return mesh
def gen_block_mesh(dims, shape, centre, mat_id=0, name='block',
coors=None, verbose=True):
"""
Generate a 2D or 3D block mesh. The dimension is determined by the
lenght of the shape argument.
Parameters
----------
dims : array of 2 or 3 floats
Dimensions of the block.
shape : array of 2 or 3 ints
Shape (counts of nodes in x, y, z) of the block mesh.
centre : array of 2 or 3 floats
Centre of the block.
mat_id : int, optional
The material id of all elements.
name : string
Mesh name.
verbose : bool
If True, show progress of the mesh generation.
Returns
-------
mesh : Mesh instance
"""
dims = nm.asarray(dims, dtype=nm.float64)
shape = nm.asarray(shape, dtype=nm.int32)
centre = nm.asarray(centre, dtype=nm.float64)
dim = shape.shape[0]
centre = centre[:dim]
dims = dims[:dim]
n_nod = nm.prod(shape)
output('generating %d vertices...' % n_nod, verbose=verbose)
x0 = centre - 0.5 * dims
dd = dims / (shape - 1)
ngrid = nm.mgrid[[slice(ii) for ii in shape]]
ngrid.shape = (dim, n_nod)
coors = x0 + ngrid.T * dd
output('...done', verbose=verbose)
n_el = nm.prod(shape - 1)
output('generating %d cells...' % n_el, verbose=verbose)
mat_ids = nm.empty((n_el,), dtype=nm.int32)
mat_ids.fill(mat_id)
conn, desc = get_tensor_product_conn(shape)
output('...done', verbose=verbose)
mesh = Mesh.from_data(name, coors, None, [conn], [mat_ids], [desc])
return mesh
def gen_cylinder_mesh(dims, shape, centre, axis='x', force_hollow=False,
is_open=False, open_angle=0.0, non_uniform=False,
name='cylinder', verbose=True):
"""
Generate a cylindrical mesh along an axis. Its cross-section can be
ellipsoidal.
Parameters
----------
dims : array of 5 floats
Dimensions of the cylinder: inner surface semi-axes a1, b1, outer
surface semi-axes a2, b2, length.
shape : array of 3 ints
Shape (counts of nodes in radial, circumferential and longitudinal
directions) of the cylinder mesh.
centre : array of 3 floats
Centre of the cylinder.
axis: one of 'x', 'y', 'z'
The axis of the cylinder.
force_hollow : boolean
Force hollow mesh even if inner radii a1 = b1 = 0.
is_open : boolean
Generate an open cylinder segment.
open_angle : float
Opening angle in radians.
non_uniform : boolean
If True, space the mesh nodes in radial direction so that the element
volumes are (approximately) the same, making thus the elements towards
the outer surface thinner.
name : string
Mesh name.
verbose : bool
If True, show progress of the mesh generation.
Returns
-------
mesh : Mesh instance
"""
dims = nm.asarray(dims, dtype=nm.float64)
shape = nm.asarray(shape, dtype=nm.int32)
centre = nm.asarray(centre, dtype=nm.float64)
a1, b1, a2, b2, length = dims
nr, nfi, nl = shape
origin = centre - nm.array([0.5 * length, 0.0, 0.0])
dfi = 2.0 * (nm.pi - open_angle) / nfi
if is_open:
nnfi = nfi + 1
else:
nnfi = nfi
is_hollow = force_hollow or not (max(abs(a1), abs(b1)) < 1e-15)
if is_hollow:
mr = 0
else:
mr = (nnfi - 1) * nl
grid = nm.zeros((nr, nnfi, nl), dtype=nm.int32)
n_nod = nr * nnfi * nl - mr
coors = nm.zeros((n_nod, 3), dtype=nm.float64)
angles = nm.linspace(open_angle, open_angle+(nfi)*dfi, nfi+1)
xs = nm.linspace(0.0, length, nl)
if non_uniform:
ras = nm.zeros((nr,), dtype=nm.float64)
rbs = nm.zeros_like(ras)
advol = (a2**2 - a1**2) / (nr - 1)
bdvol = (b2**2 - b1**2) / (nr - 1)
ras[0], rbs[0] = a1, b1
for ii in range(1, nr):
ras[ii] = nm.sqrt(advol + ras[ii-1]**2)
rbs[ii] = nm.sqrt(bdvol + rbs[ii-1]**2)
else:
ras = nm.linspace(a1, a2, nr)
rbs = nm.linspace(b1, b2, nr)
# This is 3D only...
output('generating %d vertices...' % n_nod, verbose=verbose)
ii = 0
for ix in range(nr):
a, b = ras[ix], rbs[ix]
for iy, fi in enumerate(angles[:nnfi]):
for iz, x in enumerate(xs):
grid[ix,iy,iz] = ii
coors[ii] = origin + [x, a * nm.cos(fi), b * nm.sin(fi)]
ii += 1
if not is_hollow and (ix == 0):
if iy > 0:
grid[ix,iy,iz] = grid[ix,0,iz]
ii -= 1
assert_(ii == n_nod)
output('...done', verbose=verbose)
n_el = (nr - 1) * nnfi * (nl - 1)
conn = nm.zeros((n_el, 8), dtype=nm.int32)
#.........这里部分代码省略.........
def dump_to_vtk(filename, output_filename_trunk=None, step0=0, steps=None,
fields=None, linearization=None):
"""Dump a multi-time-step results file into a sequence of VTK files."""
def _save_step(suffix, out, mesh):
if linearization is not None:
output('linearizing...')
out = _linearize(out, fields, linearization)
output('...done')
for key, val in out.iteritems():
lmesh = val.get('mesh', mesh)
lmesh.write(output_filename_trunk + '_' + key + suffix,
io='auto', out={key : val})
if hasattr(val, 'levels'):
output('max. refinement per group:', val.levels)
else:
mesh.write(output_filename_trunk + suffix, io='auto', out=out)
output('dumping to VTK...')
io = MeshIO.any_from_filename(filename)
mesh = Mesh.from_file(filename, io=io)
if output_filename_trunk is None:
output_filename_trunk = get_trunk(filename)
try:
ts = TimeStepper(*io.read_time_stepper())
all_steps, times, nts, dts = extract_times(filename)
except ValueError:
output('no time stepping info found, assuming single step')
out = io.read_data(0)
if out is not None:
_save_step('.vtk', out, mesh)
ret = None
else:
ts.times = times
ts.n_step = times.shape[0]
if steps is None:
ii0 = nm.searchsorted(all_steps, step0)
iterator = ((all_steps[ii], times[ii])
for ii in xrange(ii0, len(times)))
else:
iterator = [(step, ts.times[step]) for step in steps]
max_step = all_steps.max()
for step, time in iterator:
output(ts.format % (step, max_step))
out = io.read_data(step)
if out is None: break
_save_step('.' + ts.suffix % step + '.vtk', out, mesh)
ret = ts.suffix
output('...done')
return ret
def create_expression_output(expression, name, primary_field_name,
fields, materials, variables,
functions=None, mode='eval', term_mode=None,
extra_args=None, verbose=True, kwargs=None,
min_level=0, max_level=1, eps=1e-4):
"""
Create output mesh and data for the expression using the adaptive
linearizer.
Parameters
----------
expression : str
The expression to evaluate.
name : str
The name of the data.
primary_field_name : str
The name of field that defines the element groups and polynomial
spaces.
fields : dict
The dictionary of fields used in `variables`.
materials : Materials instance
The materials used in the expression.
variables : Variables instance
The variables used in the expression.
functions : Functions instance, optional
The user functions for materials etc.
mode : one of 'eval', 'el_avg', 'qp'
The evaluation mode - 'qp' requests the values in quadrature points,
'el_avg' element averages and 'eval' means integration over
each term region.
term_mode : str
The term call mode - some terms support different call modes
and depending on the call mode different values are
returned.
extra_args : dict, optional
Extra arguments to be passed to terms in the expression.
verbose : bool
If False, reduce verbosity.
kwargs : dict, optional
The variables (dictionary of (variable name) : (Variable
instance)) to be used in the expression.
min_level : int
The minimum required level of mesh refinement.
max_level : int
The maximum level of mesh refinement.
eps : float
The relative tolerance parameter of mesh adaptivity.
Returns
-------
out : dict
The output dictionary.
"""
field = fields[primary_field_name]
vertex_coors = field.coors[:field.n_vertex_dof, :]
ps = field.poly_space
gps = field.gel.poly_space
vertex_conn = field.econn[:, :field.gel.n_vertex]
eval_dofs = get_eval_expression(expression,
fields, materials, variables,
functions=functions,
mode=mode, extra_args=extra_args,
verbose=verbose, kwargs=kwargs)
eval_coors = get_eval_coors(vertex_coors, vertex_conn, gps)
(level, coors, conn,
vdofs, mat_ids) = create_output(eval_dofs, eval_coors,
vertex_conn.shape[0], ps,
min_level=min_level,
max_level=max_level, eps=eps)
mesh = Mesh.from_data('linearized_mesh', coors, None, [conn], [mat_ids],
field.domain.mesh.descs)
out = {}
out[name] = Struct(name='output_data', mode='vertex',
data=vdofs, var_name=name, dofs=None,
mesh=mesh, level=level)
out = convert_complex_output(out)
return out
def extract_time_history(filename, extract, verbose=True):
"""Extract time history of a variable from a multi-time-step results file.
Parameters
----------
filename : str
The name of file to extract from.
extract : str
The description of what to extract in a string of comma-separated
description items. A description item consists of: name of the variable
to extract, mode ('e' for elements, 'n' for nodes), ids of the nodes or
elements (given by the mode). Example: 'u n 10 15, p e 0' means
variable 'u' in nodes 10, 15 and variable 'p' in element 0.
verbose : bool
Verbosity control.
Returns
-------
ths : dict
The time histories in a dict with variable names as keys. If a
nodal variable is requested in elements, its value is a dict of histories
in the element nodes.
ts : TimeStepper instance
The time stepping information.
"""
output('extracting selected data...', verbose=verbose)
output('selection:', extract, verbose=verbose)
##
# Parse extractions.
pes = OneTypeList(Struct)
for chunk in extract.split(','):
aux = chunk.strip().split()
pes.append(Struct(var = aux[0],
mode = aux[1],
indx = map(int, aux[2:]),
igs = None))
##
# Verify array limits, set igs for element data, shift indx.
mesh = Mesh.from_file(filename)
n_el, n_els, offs = mesh.n_el, mesh.n_els, mesh.el_offsets
for pe in pes:
if pe.mode == 'n':
for ii in pe.indx:
if (ii < 0) or (ii >= mesh.n_nod):
raise ValueError('node index 0 <= %d < %d!'
% (ii, mesh.n_nod))
if pe.mode == 'e':
pe.igs = []
for ii, ie in enumerate(pe.indx[:]):
if (ie < 0) or (ie >= n_el):
raise ValueError('element index 0 <= %d < %d!'
% (ie, n_el))
ig = (ie < n_els).argmax()
pe.igs.append(ig)
pe.indx[ii] = ie - offs[ig]
## print pes
##
# Extract data.
# Assumes only one element group (ignores igs)!
io = MeshIO.any_from_filename(filename)
ths = {}
for pe in pes:
mode, nname = io.read_data_header(pe.var)
output(mode, nname, verbose=verbose)
if ((pe.mode == 'n' and mode == 'vertex') or
(pe.mode == 'e' and mode == 'cell')):
th = io.read_time_history(nname, pe.indx)
elif pe.mode == 'e' and mode == 'vertex':
conn = mesh.conns[0]
th = {}
for iel in pe.indx:
ips = conn[iel]
th[iel] = io.read_time_history(nname, ips)
else:
raise ValueError('cannot extract cell data %s in nodes!' % pe.var)
ths[pe.var] = th
output('...done', verbose=verbose)
ts = TimeStepper(*io.read_time_stepper())
return ths, ts
客服电话
APP下载
官方微信
请发表评论