def METHOD_NAME():
    """Define the problem to solve."""
    from sfepy.discrete.fem.meshio import UserMeshIO
    from sfepy.mesh.mesh_generators import gen_block_mesh
    from sfepy.mechanics.matcoefs import stiffness_from_lame
    def mesh_hook(mesh, mode):
        """
        Generate the block mesh.
        """
        if mode == 'read':
            mesh = gen_block_mesh([2, 2, 3], [2, 2, 4], [0, 0, 1.5], name='el3',
                                  verbose=False)
            return mesh
        elif mode == 'write':
            pass
    filename_mesh = UserMeshIO(mesh_hook)
    options = {
        'nls' : 'newton',
        'ls' : 'ls',
        'ts' : 'ts',
        'save_times' : 'all',
    }
    functions = {
        'linear_tension' : (linear_tension,),
        'linear_compression' : (linear_compression,),
        'empty' : (lambda ts, coor, mode, region, ig: None,),
    }
    fields = {
        'displacement' : ('real', 3, 'Omega', 1),
    }
    # Coefficients are chosen so that the tangent stiffness is the same for all
    # material for zero strains.
    # Young modulus = 10 kPa, Poisson's ratio = 0.3
    materials = {
        'solid' : ({
            'K'  : 8.333, # bulk modulus
            'mu_nh' : 3.846, # shear modulus of neoHookean term
            'mu_mr' : 1.923, # shear modulus of Mooney-Rivlin term
            'kappa' : 1.923, # second modulus of Mooney-Rivlin term
            # elasticity for LE term
            'D' : stiffness_from_lame(dim=3, lam=5.769, mu=3.846),
        },),
        'load' : 'empty',
    }
    variables = {
        'u' : ('unknown field', 'displacement', 0),
        'v' : ('test field', 'displacement', 'u'),
    }
    regions = {
        'Omega' : 'all',
        'Bottom' : ('vertices in (z < 0.1)', 'facet'),
        'Top' : ('vertices in (z > 2.9)', 'facet'),
    }
    ebcs = {
        'fixb' : ('Bottom', {'u.all' : 0.0}),
        'fixt' : ('Top', {'u.[0,1]' : 0.0}),
    }
    integrals = {
        'i' : 1,
        'isurf' : 2,
    }
    equations = {
        'linear' : """dw_lin_elastic.i.Omega(solid.D, v, u)
                    = dw_surface_ltr.isurf.Top(load.val, v)""",
        'neo-Hookean' : """dw_tl_he_neohook.i.Omega(solid.mu_nh, v, u)
                         + dw_tl_bulk_penalty.i.Omega(solid.K, v, u)
                         = dw_surface_ltr.isurf.Top(load.val, v)""",
        'Mooney-Rivlin' : """dw_tl_he_neohook.i.Omega(solid.mu_mr, v, u)
                           + dw_tl_he_mooney_rivlin.i.Omega(solid.kappa, v, u)
                           + dw_tl_bulk_penalty.i.Omega(solid.K, v, u)
                           = dw_surface_ltr.isurf.Top(load.val, v)""",
    }
    solvers = {
        'ls' : ('ls.scipy_direct', {}),
        'newton' : ('nls.newton', {
            'i_max'      : 5,
            'eps_a'      : 1e-10,
            'eps_r'      : 1.0,
        }),
        'ts' : ('ts.simple', {
            't0'    : 0,
            't1'    : 1,
            'dt'    : None,
            'n_step' : 101, # has precedence over dt!
            'verbose' : 1,
        }),
    }
    return locals()