def METHOD_NAME(nelems: int = 24,
         etype: str = 'triangle',
         btype: str = 'std',
         degree: int = 2,
         poisson: float = .3,
         direct: bool = False):
    '''Plane strain plate under gravitational pull
    Solves the linear elasticity problem on a unit square domain, clamped at
    the top boundary, and stretched under the influence of a vertical
    distributed load.
    Parameters
    ----------
    nelems
        Number of elements along edge.
    etype
        Type of elements (square/triangle/mixed).
    btype
        Type of basis function (std/spline), with availability depending on the
        configured element type.
    degree
        Polynomial degree.
    poisson
        Poisson's ratio, nonnegative and strictly smaller than 1/2.
    direct
        Use direct traction evaluation.
    '''
    domain, geom = mesh.unitsquare(nelems, etype)
    ns = Namespace()
    ns.δ = function.eye(domain.ndims)
    ns.x = geom
    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
    ns.add_field(('u', 't'), domain.basis(btype, degree=degree), shape=(2,))
    ns.X_i = 'x_i + u_i'
    ns.λ = 1
    ns.μ = .5/poisson - 1
    ns.ε_ij = '.5 (∇_i(u_j) + ∇_j(u_i))'
    ns.σ_ij = 'λ ε_kk δ_ij + 2 μ ε_ij'
    ns.E = 'ε_ij σ_ij'
    ns.q_i = '-δ_i1'
    sqr = domain.boundary['top'].integral('u_k u_k dS' @ ns, degree=degree*2)
    cons = solver.optimize(('u',), sqr, droptol=1e-15)
    # solve for equilibrium configuration
    internal = domain.integral('E dV' @ ns, degree=degree*2)
    external = domain.integral('u_i q_i dV' @ ns, degree=degree*2)
    args = solver.optimize('u,', internal - external, constrain=cons)
    # evaluate tractions and net force
    if direct:
        ns.t_i = 'σ_ij n_j' # <-- this is an inadmissible boundary term
    else:
        external += domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2)
        invcons = dict(t=numpy.choose(numpy.isnan(cons['u']), [numpy.nan, 0.]))
        args = solver.solve_linear(('t',), [(internal - external).derivative('u')], constrain=invcons, arguments=args)
    F = domain.boundary['top'].integrate('t_i dS' @ ns, degree=degree*2, arguments=args)
    log.user('total clamping force:', F)
    # visualize solution
    bezier = domain.sample('bezier', 3)
    X, E = bezier.eval(['X_i', 'E'] @ ns, **args)
    Xt, t = domain.boundary['top'].sample('bezier', 2).eval(['X_i', 't_i'] @ ns, **args)
    with export.mplfigure('energy.png') as fig:
        ax = fig.add_subplot(111, ylim=(-.2,1), aspect='equal')
        im = ax.tripcolor(*X.T, bezier.tri, E, shading='gouraud', rasterized=True, cmap='turbo')
        export.plotlines_(ax, X.T, bezier.hull, colors='k', linewidths=.1, alpha=.5)
        ax.quiver(*Xt.T, *t.T, clip_on=False)
        fig.colorbar(im)
    return cons, args