Source code for sfepy.tests.test_mesh_interp

import os.path as op

import numpy as nm

from sfepy.base.conf import transform_variables
import sfepy.base.testing as tst

variables = {
    'u'       : ('unknown field', 'f', 0),
    'v'       : ('test field',    'f', 'u'),
}

[docs] def in_dir(adir): return lambda x: op.join(adir, x)
[docs] def gen_datas(meshes): datas = {} for key, mesh in meshes.items(): bbox = mesh.get_bounding_box() nx = bbox[1,0] - bbox[0,0] centre = 0.5 * bbox.sum(axis=0) mesh.coors[:] -= centre data = nm.sin(4.0 * nm.pi * mesh.coors[:,0:1] / nx) datas['scalar_' + key] = data data = nm.zeros_like(mesh.coors) data[:,0] = 0.05 * nx * nm.sin(4.0 * nm.pi * mesh.coors[:,0] / nx) data[:,2] = 0.05 * nx * nm.cos(4.0 * nm.pi * mesh.coors[:,0] / nx) datas['vector_' + key] = data return datas
[docs] def do_interpolation(m2, m1, data, field_name, force=False): """Interpolate data from m1 to m2. """ from sfepy.discrete import Variables from sfepy.discrete.fem import FEDomain, Field fields = { 'scalar_si' : ((1,1), 'Omega', 2), 'vector_si' : ((3,1), 'Omega', 2), 'scalar_tp' : ((1,1), 'Omega', 1), 'vector_tp' : ((3,1), 'Omega', 1), } d1 = FEDomain('d1', m1) d1.create_region('Omega', 'all') f = fields[field_name] field1 = Field.from_args('f', nm.float64, f[0], d1.regions[f[1]], approx_order=f[2]) ff = {field1.name : field1} vv = Variables.from_conf(transform_variables(variables), ff) u1 = vv['u'] u1.set_from_mesh_vertices(data) d2 = FEDomain('d2', m2) d2.create_region('Omega', 'all') field2 = Field.from_args('f', nm.float64, f[0], d2.regions[f[1]], approx_order=f[2]) ff2 = {field2.name : field2} vv2 = Variables.from_conf(transform_variables(variables), ff2) u2 = vv2['u'] if not force: # Performs interpolation, if other field differs from self.field # or, in particular, is defined on a different mesh. u2.set_from_other(u1, strategy='interpolation', close_limit=0.5) else: coors = u2.field.get_coor() vals = u1.evaluate_at(coors, close_limit=0.5) u2.set_data(vals) return u1, u2
[docs] def prepare_variable(filename, n_components): from sfepy.discrete import FieldVariable from sfepy.discrete.fem import Mesh, FEDomain, Field mesh = Mesh.from_file(filename) bbox = mesh.get_bounding_box() dd = bbox[1,:] - bbox[0,:] data = (nm.sin(4.0 * nm.pi * mesh.coors[:,0:1] / dd[0]) * nm.cos(4.0 * nm.pi * mesh.coors[:,1:2] / dd[1])) domain = FEDomain('domain', mesh) omega = domain.create_region('Omega', 'all') field = Field.from_args('field', nm.float64, n_components, omega, approx_order=2) u = FieldVariable('u', 'parameter', field, primary_var_name='(set-to-None)') u.set_from_mesh_vertices(data * nm.arange(1, n_components + 1)[None, :]) return u
[docs] def test_interpolation(output_dir): from sfepy import data_dir from sfepy.discrete.fem import Mesh from sfepy.linalg import make_axis_rotation_matrix fname = in_dir(output_dir) meshes = { 'tp' : Mesh.from_file(data_dir + '/meshes/3d/block.mesh'), 'si' : Mesh.from_file(data_dir + '/meshes/3d/cylinder.mesh'), } datas = gen_datas(meshes) for field_name in ['scalar_si', 'vector_si', 'scalar_tp', 'vector_tp']: m1 = meshes[field_name[-2:]] for ia, angle in enumerate(nm.linspace(0.0, nm.pi, 11)): tst.report('%s: %d. angle: %f' % (field_name, ia, angle)) mtx = make_axis_rotation_matrix([0, 1, 0], angle) m2 = m1.copy('rotated mesh') m2.transform_coors(mtx) data = datas[field_name] u1, u2 = do_interpolation(m2, m1, data, field_name) if ia == 0: u1.save_as_mesh(fname('test_mesh_interp_%s_u1.vtk' % field_name)) u2.save_as_mesh(fname('test_mesh_interp_%s_u2.%03d.vtk' % (field_name, ia)))
[docs] def test_interpolation_two_meshes(output_dir): from sfepy import data_dir from sfepy.discrete import Variables from sfepy.discrete.fem import Mesh, FEDomain, Field m1 = Mesh.from_file(data_dir + '/meshes/3d/block.mesh') m2 = Mesh.from_file(data_dir + '/meshes/3d/cube_medium_tetra.mesh') m2.coors[:] *= 2.0 bbox = m1.get_bounding_box() dd = bbox[1,:] - bbox[0,:] data = nm.sin(4.0 * nm.pi * m1.coors[:,0:1] / dd[0]) \ * nm.cos(4.0 * nm.pi * m1.coors[:,1:2] / dd[1]) variables1 = { 'u' : ('unknown field', 'scalar_tp', 0), 'v' : ('test field', 'scalar_tp', 'u'), } variables2 = { 'u' : ('unknown field', 'scalar_si', 0), 'v' : ('test field', 'scalar_si', 'u'), } d1 = FEDomain('d1', m1) omega1 = d1.create_region('Omega', 'all') field1 = Field.from_args('scalar_tp', nm.float64, (1,1), omega1, approx_order=1) ff1 = {field1.name : field1} d2 = FEDomain('d2', m2) omega2 = d2.create_region('Omega', 'all') field2 = Field.from_args('scalar_si', nm.float64, (1,1), omega2, approx_order=0) ff2 = {field2.name : field2} vv1 = Variables.from_conf(transform_variables(variables1), ff1) u1 = vv1['u'] u1.set_from_mesh_vertices(data) vv2 = Variables.from_conf(transform_variables(variables2), ff2) u2 = vv2['u'] # Performs interpolation, if other field differs from self.field # or, in particular, is defined on a different mesh. u2.set_from_other(u1, strategy='interpolation', close_limit=0.1) fname = in_dir(output_dir) u1.save_as_mesh(fname('test_mesh_interp_block_scalar.vtk')) u2.save_as_mesh(fname('test_mesh_interp_cube_scalar.vtk'))
[docs] def test_invariance(): from sfepy import data_dir from sfepy.discrete.fem import Mesh meshes = { 'tp' : Mesh.from_file(data_dir + '/meshes/3d/block.mesh'), 'si' : Mesh.from_file(data_dir + '/meshes/3d/cylinder.mesh'), } datas = gen_datas(meshes) ok = True for field_name in ['scalar_si', 'vector_si', 'scalar_tp', 'vector_tp']: m1 = meshes[field_name[-2:]] data = datas[field_name] u1, u2 = do_interpolation(m1, m1, data, field_name, force=True) tst.report('max. difference:', nm.abs(u1() - u2()).max()) _ok = nm.allclose(u1(), u2(), rtol=0.0, atol=1e-12) tst.report('invariance for %s field: %s' % (field_name, _ok)) ok = ok and _ok assert ok
[docs] def test_invariance_qp(): from sfepy import data_dir from sfepy.discrete import Integral from sfepy.terms import Term from sfepy.discrete.common.mappings import get_physical_qps ok = True for name in ['meshes/3d/block.mesh', 'meshes/3d/cylinder.mesh', 'meshes/2d/square_quad.mesh', 'meshes/2d/square_unit_tri.mesh']: tst.report(name) u = prepare_variable(op.join(data_dir, name), n_components=3) omega = u.field.region integral = Integral('i', order=3) qps = get_physical_qps(omega, integral) coors = qps.values term = Term.new('ev_integrate(u)', integral, omega, u=u) term.setup() val1 = term.evaluate(mode='qp') val1 = val1.ravel() val2 = u.evaluate_at(coors).ravel() tst.report('value: max. difference:', nm.abs(val1 - val2).max()) ok1 = nm.allclose(val1, val2, rtol=0.0, atol=1e-12) tst.report('->', ok1) term = Term.new('ev_grad(u)', integral, omega, u=u) term.setup() val1 = term.evaluate(mode='qp') val1 = val1.transpose((0, 1, 3, 2)).ravel() val2 = u.evaluate_at(coors, mode='grad').ravel() tst.report('gradient: max. difference:', nm.abs(val1 - val2).max()) ok2 = nm.allclose(val1, val2, rtol=0.0, atol=1e-10) tst.report('->', ok2) ok = ok and ok1 and ok2 assert ok
[docs] def test_field_gradient(): from sfepy import data_dir ok = True for name in ['meshes/3d/block.mesh', 'meshes/3d/cylinder.mesh', 'meshes/2d/square_quad.mesh', 'meshes/2d/square_unit_tri.mesh']: tst.report(name) u = prepare_variable(op.join(data_dir, name), n_components=5) bbox = u.field.domain.get_mesh_bounding_box() coors = nm.c_[tuple([nm.linspace(ii[0] + 1e-3, ii[1] - 1e-3, 100) for ii in bbox.T])] grad, cells, status = u.evaluate_at(coors, mode='grad', close_limit=0.0, ret_status=True) agrad = nm.dot(grad[:, :, :], nm.ones((grad.shape[2], 1)))[..., 0] eps = 1e-5 val0 = u.evaluate_at(coors - eps, close_limit=0.0)[..., 0] val1 = u.evaluate_at(coors + eps, close_limit=0.0)[..., 0] ngrad = 0.5 * (val1 - val0) / eps ii = nm.where(status == 0) tst.report('max. difference:', nm.abs(agrad[ii] - ngrad[ii]).max()) _ok = nm.allclose(agrad[ii], ngrad[ii], rtol=0.0, atol=10 * eps) tst.report('->', _ok) ok = ok and _ok for ic in range(1, u.n_components): _ok = nm.allclose((ic + 1) * agrad[ii, 0], agrad[ii, ic], rtol=0.0, atol=1e-12) tst.report('component %d / component 0: mean: %.2f' % (ic, (agrad[ii, ic] / agrad[ii, 0]).mean())) tst.report('->', _ok) ok = ok and _ok assert ok
[docs] def test_evaluate_at(): from sfepy import data_dir from sfepy.discrete.fem import Mesh from sfepy.discrete import Variables from sfepy.discrete.fem import FEDomain, Field meshes = { 'tp' : Mesh.from_file(data_dir + '/meshes/3d/block.mesh'), } datas = gen_datas(meshes) fields = { 'scalar_tp' : ((1,1), 'Omega', 1), 'vector_tp' : ((3,1), 'Omega', 1), } ok = True for field_name in ['scalar_tp', 'vector_tp']: d = FEDomain('d', meshes['tp']) d.create_region('Omega', 'all') f = fields[field_name] field = Field.from_args('f', nm.complex128, f[0], d.regions[f[1]], approx_order=f[2]) ff = {field.name : field} vv = Variables.from_conf(transform_variables(variables), ff) u = vv['u'] bbox = d.get_mesh_bounding_box() t = nm.expand_dims(nm.linspace(0, 1, 100), 1) coors = nm.expand_dims(bbox[1] - bbox[0], 0) * t + bbox[0] data_r = datas[field_name] data_i = 2. / (1 + datas[field_name]) u.set_from_mesh_vertices(data_r) vals_r = u.evaluate_at(coors) u.set_from_mesh_vertices(data_i) vals_i = u.evaluate_at(coors) u.set_from_mesh_vertices(data_r + data_i * 1j) vals = u.evaluate_at(coors) _ok = nm.allclose(vals_r + vals_i * 1j, vals, rtol=0.0, atol=1e-12) _ok = _ok and nm.abs(vals).sum() > 1 tst.report('evaluating complex field %s: %s' % (field_name, _ok)) ok = ok and _ok assert ok