Order of boundary points in on_boundary - Does it match the order in boundary_points argument of PointCloud?

[DeepaMahm]

Hi Prof. @lululxvi ,

This is a follow-up to my previous question https://github.com/lululxvi/deepxde/issues/1679.

I've a 2D geometry created using PointCloud.

geom = dde.geometry.PointCloud(
     points=np.asarray(domain_points),
     boundary_points=np.asarray(boundary_points),
 )

 timedomain = dde.geometry.TimeDomain(0, 60)
 geomtime = dde.geometry.GeometryXTime(geom, timedomain)

The geometry is shown below

Question:

For instance, when I want to specify left boundary condition for the 2D geometry (marked blue in the above figure)

Referring to https://github.com/lululxvi/deepxde/issues/22#issuecomment-605493081,

I'm doing the following for this simple geometry

def boundary_left(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)  # x coordinate x[0] is 0

def boundary_right(x, on_boundary):
    return on_boundary and np.isclose(x[0], 20)  # x coordinate x[0] is 20

I would like to know what is the order of points in the bool variable on_boundary. Does it match the order of boundary_points in PointCloud?

I am asking this because I know the indices of left and right boundary points (colored blue) in boundary_points (which also includes top and bottom boundary points).

So ìn the boundary_left function, can I directly return the bool and avoid doing np.isclose(x[0], 0)?

i.e. instead of
return on_boundary and np.isclose(x[0], 0) # x coordinate x[0] is 0 can I return?

return np.ones(len(lbc_points), dtype=bool) + np.zeros(len(other_boundary_points), dtype=bool)

Could you please clarify?

@praksharma Could you please have a look at this post when you have time? Thanks a lot

[praksharma]

on_boundary will only work when you define the geometry using modules available in dde.geometry . Since you know the coordinates of the boundaries and from the plot i can see they are exactly, at x=0,20, you can directly specify the nodal coordinates.

def boundary_left(x, on_boundary):
    return x[0]== 0

def boundary_right(x, on_boundary):
    return x[0]==20

Coming to your second query. I could not understand the meaning of "order of points". Are you talking about the magnitude of the coodinates at those points?

No you can't do this.

  return np.ones(len(lbc_points), dtype=bool) + np.zeros(len(other_boundary_points), dtype=bool)

The boundaries of a particular BCs is defined using lambda function. Here is an example:

bc = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)

