''' define interpolator from mesh field '''
from warnings import warn
import numpy as np
from matplotlib.tri import Triangulation
from pyxfem.core.config import npfloat
from pyxfem.core.mapping import MappingT3
[docs]
def get_triangulation(xy, t2v):
''' return a matplotlib.tri.Triangulation from
- a set of point xy and triangle connectivity t2xy '''
return Triangulation(xy[:,0], xy[:,1], t2v)
[docs]
def get_trifinder( xy, t2v):
''' return a matplotlib "trifinder " : a function that return the triangle number in
which each point it is called on is in '''
tri = get_triangulation(xy, t2v)
trif = tri.get_trifinder()
def finder(xyp):
return trif(xyp[...,0], xyp[...,1])
return finder
[docs]
class ConstantPerTriInterpolator():
''' Constant per Element field interpolator : find in which element is a point and return
the corresponding element consant value'''
def __init__(self, xy, t2v, values, trifinder = None):
''' values are the per element values, mus also pass either a trifinder objet or a sMesh'''
if trifinder is None:
self.trifinder = get_trifinder(xy, t2v)
else :
self.trifinder = trifinder
self.values = values
def __call__(self, xy):
elems = self.trifinder(xy)
return self.values[elems]
[docs]
class LinearPerTriInterpolator():
''' linear per Element field interpolator field is supposed C0 on the input mesh,
and linear per element input : value at node and mesh.
'''
def __init__(self, xy, t2v, values, trifinder = None):
''' values are the per node values, must also pass either a trifinder objet or a sMesh'''
if trifinder is None:
self.trifinder = get_trifinder(xy, t2v)
else :
self.trifinder = trifinder
self.t2v = t2v
t3xy = xy[t2v]
self.map = MappingT3(t3xy)
self.values = values
assert(values.shape[0] == xy.shape[0])
self.fieldshape = values.shape[1:]
def __call__(self, xy):
xy = np.atleast_2d(xy)
elems = np.asarray(self.trifinder(xy), dtype =int)
xy0 = self.map.xy[elems, 0, :]
xy0s = xy - xy0
uv = np.einsum('...jl,...l->...j', self.map.inv_F[elems], xy0s)
valelems = self.values[self.t2v[elems]]
u = uv[...,0]
v = uv[...,1]
shape = u.shape + tuple(1 for i in range(len(self.fieldshape)))
u = u.reshape(shape)
v = v.reshape(shape)
nodeaxis = -1-len(self.fieldshape)
valelemsv0 = valelems.take(0, axis = nodeaxis )
valelemsv1 = valelems.take(1, axis = nodeaxis )
valelemsv2 = valelems.take(2, axis = nodeaxis )
val = valelemsv0 + (valelemsv1-valelemsv0)*u + (valelemsv2-valelemsv0)*v
if np.any(elems < 0) :
warn("some point are outside of the mesh linearPerTriInterpolator," \
" corresponding values are set to np.inf")
val[elems < 0] = np.inf
return np.asarray(val, dtype = npfloat)
if __name__ == '__main__':
import gmsh
import pyxfem.core.gmsh_mesh_generator as mg
from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh
import pyxfem.postpro.gmsh_post_pro_interface as gp
gmsh.initialize()
gmsh.option.setNumber("PostProcessing.AnimationCycle", 1)
gmsh.option.setNumber('General.Terminal', 0)
#model_name1, ____ = mg.square(length=1., relh=0.1, transfinite = True)
model_name1, ____ = mg.crack_in_square_contour(relh = 0.1)
gmsh.fltk.run()
mesh1 = gmshmodel2smesh(model_name1)
model_name2, ____ = mg.l_shape(relh=0.05)
mesh2 = gmshmodel2smesh(model_name2)
def fun_scal(xy):
''' a simple scalar function of coordinate '''
return xy[...,0]**2
def fun_vec(xy):
''' a simple vector function of coordinate '''
return xy
fun = fun_scal
mesh1_trixy = mesh1.t2v_coord()
gp.listPlotFieldTri(mesh1_trixy, fun(mesh1_trixy),
viewname='base scalar field on mesh1', Range =[0.,0.25])
interp = LinearPerTriInterpolator(mesh1.xy, mesh1.t2v, fun(mesh1.xy))
val2 = interp(mesh2.xy)
gp.listPlotFieldTri(mesh2.t2v_coord(), val2[np.asarray(mesh2.t2v, dtype=int)],
viewname='interp scalar field on mesh2', Range =[0.,0.25])
fun = fun_vec
mesh1_trixy = mesh1.t2v_coord()
gp.listPlotFieldTri(mesh1_trixy, fun(mesh1_trixy),
viewname='base vector field on mesh1', Range = [0., np.sqrt(2)])
interp = LinearPerTriInterpolator(mesh1.xy, mesh1.t2v, fun(mesh1.xy), trifinder=interp.trifinder)
val2 = interp(mesh2.xy)
gp.listPlotFieldTri(mesh2.t2v_coord(), val2[np.asarray(mesh2.t2v, dtype=int)],
viewname='interp vector field on mesh2', Range = [0., np.sqrt(2)])
gmsh.fltk.run()