Hard boundary for analytic but nontrivial boundaries

[akaptano]

Hi all,

My impression is that complex geometries are usually not implemented with a hard boundary condition. However, I am not sure why this should be impossible, so I am looking at implementing hard boundary constraints for a 2D PDE solve on a toroidal cross section, which can be parametrized by an angle $\theta$ and a set of constants $\epsilon$, $\delta$, $\kappa$, such that:
$x = 1 + \epsilon cos(\theta + arcsin(\delta) sin(\theta))$
$y = \epsilon \kappa sin(\theta)$

The idea now is to use the normal apply_output_transform() function and simply use a function F(x) such that F(x) = 0 on the cross sectional boundary and F(x) = 1 otherwise.
The relevant code is:

spatial_domain = dde.geometry.Ellipse(eps, kappa, delta)   # user-defined geometry object corresponding to the toroidal parametrization
def modify_output(x, y):
    return spatial_domain.strictly_inside(x) * y

where the strictly inside function is defined through:

def strictly_inside(self, x):
        bool_x = tf.map_fn(self.is_point_in_poly, x, fn_output_signature=tf.float32)  # check if each point is in the ellipse
        bool_x = tf.expand_dims(bool_x, 1)  # match shape of output vector y
        return bool_x

def is_point_in_poly(self, x_full):
        x = x_full[0:1]
        y = x_full[1:2]
        poly = self.x_ellipse
        num = len(poly)
        j = num - 1
        c = False
        for i in range(num):
            if (x == poly[i][0]) and (y == poly[i][1]):
                # point is a corner
                return 0.0
            if ((poly[i][1] > y) != (poly[j][1] > y)):
                slope = (x-poly[i][0])*(poly[j][1]-poly[i][1])-(poly[j][0]-poly[i][0])*(y-poly[i][1])
                if slope == 0:
                    # point is on boundary
                    return 0.0
                if (slope < 0) != (poly[j][1] < poly[i][1]):
                    c = not c
            j = i
        return float(c)

Here self.x_ellipse is the parametrized cross section shape.
The code seems to run without error:

bc_dirichlet = dde.DirichletBC(spatial_domain, lambda x: 0, lambda _, on_boundary: on_boundary)

data = dde.data.PDE(
    spatial_domain,
    pde,
    bc_dirichlet,
    num_domain=256,
    num_boundary=2000,
    num_test=n_test,
    train_distribution="LHS"
)
net = dde.maps.FNN(
    [2] + 3 * [20] + [1], 
    "swish", "Glorot normal")

net.apply_output_transform(modify_output)  # apply the hard boundary

model = dde.Model(data, net)
model.compile(
      "adam", 
      loss_weights=[1,1], 
)
loss_history, train_state = model.train(epochs=100, display_every=1)  # trains without error, but boundary loss is nonzero

but the boundary loss term appears nonzero, so I appear to be missing something. I would provide a full working example but there is significant code behind the scenes for the new geometry.

[akaptano]

For anyone who finds this later -- looks like the implementation above doesn't make a lot of sense since tensorflow needs to be able to back-propagate the gradients of this thing. Also, the ellipse cross section is invertible anyways (it just looks complicated) and I was able to implement the hard boundary conditions.