"""
Global interpolation functions.
"""
import numpy as nm
from sfepy.base.base import assert_, output, get_default_attr
from sfepy.base.timing import Timer
from sfepy.discrete.fem.geometry_element import create_geometry_elements
import sfepy.discrete.common.extmods.crefcoors as crc
[docs]def get_ref_coors_convex(field, coors, close_limit=0.1, cache=None,
verbose=False):
"""
Get reference element coordinates and elements corresponding to given
physical coordinates.
Parameters
----------
field : Field instance
The field defining the approximation.
coors : array
The physical coordinates.
close_limit : float, optional
The maximum limit distance of a point from the closest
element allowed for extrapolation.
cache : Struct, optional
To speed up a sequence of evaluations, the field mesh and other data
can be cached. Optionally, the cache can also contain the reference
element coordinates as `cache.ref_coors`, `cache.cells` and
`cache.status`, if the evaluation occurs in the same coordinates
repeatedly. In that case the mesh related data are ignored.
verbose : bool
If False, reduce verbosity.
Returns
-------
ref_coors : array
The reference coordinates.
cells : array
The cell indices corresponding to the reference coordinates.
status : array
The status: 0 is success, 1 is extrapolation within `close_limit`, 2 is
extrapolation outside `close_limit`, 3 is failure, 4 is failure due to
non-convergence of the Newton iteration in tensor product cells.
Notes
-----
Outline of the algorithm for finding xi such that X(xi) = P:
1. make inverse connectivity - for each vertex have cells it is in.
2. find the closest vertex V.
3. choose initial cell: i0 = first from cells incident to V.
4. while not P in C_i, change C_i towards P, check if P in new C_i.
"""
timer = Timer()
ref_coors = get_default_attr(cache, 'ref_coors', None)
if ref_coors is None:
extrapolate = close_limit > 0.0
ref_coors = nm.empty_like(coors)
cells = nm.empty((coors.shape[0],), dtype=nm.int32)
status = nm.empty((coors.shape[0],), dtype=nm.int32)
cmesh = get_default_attr(cache, 'cmesh', None)
if cmesh is None:
timer.start()
mesh = field.create_mesh(extra_nodes=False)
cmesh = mesh.cmesh
gels = create_geometry_elements()
cmesh.set_local_entities(gels)
cmesh.setup_entities()
centroids = cmesh.get_centroids(cmesh.tdim)
if field.gel.name != '3_8':
normals0 = cmesh.get_facet_normals()
normals1 = None
else:
normals0 = cmesh.get_facet_normals(0)
normals1 = cmesh.get_facet_normals(1)
output('cmesh setup: %f s' % timer.stop(), verbose=verbose)
else:
centroids = cache.centroids
normals0 = cache.normals0
normals1 = cache.normals1
kdtree = get_default_attr(cache, 'kdtree', None)
if kdtree is None:
from scipy.spatial import cKDTree as KDTree
timer.start()
kdtree = KDTree(cmesh.coors)
output('kdtree: %f s' % timer.stop(), verbose=verbose)
timer.start()
ics = kdtree.query(coors)[1]
output('kdtree query: %f s' % timer.stop(), verbose=verbose)
ics = nm.asarray(ics, dtype=nm.int32)
coors = nm.ascontiguousarray(coors)
ctx = field.create_basis_context()
timer.start()
crc.find_ref_coors_convex(ref_coors, cells, status, coors, cmesh,
centroids, normals0, normals1, ics,
extrapolate, 1e-15, close_limit, ctx)
output('ref. coordinates: %f s' % timer.stop(), verbose=verbose)
else:
cells = cache.cells
status = cache.status
return ref_coors, cells, status
[docs]def get_potential_cells(coors, cmesh, centroids=None, extrapolate=True):
"""
Get cells that potentially contain points with the given physical
coordinates.
Parameters
----------
coors : array
The physical coordinates.
cmesh : CMesh instance
The cmesh defining the cells.
centroids : array, optional
The centroids of the cells.
extrapolate : bool
If True, even the points that are surely outside of the
cmesh are considered and assigned potential cells.
Returns
-------
potential_cells : array
The indices of the cells that potentially contain the points.
offsets : array
The offsets into `potential_cells` for each point: a point ``ip`` is
potentially in cells ``potential_cells[offsets[ip]:offsets[ip+1]]``.
"""
from scipy.spatial import cKDTree as KDTree
if centroids is None:
centroids = cmesh.get_centroids(cmesh.tdim)
kdtree = KDTree(coors)
conn = cmesh.get_cell_conn()
cc = conn.indices.reshape(cmesh.n_el, -1)
cell_coors = cmesh.coors[cc]
rays = cell_coors - centroids[:, None]
radii = nm.linalg.norm(rays, ord=nm.inf, axis=2).max(axis=1)
potential_cells = [[]] * coors.shape[0]
for ic, centroid in enumerate(centroids):
ips = kdtree.query_ball_point(centroid, radii[ic], p=nm.inf)
if len(ips):
for ip in ips:
if not len(potential_cells[ip]):
potential_cells[ip] = []
potential_cells[ip].append(ic)
lens = nm.array([0] + [len(ii) for ii in potential_cells], dtype=nm.int32)
if extrapolate:
# Deal with the points outside of the field domain - insert elements
# incident to the closest mesh vertex.
iin = nm.where(lens[1:] == 0)[0]
if len(iin):
kdtree = KDTree(cmesh.coors)
ics = kdtree.query(coors[iin])[1]
cmesh.setup_connectivity(0, cmesh.tdim)
conn = cmesh.get_conn(0, cmesh.tdim)
oo = conn.offsets
for ii, ip in enumerate(iin):
ik = ics[ii]
potential_cells[ip] = conn.indices[oo[ik]:oo[ik+1]]
lens[ip+1] = len(potential_cells[ip])
offsets = nm.cumsum(lens, dtype=nm.int32)
potential_cells = nm.concatenate(potential_cells).astype(nm.int32)
return potential_cells, offsets
[docs]def get_ref_coors_general(field, coors, close_limit=0.1, get_cells_fun=None,
cache=None, verbose=False):
"""
Get reference element coordinates and elements corresponding to given
physical coordinates.
Parameters
----------
field : Field instance
The field defining the approximation.
coors : array
The physical coordinates.
close_limit : float, optional
The maximum limit distance of a point from the closest
element allowed for extrapolation.
get_cells_fun : callable, optional
If given, a function with signature ``get_cells_fun(coors, cmesh,
**kwargs)`` returning cells and offsets that potentially contain points
with the coordinates `coors`. When not given,
:func:`get_potential_cells()` is used.
cache : Struct, optional
To speed up a sequence of evaluations, the field mesh and other data
can be cached. Optionally, the cache can also contain the reference
element coordinates as `cache.ref_coors`, `cache.cells` and
`cache.status`, if the evaluation occurs in the same coordinates
repeatedly. In that case the mesh related data are ignored.
verbose : bool
If False, reduce verbosity.
Returns
-------
ref_coors : array
The reference coordinates.
cells : array
The cell indices corresponding to the reference coordinates.
status : array
The status: 0 is success, 1 is extrapolation within `close_limit`, 2 is
extrapolation outside `close_limit`, 3 is failure, 4 is failure due to
non-convergence of the Newton iteration in tensor product cells. If
close_limit is 0, then status 5 indicates points outside of the field
domain that had no potential cells.
"""
timer = Timer()
ref_coors = get_default_attr(cache, 'ref_coors', None)
if ref_coors is None:
extrapolate = close_limit > 0.0
get = get_potential_cells if get_cells_fun is None else get_cells_fun
ref_coors = nm.empty_like(coors)
cells = nm.empty((coors.shape[0],), dtype=nm.int32)
status = nm.empty((coors.shape[0],), dtype=nm.int32)
cmesh = get_default_attr(cache, 'cmesh', None)
if cmesh is None:
timer.start()
mesh = field.create_mesh(extra_nodes=False)
cmesh = mesh.cmesh
if get_cells_fun is None:
centroids = cmesh.get_centroids(cmesh.tdim)
else:
centroids = None
output('cmesh setup: %f s' % timer.stop(), verbose=verbose)
else:
centroids = cache.centroids
timer.start()
potential_cells, offsets = get(coors, cmesh, centroids=centroids,
extrapolate=extrapolate)
output('potential cells: %f s' % timer.stop(), verbose=verbose)
coors = nm.ascontiguousarray(coors)
ctx = field.create_basis_context()
eval_cmesh = get_default_attr(cache, 'eval_cmesh', None)
if eval_cmesh is None:
timer.start()
mesh = field.create_eval_mesh()
if mesh is None:
eval_cmesh = cmesh
else:
eval_cmesh = mesh.cmesh
output('eval_cmesh setup: %f s'
% timer.stop(), verbose=verbose)
timer.start()
crc.find_ref_coors(ref_coors, cells, status, coors, eval_cmesh,
potential_cells, offsets, extrapolate,
1e-15, close_limit, ctx)
if extrapolate:
assert_(nm.all(status < 5))
output('ref. coordinates: %f s' % timer.stop(), verbose=verbose)
else:
cells = cache.cells
status = cache.status
return ref_coors, cells, status
[docs]def get_ref_coors(field, coors, strategy='general', close_limit=0.1,
get_cells_fun=None, cache=None, verbose=False):
"""
Get reference element coordinates and elements corresponding to given
physical coordinates.
Parameters
----------
field : Field instance
The field defining the approximation.
coors : array
The physical coordinates.
strategy : {'general', 'convex'}, optional
The strategy for finding the elements that contain the coordinates. For
convex meshes, the 'convex' strategy might be faster than the 'general'
one.
close_limit : float, optional
The maximum limit distance of a point from the closest
element allowed for extrapolation.
get_cells_fun : callable, optional
If given, a function with signature ``get_cells_fun(coors, cmesh,
**kwargs)`` returning cells and offsets that potentially contain points
with the coordinates `coors`. Applicable only when `strategy` is
'general'. When not given, :func:`get_potential_cells()` is used.
cache : Struct, optional
To speed up a sequence of evaluations, the field mesh and other data
can be cached. Optionally, the cache can also contain the reference
element coordinates as `cache.ref_coors`, `cache.cells` and
`cache.status`, if the evaluation occurs in the same coordinates
repeatedly. In that case the mesh related data are ignored.
verbose : bool
If False, reduce verbosity.
Returns
-------
ref_coors : array
The reference coordinates.
cells : array
The cell indices corresponding to the reference coordinates.
status : array
The status: 0 is success, 1 is extrapolation within `close_limit`, 2 is
extrapolation outside `close_limit`, 3 is failure, 4 is failure due to
non-convergence of the Newton iteration in tensor product cells. If
close_limit is 0, then for the 'general' strategy the status 5
indicates points outside of the field domain that had no potential
cells.
"""
if strategy == 'general':
return get_ref_coors_general(field, coors, close_limit=close_limit,
get_cells_fun=get_cells_fun,
cache=cache, verbose=verbose)
elif strategy == 'convex':
return get_ref_coors_convex(field, coors, close_limit=close_limit,
cache=cache, verbose=verbose)
else:
raise ValueError('unsupported strategy! (%s)' % strategy)