Source code for sfepy.tests.test_laplace_unit_square

import pytest

from sfepy import data_dir
import sfepy.base.testing as tst

filename_mesh = data_dir + '/meshes/2d/square_unit_tri.mesh'

material_1 = {
    'name' : 'coef',
    'values' : {
        'val' : 1.0,
    },
}
material_2 = {
    'name' : 'm',
    'values' : {
        'K' : [[1.0, 0.0], [0.0, 1.0]],
    },
}

field_1 = {
    'name' : 'a_harmonic_field',
    'dtype' : 'real',
    'shape' : 'scalar',
    'region' : 'Omega',
    'approx_order' : 2,
}

variable_1 = {
    'name' : 't',
    'kind' : 'unknown field',
    'field' : 'a_harmonic_field',
    'order' : 0,
}
variable_2 = {
    'name' : 's',
    'kind' : 'test field',
    'field' : 'a_harmonic_field',
    'dual' : 't',
}

region_1000 = {
    'name' : 'Omega',
    'select' : 'all',
}

region_1 = {
    'name' : 'Left',
    'select' : 'vertices in (x < -0.499)',
    'kind' : 'facet',
}
region_2 = {
    'name' : 'Right',
    'select' : 'vertices in (x > 0.499)',
    'kind' : 'facet',
}
region_3 = {
    'name' : 'Gamma',
    'select' : 'vertices of surface',
    'kind' : 'facet',
}

ebc_1 = {
    'name' : 't_left',
    'region' : 'Left',
    'dofs' : {'t.0' : 5.0},
}
ebc_2 = {
    'name' : 't_right',
    'region' : 'Right',
    'dofs' : {'t.0' : 0.0},
}

integral_1 = {
    'name' : 'i',
    'order' : 2,
}

equations = {
    'Temperature' : """dw_laplace.i.Omega(coef.val, s, t) = 0"""
}

solution = {
    't' : '- 5.0 * (x - 0.5)',
}

solver_0 = {
    'name' : 'ls',
    'kind' : 'ls.scipy_direct',
}

solver_1 = {
    'name' : 'newton',
    'kind' : 'nls.newton',

    'i_max'      : 1,
    'eps_a'      : 1e-10,
}

lin_min, lin_max = 0.0, 2.0

[docs] def linear(bc, ts, coor, which): vals = coor[:,which] min_val, max_val = vals.min(), vals.max() vals = (vals - min_val) / (max_val - min_val) * (lin_max - lin_min) + lin_min return vals
[docs] def linear_x(bc, ts, coor): return linear(bc, ts, coor, 0)
[docs] def linear_y(bc, ts, coor): return linear(bc, ts, coor, 1)
[docs] def linear_z(bc, ts, coor): return linear(bc, ts, coor, 2)
[docs] @pytest.fixture(scope='module') def data(): import sys from sfepy.applications import solve_pde from sfepy.base.conf import ProblemConf from sfepy.base.base import Struct conf = ProblemConf.from_dict(globals(), sys.modules[__name__]) problem, state = solve_pde(conf, save_results=False) return Struct(problem=problem, state=state)
[docs] def test_solution(data): problem = data.problem sol = problem.conf.solution vec = data.state() variables = problem.get_variables() ok = True for var_name, expression in sol.items(): coor = variables[var_name].field.get_coor() ana_sol = tst.eval_coor_expression(expression, coor) num_sol = variables.get_vec_part(vec, var_name) ret = tst.compare_vectors(ana_sol, num_sol, label1='analytical %s' % var_name, label2='numerical %s' % var_name) if not ret: tst.report('variable %s: failed' % var_name) ok = ok and ret assert ok
[docs] def test_boundary_fluxes(data, output_dir): import os.path as op from sfepy.linalg import rotation_matrix2d from sfepy.discrete import Material problem = data.problem angles = [0, 30, 45] region_names = ['Left', 'Right', 'Gamma'] values = [5.0, -5.0, 0.0] variables = problem.get_variables() get_state = variables.get_vec_part state = data.state.copy() problem.time_update(ebcs={}, epbcs={}) state.apply_ebc() nls = problem.get_nls() aux = nls.fun(state()) field = variables['t'].field conf_m = problem.conf.get_item_by_name('materials', 'm') m = Material.from_conf(conf_m, problem.functions) name = op.join(output_dir, op.split(problem.domain.mesh.name)[1] + '_%02d.mesh') orig_coors = problem.get_mesh_coors().copy() ok = True for ia, angle in enumerate(angles): tst.report('%d: mesh rotation %d degrees' % (ia, angle)) problem.domain.mesh.transform_coors(rotation_matrix2d(angle), ref_coors = orig_coors) problem.set_mesh_coors(problem.domain.mesh.coors, update_fields=True) problem.domain.mesh.write(name % angle, io='auto') for ii, region_name in enumerate(region_names): flux_term = 'ev_surface_flux.i.%s(m.K, t)' % region_name val1 = problem.evaluate(flux_term, t=variables['t'], m=m) rvec = get_state(aux, 't', True) reg = problem.domain.regions[region_name] nods = field.get_dofs_in_region(reg, merge=True) val2 = rvec[nods].sum() # Assume 1 dof per node. ok = ok and ((abs(val1 - values[ii]) < 1e-10) and (abs(val2 - values[ii]) < 1e-10)) tst.report(' %d. %s: %e == %e == %e' % (ii, region_name, val1, val2, values[ii])) # Restore original coordinates. problem.domain.mesh.transform_coors(rotation_matrix2d(0), ref_coors=orig_coors) problem.set_mesh_coors(problem.domain.mesh.coors, update_fields=True) assert ok