Source code for sfepy.tests.test_term_sensitivity

import gc
import numpy as nm
import pytest

from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy import data_dir
import sfepy.base.testing as tst

filename_mesh = data_dir + '/meshes/3d/special/cube_sphere.mesh'

fields = {
    'scalar_field' : ('real', 'scalar', 'Omega', 1),
    'vector_field' : ('real', 'vector', 'Omega', 1),
}

variables = {
    'U1': ('parameter field', 'vector_field', '(set-to-None)'),
    'U2': ('parameter field', 'vector_field', '(set-to-None)'),
    'P1': ('parameter field', 'scalar_field', '(set-to-None)'),
    'P2': ('parameter field', 'scalar_field', '(set-to-None)'),
    'V': ('parameter field', 'vector_field', '(set-to-None)'),
}

regions = {
    'Omega': 'all',
    'Omega1': 'cells of group 1',
    'Omega2': 'cells of group 2',
    'Interface': ('r.Omega1 *v r.Omega2', 'facet', 'Omega1'),
}

integrals = {
    'i': 2,
}

piezo_g = nm.array([[0, 0, 0, 0, 10, 0],
                    [0, 0, 0, 0, 0, 10],
                    [-5, -5, 15, 0, 0, 0]])

materials = {
    'mat': ({
        'D': {
            'Omega1': stiffness_from_youngpoisson(3, 100, 0.3),
            'Omega2': stiffness_from_youngpoisson(3, 1, 0.49),
        },
        'K': {'Omega1': nm.eye(3), 'Omega2': nm.eye(3) * 1e-3},
        'c': {'Omega1': 1, 'Omega2': 1000},
        'g': {'Omega1': piezo_g * 1e2, 'Omega2': piezo_g},
        's': {
            'Omega1': nm.array([[1, 2, 3, 0, 0, 0]]).T,
            'Omega2': nm.array([[8, 4, 2, 0, 0, 0]]).T,
        },
    },),
}

equations = {}

test_terms = [
    ('ev_sd_lin_elastic', 'dw_lin_elastic', 'Omega', 'mat.D', 'U1', 'U2'),
    ('de_sd_lin_elastic', 'dw_lin_elastic', 'Omega', 'mat.D', 'U1', 'U2'),
    ('ev_sd_dot', 'dw_dot', 'Omega', None, 'U1', 'U2'),
    ('de_sd_dot', 'de_dot', 'Omega', 'mat.K', 'U1', 'U2'),
    ('de_sd_dot', 'de_dot', 'Omega', 'mat.c', 'P1', 'P2'),
    ('ev_sd_diffusion', 'dw_diffusion', 'Omega', 'mat.K', 'P1', 'P2'),
    ('de_sd_diffusion', 'de_diffusion', 'Omega', 'mat.K', 'P1', 'P2'),
    ('ev_sd_div', 'dw_stokes', 'Omega', None, 'U1', 'P1'),
    ('de_sd_stokes', 'de_stokes', 'Omega', None, 'U1', 'P1'),
    ('ev_sd_div_grad', 'dw_div_grad', 'Omega', 'mat.c', 'U1', 'U2'),
    ('de_sd_div_grad', 'de_div_grad', 'Omega', 'mat.c', 'U1', 'U2'),
    ('ev_sd_piezo_coupling', 'dw_piezo_coupling', 'Omega', 'mat.g', 'U1', 'P1'),
    ('de_sd_piezo_coupling', 'dw_piezo_coupling', 'Omega', 'mat.g', 'U1', 'P1'),
    ('ev_sd_convect', 'de_convect', 'Omega', None, 'U1', 'U2'),
    ('ev_sd_surface_ltr', 'dw_surface_ltr', 'Interface', 'mat.s', 'U1', None),
    ('de_sd_surface_ltr', 'de_surface_ltr', 'Interface', 'mat.s', 'U1', None),
    ('de_sd_surface_ltr', 'de_surface_ltr', 'Interface', 'mat.c', 'U1', None),
    ('de_sd_v_dot_grad_s', 'dw_v_dot_grad_s', 'Omega', 'mat.c', 'U1', 'P1'),
]


