Creating multivariate Piecewise functions using pyomo.environ

[TD]

I am trying to create a piecewise linear objective from a multivariate nonlinear function. Using the pyomo kernel library, I came up with the following working example:

import pyomo.environ as pyo
import pyomo.kernel as pmo
import numpy as np
from pyomo.core.kernel.piecewise_library.transforms_nd import piecewise_nd

model = pmo.block()
model.x = pmo.variable(lb=0, ub=3)
model.y = pmo.variable(lb=0, ub=3)
model.z = pmo.variable()

def f(x, y):
    return x*y**2 + 20*np.sin(x*y)

var_list = [model.x, model.y]
tri = pmo.piecewise_util.generate_delaunay(var_list, num=20)
values = f(*tri.points.T)

model.pw = piecewise_nd(tri, values, input=var_list, output=model.z)
model.obj = pmo.objective(expr=model.z)

pyo.SolverFactory('glpk').solve(model, tee=True)

Since I would like to work with pyomo models created using the pyomo.environ components, I would prefer to use pyomo.environ.Piecewise. In the documentation, no examples using multivariate functions are given.

Question

Is it possible to recreate the above using only features from pyomo.environ?

(I also posted this question on StackOverflow here)

Thanks for any help!