Complex 3D, Boundary normals for point cloud, traction condition

[ruchitkini]

Hello @lululxvi and deepxde community,

I would like to know if there is any way to calculate boundary normal for arbitrary 3D geometry implemented using PointCloud in deepxde. I am trying to do linear elasticity analysis and for traction conditions, I need the normal vector. Could anyone suggest how can that be done?

Currently, I am trying the finite difference approach, where the tangent can be calculated using the nearby points and then a normal vector from the tangent. It works well for 2D with a sufficient number of points but I am facing problems to generalise it on 3D geometry. Any help would be appreciated.

Thank you.

[lululxvi]

You need to compute using some method by yourself, and then pass into DeepXDE.

We will have a paper on arXiv in this or next week about 3D complex boundary.

[ruchitkini]

That's great. Looking forward to the paper. Thanks.

[ruchitkini]

Hello @lululxvi ,

Is the paper already published? If yes, please provide the link.

For my problem of 3D airfoil, I am creating a mesh in Gmsh and then giving the coordinates as the training points. I am also importing boundary points on the top and bottom curved surfaces of the airfoil to apply traction conditions. For the fixed end I am using hard constraints for displacement conditions and for the free end, I am applying traction conditions through Dirichlet bc as it is 0.
I am then using OperatorBC to apply traction conditions on curved surfaces, the question here is how to provide a boundary for this function as the surface is curved. For that, I called the boundary points also, through PointCloud but I am not sure if it is the correct approach. Or will PointSetOperatorBC be a better option?

In the below snippet, I am calling the points that I need from a CSV file.

geom = dde.geometry.PointCloud(coords, boundary_points=coords_boundary)
geom_top = dde.geometry.PointCloud(coords_boundary_top, boundary_points=coords_boundary_top)
geom_bottom = dde.geometry.PointCloud(coords_boundary_bottom, boundary_points=coords_boundary_bottom)

Then I am defining the boundary function as:

def boundary_top(x, on_boundary):
    return on_boundary

def boundary_bottom(x, on_boundary):
    return on_boundary

I am just returning on_boundary, and then while defining the condition I give the geometry as the set of those boundary points as shown below:

t1_top_bc = dde.icbc.OperatorBC(geom_top, traction_top1, boundary_top)
t2_top_bc = dde.icbc.OperatorBC(geom_top, traction_top2, boundary_top)
t3_top_bc = dde.icbc.OperatorBC(geom_top, traction_top3, boundary_top)

Is this the right approach? My code works, but I am not getting the correct results, I use ADAM for a few iterations and then LBFGS, the loss drops very rapidly when using ADAM and when the iteration for LBFGS starts the loss is again very high as it was during the start of the training. Even after running for more than 100,000 iterations the loss is still way higher than the loss from ADAM.
And if I use only ADAM, I get low losses but the results are completely incorrect.

Could you please provide a suggestion on this?

Thank you.

[lululxvi]

Our paper is not published yet.

[ruchitkini]

Hello @lululxvi ,

For this issue, I did some tests on a 2D plate with a simple circular curve. In the first case, I defined the normals through the torch.atan2 function to define the angles as shown:

def traction_arc11(x, f, _):
    theta = torch.atan2(x[:, 1:2], x[:, 0:1])
    n1, n2 = torch.cos(theta), torch.sin(theta)
    # n1, n2 = torch.tensor(normals[:, 0]), torch.tensor(normals[:, 1])
    s11 = (2 * mu + lmbd) * (dde.grad.jacobian(f, x, i=0, j=0)) + lmbd * (dde.grad.jacobian(f, x, i=1, j=1))
    s12 = mu * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))
    value11 = torch.mul(s11, n1) + torch.mul(s12, n2) + torch.div(n1, E)

    return value11

Here n1 and n2 are calculated inside the function using simple cos and sines. Second, I am calculating the normals for my boundary and saved in a 'CSV' file and here I am calling it and defining it through tensor for multiplication (the commented line). I ran both the cases for 50,000 iterations and found out that the loss for the second case converges very slowly compared to first case and also takes more time. After 50000, for first case: loss = 1.6 x 10e-5, time = 92.4 min, for second case: loss = 8.4 x 10e-4, time = 106.5 min. I am using ADAM optimiser. Moreover, if I keeping running the second case for like more 100,000 iteration the loss reduces to 7.2 e-6, but there are irregular results on that boundary.