[docs] def modify_mesh(val, spbox, dv_mode, cp_pos): cpoints0 = spbox.get_control_points().copy() for pts, dir in dv_mode: for pt in pts: spbox.move_control_point(cp_pos[pt], nm.array(dir) * val) new_coors = spbox.evaluate() spbox.set_control_points(cpoints0) return new_coors
[docs] @pytest.fixture(scope='module') def problem(): import sys from sfepy.discrete import Problem from sfepy.base.conf import ProblemConf conf = ProblemConf.from_dict(globals(), sys.modules[__name__]) problem = Problem.from_conf(conf, init_equations=False) return problem
[docs] @pytest.mark.slow def test_sensitivity(problem): from sfepy.discrete import Variables from sfepy.mesh.splinebox import SplineBox tolerance = 1e-4 ok = True variables = Variables.from_conf(problem.conf.variables, problem.fields) for var_name in variables.names: var = variables[var_name] n_dof = var.field.n_nod * var.field.shape[0] aux = nm.arange(n_dof, dtype=nm.float64) var.set_data(aux) mesh = problem.domain.mesh bbox = nm.array(mesh.get_bounding_box()).T spbox = SplineBox(bbox, mesh.coors) dvel_modes = [ # expand inner cylinder, no volume change [([20, 21, 22, 23], (-1, -1, 0)), ([24, 25, 26, 27], (-1, 1, 0)), ([36, 37, 38, 39], (1, -1, 0)), ([40, 41, 42, 43], (1, 1, 0))], # volume change [(range(16, 32), (0.2, 0, 0)), (range(32, 48), (0.4, 0, 0)), (range(48, 52), (0.6, 0.2, 0.2)), (range(52, 56), (0.8, 0.2, 0.3)), (range(56, 60), (1.0, 0.2, 0.4)), (range(60, 64), (1.2, 0.2, 0.5))], ] r4 = range(4) cp_pos = {i*16 + j*4 + k: (i, j, k) for k in r4 for j in r4 for i in r4} # compute design velocities dvels = [] for dv_mode in dvel_modes: dvel = 0 for pts, dir in dv_mode: for pt in pts: dvel += spbox.evaluate_derivative(cp_pos[pt], dir) dvels.append(dvel) for tname_sa, tname, rname, mat, var1, var2 in test_terms: args = [] if mat is None else [mat] args += [var1] if var2 is None else [var1, var2] term = '%s.i.%s(%s)' % (tname, rname, ', '.join(args)) term_sa = '%s.i.%s(%s)' % (tname_sa, rname, ', '.join(args + ['V'])) val = problem.evaluate(term, var_dict=variables.as_dict()) tst.report('%s: %s' % (tname, val)) dt = 1e-6 for ii, dvel in enumerate(dvels): val = problem.evaluate(term, var_dict=variables.as_dict()) variables['V'].set_data(dvel) val_sa = problem.evaluate(term_sa, var_dict=variables.as_dict()) tst.report('%s - mesh_velocity mode %d' % (tname_sa, ii)) # mesh perturbation + new_coors = modify_mesh(dt/2., spbox, dvel_modes[ii], cp_pos) problem.set_mesh_coors(new_coors, update_fields=True) val1 = problem.evaluate(term, var_dict=variables.as_dict()) # mesh perturbation - new_coors = modify_mesh(-dt/2., spbox, dvel_modes[ii], cp_pos) problem.set_mesh_coors(new_coors, update_fields=True) val2 = problem.evaluate(term, var_dict=variables.as_dict()) val_fd = (val1 - val2) / dt err = nm.abs(val_sa - val_fd) / nm.linalg.norm(val_sa) tst.report('term: %s' % val) tst.report('sensitivity term: %s' % val_sa) tst.report('finite differences: %s' % val_fd) tst.report('relative error: %s' % err) _ok = err < tolerance ok = ok and _ok gc.collect() assert ok