You are returning the same array each time the lambda function is called for a particular value of x. This will give you very unexpected results (if DeepXDE doesn't throw an error) as you are trying to assign same BC to the same input each time irrespective of your x.

[DeepaMahm]

@praksharma Thanks a lot for the clarification and thanks for always guiding me towards the right direction.

I was actually trying to figure out if I can still use the dde.icbc.NeumannBC() for geometries that are not created using modules available in dde.geometry (i.e. for geometries created using pointcould).

From your response, I understand that it wouldn't be possible since on_boundary works only when geom is created using dde.geometry. I had a look at the NeumannBC class to try and implement PointSetBCNeumann for pointcloud geometries. Unfortunately, I am not sure how to proceed. Do we need the boundary_normals (which is an input argument in ) to create a PointSetBCNeumann function ?

n.dy/dn = 0 (https://github.com/lululxvi/deepxde/issues/22#issuecomment-607704951)

class NeumannBC(BC):
    """Neumann boundary conditions: dy/dn(x) = func(x)."""

    def __init__(self, geom, func, on_boundary, component=0):
        super().__init__(geom, on_boundary, component)
        self.func = npfunc_range_autocache(utils.return_tensor(func))

    def error(self, X, inputs, outputs, beg, end, aux_var=None):
        values = self.func(X, beg, end, aux_var)
        return self.normal_derivative(X, inputs, outputs, beg, end) - values

I have geometries which include sections like the below (and Neumann needs to be specified on the edge points)

and I am not sure if boundary_normals are required.

If your implementation of PointSetBCNeumann is available elsewhere, could you please direct me to it?

[praksharma]

You can only use Neumann BC when you provide the normal vectors. In DeepXDE, the base geometry class has a function named boundary_normal, which stores the normal vectors.

https://github.com/lululxvi/deepxde/blob/44f9324236dc45f865af73430b46d0dd5b2b9ebc/deepxde/geometry/geometry.py#L77

The point cloud class allows you to manually input the normal vectors for each boundary, as it inherits the Geometry base class.

https://github.com/lululxvi/deepxde/blob/44f9324236dc45f865af73430b46d0dd5b2b9ebc/deepxde/geometry/pointcloud.py#L20

So, here when you define your pointcloud. You can easily add normal vectors for each boundary_points you have specified.

geom = dde.geometry.PointCloud(
     points=np.asarray(domain_points),
     boundary_points=np.asarray(boundary_points),
 )

You can use if-else condition to assign normal vectors. For your left boundary, i.e. x=0, all entries in boundary_normals must be -1 and for right boundary it must be +1. This is because boundary_normals indicates normally outward vector to the boundary. I hope this is the correct convention, I haven't used FEM for a long time.

Once you have created a vector of boundary_normals using if-else condition. You can add them in your geom and use dde.icbc.NeumannBC just like you'd use with a circle or a rectangle.

geom = dde.geometry.PointCloud(
     points=np.asarray(domain_points),
     boundary_points=np.asarray(boundary_points),
     boundary_normals=<your array here>
 )

I never implemented the PointSetBCNeumann because calculating normal vectors requires external libraries that implement k-nearest neighbors.

In my case, I had a MED file of my geometry, which automatically stores the boundary normals vectors. Also, I am not using DeepXDE because I am doing inverse modelling.

[DeepaMahm]

Hi @praksharma ,

Thanks a ton!

I never implemented the PointSetBCNeumann because calculating normal vectors requires external libraries that implement k-nearest neighbors. And sorry, I got a bit confused reading this comment https://github.com/lululxvi/deepxde/discussions/1633#discussioncomment-8721874.

I managed to get the boundary normals of my geometry. For a simple 2D geometry shared in the initial post, it looks like the below .

I've specified the boundary_normals as inputs in the PointCloud function, and it works nicely. Thank you.

geom = dde.geometry.PointCloud( points=np.asarray(domain_points), boundary_points=boundary_points, boundary_normals=boundary_normals )

Now, for a regular geometry e.g. rectangle created using pointcloud, we could specify NeumannBC as follows

def func_neumann_bc(x):
    return np.zeros((len(x), 1), dtype=np.float64)

def boundary_right(x, on_boundary):
    """
    NC boundary for model_v1
    :param x:
    :param on_boundary:
    :return:
    """
    return on_boundary and np.isclose(x[0], 20)  # x[0]==20

   right_bc = dde.icbc.NeumannBC(geomtime, func_neumann_bc, boundary_right)

Question

For a point cloud geometry like the below , I am not sure how to specify the NeumannBC for the face indicated by the red arrow since x[0] is not a constant value due to the geometry tilt. This is the reason I was trying to figure out how to write a function for PointSetBCNeumann. .

Also, I am not using DeepXDE because I am doing inverse modeling.

Could you please let me know if there are other libraries for inverse modeling? I'm still a beginner in this field but I need to do inverse modeling. If you could give some suggestions, it would be really helpful.

[praksharma]

Glad it worked. I have solved similar problems, where the boundary is complicated.
Since your boundaries are straight lines, you could try using slope-intercept form of line: y=mx+c. You need to manually determine the slope (m) and intercept (c) for each boundary since they are straight line. Slope (m) is fairly easy: y/x and intercept (c) is the point where the line crosses the y-axis. You need to calculate them manually.

Once, you determine m and c. You could try this.

def boundary_left(x, on_boundary):
    return on_boundary and np.isclose(x[1], x[0] * m + c) # x[1] is y and x[0] is x

From the figure, I can say the slope is same for both left and right boundary, but the intercept will change.

DeepXDE solves inverse problem involving parameter estimation. You look in the DeepXDE's docs to learn inverse modelling. My research in not parameter estimation, thus I am using my own code. I am working on solution reconstruction without the knowldege of BCs.

I just noticed, the boundary normal in your figure is mostly correct. The normals are undefined on the corners, so just exclude those points. Mathematically, the derivative is undefined on corners.

[DeepaMahm]

Hi @praksharma ,

Really sorry for the late reply.

Thanks a lot for the detailed explanation. I'm excluding the corner points too, thanks. Happy to know about the interesting research area you are working on.

If I may ask another doubt, The actual geometry is a bit more complicated.

I now have the boundary normals of the boundary vertices for the geometry. But I am not sure how to specify the Neumann BC.

I have the coordinates of points for which Neumann has to be applied saved in an array. But not sure how to specify it

Here, can we specify single points and use multiple and or use for loop in this function below?

def boundary_neumann(x, on_boundary):
    return on_boundary and np.isclose(x[1], x[0] * m + c) # x[1] is y and x[0] is x

For instance, should I split the lines forming the boundary and compute slope and intercept for each line and set up Neumann as return on_boundary and np.isclose(x[1], x[0] * m1 + c1) and np.isclose(x[1], x[0] * m2 + c2) and ... # x[1] is y and x[0] is x

Thanks a lot for your kind attention.

[praksharma]

For such complex geometry, you can't use filter function which is a combination of on_boundary and np.isclose. You could try writing a code from scratch. Here is a code I wrote for my work, 2 years ago. You can easily modify it to include Neumann BC. There is no need to calculate the normal vectors. You can just supply coordinates and corresponding BC. The code is bit long and might take time to understand.

https://github.com/praksharma/Stiff-PDEs-and-Physics-Informed-Neural-Networks/blob/main/1.%20Problem%201/1.%20Model%201/Improved_structure.ipynb

[DeepaMahm]

Thanks a ton for sharing your code. This is of great value. I will share the updates here soon