Why is this happening? Is there any other way of defining normals? I have a point cloud, and I am using pyvista library for visualization and computing normals through its function. This loss and time is important for training 3D geometry.
Please suggest something on this issue.

[ruchitkini]

Hello again @lululxvi ,

What I found out is that in the function where I multiply s11*n1, my vector n1 is just for the points on the arc, but the s11 is calculated for all the points of the domain. I believe they were calculated just for the defined boundary --> which I am defining as PointSetOperatorBC and also tried with OperatorBC by defining the boundary function.

1. Why is that the stress values are calculated for all the points in the domain?
2. How do I set it for just the points I need i.e. boundary or is there any other way for this?
3. When I define a function and call it through Dirichlet bc, it gives me the required shape of x, say x[:, 0:1].shape = (9,1), but when I call the function defined above through OperatorBC or PointSetOperatorBC, the shape of x is of all coordinates say, (300,1). Why is this so, why not just the coordinates on the boundary?

Please provide a solution to this, how should I define the traction and the normal vectors? @tsarikahin @praksharma @tangqi

[tsarikahin]

@ruchitkini please check our repo based on DeepXDE (https://github.com/imcs-compsim/pinns_for_comp_mech). You will find what you need. It does not yet include 3d examples, but it should be very straightforward to implement since we use GMSH to generate geometries, which includes already boundary normals. @maxdanwitz

[DeepaMahm]

Hello @tsarikahin and @praksharma ,

I am sorry for the cross-post, I had some issues with the pointsetBC function which I had posted here #1548 (comment)

Unfortunately, the same issue persists and I am not able to fix the error
"File "<array_function internals>", line 200, in concatenate
ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 3 and the array at index 1 has size 4
"

Could you please have a look? I tried what was suggested in the link shared above i.e. changes the values to column vector, but it didn't help.

Please find below the complete code

"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle"""
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
from deepxde import config
import tensorflow as tf

# Backend tensorflow.compat.v1 or tensorflow
# from deepxde.backend import tf


def pde(x, y):
    dc_t = dde.grad.jacobian(y, x, i=0, j=1)
    dc_xx = dde.grad.hessian(y, x, i=0, j=0)

    d = 1
    return (
        dc_t
        - d * dc_xx
    )


def get_lboundary_points():
    """
    """
    lbc_pids = [452, 402, 352, 302, 251, 201, 152, 99, 13, 96]

    lbc_vertices = [
        [0., 0.9, 0.],
        [0., 0.8, 0.],
        [0., 0.7, 0.],
        [0., 0.6, 0.],
        [0., 0.5, 0.],
        [0., 0.4, 0.],
        [0., 0.3, 0.],
        [0., 0.2, 0.],
        [0., 0.1, 0.],
        [0., 0., 0.]
    ]

    return lbc_vertices


def get_rboundary_points():
    """
    """
    rbc_pids = [499, 446, 396, 346, 293, 295, 244, 193, 146, 144]

    rbc_vertices = [
        [4.9, 0.9, 0.],
        [4.9, 0.8, 0.],
        [4.9, 0.7, 0.],
        [4.9, 0.6, 0.],
        [4.9, 0.5, 0.],
        [4.9, 0.4, 0.],
        [4.9, 0.3, 0.],
        [4.9, 0.2, 0.],
        [4.9, 0.1, 0.],
        [4.9, 0., 0.]
    ]
    return rbc_vertices


def get_points():
   pts = [[0., 0., 0.],
          [0.1, 0., 0.],
          [0.2, 0., 0.],
          [0.3, 0., 0.],
          [0.4, 0., 0.],
          [0.5, 0., 0.],
          [0.6, 0., 0.],
          [0.7, 0., 0.],
          [0.8, 0., 0.],
          [0.9, 0., 0.],
          [1., 0., 0.],
          [1.1, 0., 0.],
          [1.2, 0., 0.],
          [1.3, 0., 0.],
          [1.4, 0., 0.],
          [1.5, 0., 0.],
          [1.6, 0., 0.],
          [1.7, 0., 0.],
          [1.8, 0., 0.],
          [1.9, 0., 0.],
          [2., 0., 0.],
          [2.1, 0., 0.],
          [2.2, 0., 0.],
          [2.3, 0., 0.],
          [2.4, 0., 0.],
          [2.5, 0., 0.],
          [2.6, 0., 0.],
          [2.7, 0., 0.],
          [2.8, 0., 0.],
          [2.9, 0., 0.],
          [3., 0., 0.],
          [3.1, 0., 0.],
          [3.2, 0., 0.],
          [3.3, 0., 0.],
          [3.4, 0., 0.],
          [3.5, 0., 0.],
          [3.6, 0., 0.],
          [3.7, 0., 0.],
          [3.8, 0., 0.],
          [3.9, 0., 0.],
          [4., 0., 0.],
          [4.1, 0., 0.],
          [4.2, 0., 0.],
          [4.3, 0., 0.],
          [4.4, 0., 0.],
          [4.5, 0., 0.],
          [4.6, 0., 0.],
          [4.7, 0., 0.],
          [4.8, 0., 0.],
          ]

   return pts


def model(gdata):
    geom = dde.geometry.PointCloud(gdata)
    timedomain = dde.geometry.TimeDomain(0, 60)
    geomtime = dde.geometry.GeometryXTime(geom, timedomain)

    # left_bc = dde.DirichletBC(geomtime, func_lbc, lambda _, on_boundary: on_boundary)
    # right_bc = dde.NeumannBC(geomtime, func_rbc, lambda _, on_boundary: on_boundary)
    # bc_left = dde.DirichletBC(geomtime, lambda x: 0, lambda x, _: np.isclose(x[0], 0))
    # bc_right = dde.DirichletBC(geomtime, lambda x: 1, lambda x, _: np.isclose(x[0], 1))

    lbc_vertices = get_lboundary_points()
    lbc_pts = np.array(lbc_vertices, dtype=config.real(np))
    # lbc_values = np.array([5]*len(lbc_vertices),  dtype=config.real(np))
    lbc_values = np.array([5] * len(lbc_vertices), dtype=config.real(np)).reshape(-1, 1)

    bc_l = dde.icbc.PointSetBC(lbc_pts, lbc_values, component=0)

    # rbc_vertices = get_rboundary_points()
    # rbc_pts = np.array(rbc_vertices, dtype=config.real(np))
    # rbc_values = [1]*len(rbc_vertices)
    # bc_r = dde.icbc.PointSetBC(rbc_pts, rbc_values, component=0)

    ic = dde.icbc.IC(geomtime, lambda x: 1, lambda _, on_initial: on_initial, component=0)

    data = dde.data.TimePDE(
        geomtime,
        pde,
        [bc_l, ic],
        num_domain=320,
        num_test=80000
    )

    layer_size = [2] + [30] * 6 + [1]
    activation = "tanh"
    initializer = "Glorot uniform"
    net = dde.nn.FNN(layer_size, activation, initializer)

    model = dde.Model(data, net)
    model.compile("adam", lr=1e-3, metrics=["l2 relative error"])
    losshistory, train_state = model.train(iterations=200)

    dde.saveplot(losshistory, train_state, issave=True, isplot=True)


if __name__ == '__main__':

    pts = get_points()
    gdata = np.array(pts)
    model(gdata=gdata)

Could you please have a look? I am not sure if this is a tensor flow issue.

[tsarikahin]

I think you should work on that. At first glance, the number of initial points in TimePDE is missing. Also, your geometry seems like 3d and you have time as an additional input, so I expect a network having 4 nodes at the input layer but you have only 2.

[DeepaMahm]

@tsarikahin ,
Thanks a lot for the response.
Yes, the geometry is 3D. I had specified 2 nodes since the pde defined is 1D I.e. 1 node for time and 1 node for x dimensions.
Could you please clarify if we still need to give 4 nodes at the input layer?

[ruchitkini]

@tsarikahin ,
yes, it is working now!

[lululxvi]

@tsarikahin You are welcome to submit a PR to add some useful functions in your repo into DeepXDE!

[praksharma]

Hey getting boundary normals is very easy. Just use stl geometry. The main challenge is to obtain the derivative for the Neumann or Robin BC. This probably requires fiddling with the source code. I am currently implementing the same which I am calling PointSetNeumannBC which follows PointSetBC.

[ruchitkini]

Hey @praksharma ,

Thank you for your response, I have normals on the required boundary points through a different source, the question is how do I multiply this with the stress? When I was doing a 2D simulation, for a circular arc I was using angles. So let's say I have 'm' training points and 'n' boundary points. When I did direct multiplication s11*n1, it gave an array of shapes (m X 1). I was under the assumption that since in operator bc, I am giving my boundary, it will give me stresses just for those points --> (n X 1), but that is not the case. So I need to somehow get just stresses on the boundary and then multiply them with normals.

I have a question, How does OperatorBC work? If the output I am getting through function is of shape (m X 1) and my boundary points are just (n X 1), with what values is the OperatorBC doing the comparison, why is it not throwing error?

[praksharma]

• deepxde.icbc.boundary_conditions.DirichletBC allows y(x) = func(x) type BCs.
• deepxde.icbc.boundary_conditions.OperatorBC allows func(inputs, outputs, X) = 0 type BCs.

I am not sure how proficient are you in programming, but if you look at the constructor of operatorBC, you will clearly see the usage.

 def __init__(self, geom, func, on_boundary):

Given your error you are mixing double and float. This error is not enough.

File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\torch\nn\modules\linear.py", line 114, in forward
    return F.linear(input, self.weight, self.bias)
RuntimeError: mat1 and mat2 must have the same dtype, but got Double and Float

Please provide the full error trackback. You have provided the last 3 lines of the error which is not very useful.

[ruchitkini]

I am new to coding, I had a look at the constructor but just understood from a layman's point of view, that OperatorBC gets the values defined by the function on the boundary points for minimizing the loss. I am not still able to understand the why it works if the shape of value I am getting from function --> (m X 1) is different from the boundary points (n X 1), does the OperatorBC takes the value of function only for the indices where the on_boundary is true? I think the [beg:end] is used for that. Is there any way to access this indices outside source code?

Regarding the error, below is the complete error message:

Traceback (most recent call last):
  File "D:\ThesisCodes\DeepXDE\2D_plateWithHole\training6.py", line 107, in <module>
    losshistory, train_state = model.train(epochs=2000, display_every=200)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\utils\internal.py", line 22, in wrapper
    result = f(*args, **kwargs)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\model.py", line 628, in train
    self._test()
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\model.py", line 817, in _test
    ) = self._outputs_losses(
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\model.py", line 538, in _outputs_losses
    outs = outputs_losses(inputs, targets, auxiliary_vars)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\model.py", line 315, in outputs_losses_train
    return outputs_losses(
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\model.py", line 299, in outputs_losses
    outputs_ = self.net(inputs)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\torch\nn\modules\module.py", line 1518, in _wrapped_call_impl
    return self._call_impl(*args, **kwargs)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\torch\nn\modules\module.py", line 1527, in _call_impl
    return forward_call(*args, **kwargs)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\deepxde\nn\pytorch\fnn.py", line 43, in forward
    else self.activation(linear(x))
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\torch\nn\modules\module.py", line 1518, in _wrapped_call_impl
    return self._call_impl(*args, **kwargs)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\torch\nn\modules\module.py", line 1527, in _call_impl
    return forward_call(*args, **kwargs)
  File "D:\ThesisCodes\DeepXDE\venv\lib\site-packages\torch\nn\modules\linear.py", line 114, in forward
    return F.linear(input, self.weight, self.bias)
RuntimeError: mat1 and mat2 must have the same dtype, but got Double and Float

Process finished with exit code 1

[ruchitkini]

Hello again, I am still not able to get the indices for the required boundary points so that I can fit my normals into traction calculation as got some reference from repo provided by @tsarikahin .
So I kind of did an amateur approach, got the indices manually by debugging (as shown below in code). Since all the training points are arranged in train_x and boundary points are stacked at the start followed by domain points. I am getting correct results. But the indices will keep changing as we change the coordinates. So I want a general approach for this. How can I do it? My entire code is in chats above, and below is the code snippet where I have inserted the indices manually

def traction_arc11(x, f, _):
    dummy_normal = np.zeros((len(x[:, 0]), 3))
    dummy_normal[7372:8109, 0], dummy_normal[7372:8109, 1], dummy_normal[7372:8109, 2] = (normals[:, 0], normals[:, 1],
                                                                                          normals[:, 2])
    n1, n2, n3 = (torch.as_tensor(dummy_normal[:, 0:1]), torch.as_tensor(dummy_normal[:, 1:2]),
                  torch.as_tensor(dummy_normal[:, 2:3]))
    # n1, n2, n3 = torch.tensor(normals[:, 0:1]), torch.tensor(normals[:, 1:2]), torch.tensor(normals[:, 2:3])
    s11 = ((2 * mu + lmbd) * (dde.grad.jacobian(f, x, i=0, j=0)) +
           lmbd * ((dde.grad.jacobian(f, x, i=1, j=1)) + (dde.grad.jacobian(f, x, i=2, j=2))))
    s12 = mu * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))
    s13 = mu * (dde.grad.jacobian(f, x, i=0, j=2) + dde.grad.jacobian(f, x, i=2, j=0))

    value11 = s11 * n1 + s12 * n2 + s13 * n3

    return value11


def traction_arc22(x, f, _):
    dummy_normal = np.zeros((len(x[:, 0]), 3))
    dummy_normal[8109:8846, 0], dummy_normal[8109:8846, 1], dummy_normal[8109:8846, 2] = (normals[:, 0], normals[:, 1],
                                                                                          normals[:, 2])
    n1, n2, n3 = (torch.as_tensor(dummy_normal[:, 0:1]), torch.as_tensor(dummy_normal[:, 1:2]),
                  torch.as_tensor(dummy_normal[:, 2:3]))
    # n1, n2, n3 = torch.tensor(normals[:, 0:1]), torch.tensor(normals[:, 1:2]), torch.tensor(normals[:, 2:3])
    s22 = ((2 * mu + lmbd) * (dde.grad.jacobian(f, x, i=1, j=1)) +
           lmbd * ((dde.grad.jacobian(f, x, i=0, j=0)) + (dde.grad.jacobian(f, x, i=2, j=2))))
    s23 = mu * (dde.grad.jacobian(f, x, i=1, j=2) + dde.grad.jacobian(f, x, i=2, j=1))
    s21 = mu * (dde.grad.jacobian(f, x, i=1, j=0) + dde.grad.jacobian(f, x, i=0, j=1))

    value22 = s21 * n1 + s22 * n2 + s23 * n3

    return value22


def traction_arc33(x, f, _):
    dummy_normal = np.zeros((len(x[:, 0]), 3))
    dummy_normal[8846:9583, 0], dummy_normal[8846:9583, 1], dummy_normal[8846:9583, 2] = (normals[:, 0], normals[:, 1],
                                                                                          normals[:, 2])
    n1, n2, n3 = (torch.as_tensor(dummy_normal[:, 0:1]), torch.as_tensor(dummy_normal[:, 1:2]),
                  torch.as_tensor(dummy_normal[:, 2:3]))
    # n1, n2, n3 = torch.tensor(normals[:, 0:1]), torch.tensor(normals[:, 1:2]), torch.tensor(normals[:, 2:3])
    s33 = ((2 * mu + lmbd) * (dde.grad.jacobian(f, x, i=2, j=2)) +
           lmbd * ((dde.grad.jacobian(f, x, i=1, j=1)) + (dde.grad.jacobian(f, x, i=0, j=0))))
    s32 = mu * (dde.grad.jacobian(f, x, i=2, j=1) + dde.grad.jacobian(f, x, i=1, j=2))
    s31 = mu * (dde.grad.jacobian(f, x, i=2, j=0) + dde.grad.jacobian(f, x, i=0, j=2))

    value33 = s31 * n1 + s32 * n2 + s33 * n3

    return value33

[DeepaMahm]

Hi @praksharma,

Could you please share the PointSetNeumannBC function that you have implemented, if you don't mind?

Thanks a lot,
Deepa