Linear Elasticity

Introduction

This tutorial illustrates the linear elasticity implementation in Modulus. Modulus offers the capability of solving the linear elasticity equations in the differential or variational form, allowing to solve a wide range of problems with a variety of boundary conditions. Three examples are presented in this chapter, namely the 3D bracket, the fuselage panel, and the plane displacement, to discuss the details of the linear elasticity in Modulus. In this tutorial, you will learn:

• How to solve linear elasticity equations using the differential and variational forms.

• How to solve linear elasticity problems for 3d and thin 2d structures (plane-stress).

• How to non-dimensionalize the elasticity equations.

Note

This tutorial assumes that you have completed Lid Driven Cavity Background tutorial and have familiarized yourself with the basics of the Modulus APIs. See Weak solution of PDEs using PINNs for more details on weak solutions of PDEs. Some of the boundary conditions are defined more elaborately in the tutorial Interface Problem by Variational Method .

The linear elasticity equations in Modulus can be found in the source code directory modulus/PDES/linear_elasticity.py.

Linear Elasticity in the Differential Form

Linear elasticity equations in the displacement form

The displacement form of the (non-transient) linear elasticity equations, known as the Navier equations, is defined as

(139)\[(\lambda + \mu) u_{j,ji} + \mu u_{i,jj} + f_i = 0,\]

where \(u_i\) is the displacement vector, \(f_i\) is the body force per unit volume, and \(\lambda, \mu\) are the Lamé parameters defined as

(140)\[\lambda = \frac{E \nu}{(1+\nu)(1-2\nu)},\]

(141)\[\mu = \frac{E}{2(1+\nu)}.\]

Here, \(E\), \(v\) are, respectively, the Young’s modulus and Poisson’s ratio.

Linear elasticity equations in the mixed form

In addition to the displacement formulation, linear elasticity can also be described by the mixed-variable formulation. In fact, based on our experiments and also the studies reported in 1, the mixed-variable formulation is easier for a neural network solver to learn, possibly due to the lower order differential terms. In the mixed-variable formulation, the equilibrium equations are defend as:

(142)\[ \sigma_{ji,j} + f_i = 0,\]

where \(\sigma_{ij}\) is the Cauchy stress tensor. the stress-displacement equations are also defined as

(143)\[\sigma_{ij} = \lambda \epsilon_{kk} \delta_{ij} + 2 \mu \epsilon_{ij},\]

where \(\delta_{ij}\) is the Kronecker delta function and \(\epsilon_{ij}\) is the strain tensor that takes the following form

(144)\[\epsilon_{ij} = \frac{1}{2}\left( u_{i,j} + u_{j,i} \right).\]

Non-dimensionalized linear elasticity equations

It is often advantageous to work with the non-dimensionalized and normalized variation of the elasticity equations to improve the training convergence and accuracy. Therefore, we define non-dimensionalized variables as:

(145)\[\hat{x}_i = \frac{x_i}{L},\]

(146)\[\hat{u}_i = \frac{u_i}{U},\]

(147)\[\hat{\lambda} = \frac{\lambda}{\mu_c},\]

(148)\[\hat{\mu} = \frac{\mu}{\mu_c}\]

Here, \(L\) is the characteristic length, \(U\) is the characteristic displacement, and \(\mu_c\) is the non-dimensionalizing shear modulus. We can obtain the non-dimensionalized Navier and equilibrium equations by multiplying both sides of equations by \(L^2/\mu_c U\), as follows:

(149)\[(\hat{\lambda} + \hat{\mu}) \hat{u}_{j,ji} + \hat{\mu} \hat{u}_{i,jj} + \hat{f}_i = 0,\]

(150)\[\hat{\sigma}_{ji,j} + \hat{f_i} = 0,\]

where the non-dimensionalized body force and stress tensor are

(151)\[\hat{f}_{i} = \frac{L^2}{\mu_c U} f_{i},\]

(152)\[\hat{\sigma}_{ij} = \frac{L}{\mu_c U}\sigma_{ij}.\]

Similarly, the non-dimensionalized form of the stress-displacement equations are obtained by multiplying both sides of equations by \(L/\mu_c U\), as follows:

(153)\[ \hat{\sigma}_{ij} = \hat{\lambda} \hat{\epsilon}_{kk} \delta_{ij} + 2 \hat{\mu} \hat{\epsilon}_{ij},\]

(154)\[\hat{\epsilon}_{ij} = \frac{1}{2}\left( \hat{u}_{i,j} + \hat{u}_{j,i} \right).\]

Plane stress equations

In a plane stress setting for thin structures, we assume that

(155)\[\hat{\sigma}_{zz} = \hat{\sigma}_{xz} = \hat{\sigma}_{yz} = 0,\]

and therefore, the following relationship holds

(156)\[\hat{\sigma}_{zz} = \hat{\lambda} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} + \frac{\partial \hat{w}}{\partial \hat{z}} \right) + 2 \hat{\mu} \frac{\partial \hat{w}}{\partial \hat{z}} = 0 \Rightarrow \frac{\partial \hat{w}}{\partial \hat{z}} = \frac{-\hat{\lambda}}{(\hat{\lambda} +2\hat{\mu})} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} \right).\]

Accordingly, the equations for \(\hat{\sigma}_{xx}\) and \(\hat{\sigma}_{yy}\) can be updated as follows

(157)\[ \hat{\sigma}_{xx} = \hat{\lambda} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} + \frac{-\hat{\lambda}}{(\hat{\lambda} +2\hat{\mu})} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} \right) \right) + 2 \hat{\mu} \frac{\partial \hat{u}}{\partial \hat{x}}\]

(158)\[ \hat{\sigma}_{yy} = \hat{\lambda} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} + \frac{-\hat{\lambda}}{(\hat{\lambda} +2\hat{\mu})} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} \right) \right) + 2 \hat{\mu} \frac{\partial \hat{v}}{\partial \hat{y}}.\]

Problem 1: Deflection of a bracket

This linear elasticity example, shows a 3D bracket in Fig. 51. This example is partially adopted from the MATLAB PDE toolbox. The back face of this bracket is clamped, and a traction of \(4 \times 10^4 \, \text{Pa}\) is applied to the front face in the negative-\(z\) direction (this face is shown in red). The rest of the surface is considered as traction-free boundaries. Set \((E, \nu) = (100 \,\text{GPa}, 0.3)\). You can non-dimensionalize the linear elasticity equations by setting \(L=1 \,\text{m}, U=0.0001 \,\text{m}, \mu_c=0.01 \mu\).

Fig. 51 Geometry of the bracket. The back face of the bracket is clamped, and a shear stress is applied to the front face in the negatize-\(z\) direction.

    # Specify parameters
    nu = 0.3
    E = 100e9
    lambda_ = nu * E / ((1 + nu) * (1 - 2 * nu))
    mu = E / (2 * (1 + nu))
    mu_c = 0.01 * mu
    lambda_ = lambda_ / mu_c
    mu = mu / mu_c
    characteristic_length = 1.0
    characteristic_displacement = 1e-4
    sigma_normalization = characteristic_length / (characteristic_displacement * mu_c)
    T = -4e4 * sigma_normalization

In general, the characteristic length can be chosen in such a way to bound the largest dimension of the geometry to \((-1,1)\). The characteristic displacement and \(\mu_c\) can also be chosen such that the maximum displacement and the applied traction are close to 1 in order. Although the maximum displacement is not known a priori, we have observed that the convergence is not sensitive to the choice of the characteristic displacement as long as it is reasonably selected based on an initial guess for the approximate order of displacement. More information on non-dimensionalizing the PDEs can be found in Scaling of Differential Equations.

Case Setup and Results

The complete python script for this problem can be found at examples/bracket/bracket.py. Two separate neural networks for displacement and stresses are used as follows

    # make list of nodes to unroll graph on
    le = LinearElasticity(lambda_=lambda_, mu=mu, dim=3)
    disp_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("z")],
        output_keys=[Key("u"), Key("v"), Key("w")],
        cfg=cfg.arch.fully_connected,
    )
    stress_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("z")],
        output_keys=[
            Key("sigma_xx"),
            Key("sigma_yy"),
            Key("sigma_zz"),
            Key("sigma_xy"),
            Key("sigma_xz"),
            Key("sigma_yz"),
        ],
        cfg=cfg.arch.fully_connected,
    )
    nodes = (
        le.make_nodes()
        + [disp_net.make_node(name="displacement_network", jit=cfg.jit)]
        + [stress_net.make_node(name="stress_network", jit=cfg.jit)]
    )

The mixed-form of the linear elasticity equations is used here in this example, and therefore, the training constraints are defined as shown below:

    # add constraints to solver
    # make geometry
    x, y, z = Symbol("x"), Symbol("y"), Symbol("z")
    support_origin = (-1, -1, -1)
    support_dim = (0.25, 2, 2)
    bracket_origin = (-0.75, -1, -0.1)
    bracket_dim = (1.75, 2, 0.2)
    cylinder_radius = 0.1
    cylinder_height = 2.0
    aux_lower_origin = (-0.75, -1, -0.1 - cylinder_radius)
    aux_lower_dim = (cylinder_radius, 2, cylinder_radius)
    aux_upper_origin = (-0.75, -1, 0.1)
    aux_upper_dim = (cylinder_radius, 2, cylinder_radius)
    cylinder_lower_center = (-0.75 + cylinder_radius, 0, -0.1 - cylinder_radius)
    cylinder_upper_center = (-0.75 + cylinder_radius, 0, 0.1 + cylinder_radius)
    cylinder_hole_radius = 0.7
    cylinder_hole_height = 0.5
    cylinder_hole_center = (0.125, 0, 0)

    support = Box(
        support_origin,
        (
            support_origin[0] + support_dim[0],
            support_origin[1] + support_dim[1],
            support_origin[2] + support_dim[2],
        ),
    )
    bracket = Box(
        bracket_origin,
        (
            bracket_origin[0] + bracket_dim[0],
            bracket_origin[1] + bracket_dim[1],
            bracket_origin[2] + bracket_dim[2],
        ),
    )
    aux_lower = Box(
        aux_lower_origin,
        (
            aux_lower_origin[0] + aux_lower_dim[0],
            aux_lower_origin[1] + aux_lower_dim[1],
            aux_lower_origin[2] + aux_lower_dim[2],
        ),
    )
    aux_upper = Box(
        aux_upper_origin,
        (
            aux_upper_origin[0] + aux_upper_dim[0],
            aux_upper_origin[1] + aux_upper_dim[1],
            aux_upper_origin[2] + aux_upper_dim[2],
        ),
    )
    cylinder_lower = Cylinder(cylinder_lower_center, cylinder_radius, cylinder_height)
    cylinder_upper = Cylinder(cylinder_upper_center, cylinder_radius, cylinder_height)
    cylinder_hole = Cylinder(
        cylinder_hole_center, cylinder_hole_radius, cylinder_hole_height
    )
    cylinder_lower.rotate(np.pi / 2, "x", cylinder_lower_center)
    cylinder_upper.rotate(np.pi / 2, "x", cylinder_upper_center)

    curve_lower = aux_lower - cylinder_lower
    curve_upper = aux_upper - cylinder_upper
    geo = support + bracket + curve_lower + curve_upper - cylinder_hole

    # Doamin bounds
    bounds_x = (-1, 1)
    bounds_y = (-1, 1)
    bounds_z = (-1, 1)
    bounds_support_x = (-1, -0.65)
    bounds_support_y = (-1, 1)
    bounds_support_z = (-1, 1)
    bounds_bracket_x = (-0.65, 1)
    bounds_bracket_y = (-1, 1)
    bounds_bracket_z = (-0.1, 0.1)

    # make domain
    domain = Domain()

    # back BC
    backBC = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0, "v": 0, "w": 0},
        batch_size=cfg.batch_size.backBC,
        lambda_weighting={"u": 10, "v": 10, "w": 10},
        criteria=Eq(x, support_origin[0]),
    )
    domain.add_constraint(backBC, "backBC")

    # front BC
    frontBC = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"traction_x": 0, "traction_y": 0, "traction_z": T},
        batch_size=cfg.batch_size.frontBC,
        criteria=Eq(x, bracket_origin[0] + bracket_dim[0]),
    )
    domain.add_constraint(frontBC, "frontBC")

    # surface BC
    surfaceBC = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"traction_x": 0, "traction_y": 0, "traction_z": 0},
        batch_size=cfg.batch_size.surfaceBC,
        criteria=(x > support_origin[0]) & (x < bracket_origin[0] + bracket_dim[0]),
    )
    domain.add_constraint(surfaceBC, "surfaceBC")

    # support interior
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={
            "equilibrium_x": 0.0,
            "equilibrium_y": 0.0,
            "equilibrium_z": 0.0,
            "stress_disp_xx": 0.0,
            "stress_disp_yy": 0.0,
            "stress_disp_zz": 0.0,
            "stress_disp_xy": 0.0,
            "stress_disp_xz": 0.0,
            "stress_disp_yz": 0.0,
        },
        batch_size=cfg.batch_size.interior_support,
        bounds={x: bounds_support_x, y: bounds_support_y, z: bounds_support_z},
        lambda_weighting={
            "equilibrium_x": geo.sdf,
            "equilibrium_y": geo.sdf,
            "equilibrium_z": geo.sdf,
            "stress_disp_xx": geo.sdf,
            "stress_disp_yy": geo.sdf,
            "stress_disp_zz": geo.sdf,
            "stress_disp_xy": geo.sdf,
            "stress_disp_xz": geo.sdf,
            "stress_disp_yz": geo.sdf,
        },
    )
    domain.add_constraint(interior, "interior_support")

    # bracket interior
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={
            "equilibrium_x": 0.0,
            "equilibrium_y": 0.0,
            "equilibrium_z": 0.0,
            "stress_disp_xx": 0.0,
            "stress_disp_yy": 0.0,
            "stress_disp_zz": 0.0,
            "stress_disp_xy": 0.0,
            "stress_disp_xz": 0.0,
            "stress_disp_yz": 0.0,
        },
        batch_size=cfg.batch_size.interior_bracket,
        bounds={x: bounds_bracket_x, y: bounds_bracket_y, z: bounds_bracket_z},
        lambda_weighting={
            "equilibrium_x": geo.sdf,
            "equilibrium_y": geo.sdf,
            "equilibrium_z": geo.sdf,
            "stress_disp_xx": geo.sdf,
            "stress_disp_yy": geo.sdf,
            "stress_disp_zz": geo.sdf,
            "stress_disp_xy": geo.sdf,
            "stress_disp_xz": geo.sdf,
            "stress_disp_yz": geo.sdf,
        },
    )
    domain.add_constraint(interior, "interior_bracket")

The training constraints consists of two different sets of interior points (i.e., interior_support and interior_bracket). This is done only to generate the interior points more efficiently.

Fig. 52 shows the Modulus results and also a comparison with a commercial solver results. The results of these two solvers show good agreement, with only a 8% difference in maximum bracket displacement.

Fig. 52 Modulus linear elasticity results for the bracket deflection example

Problem 2: Stress analysis for aircraft fuselage panel

This example uses Modulus for the analysis of stress concentration in an aircraft fuselage panel. Depending on the altitude of the flying plane, the fuselage panel is exposed to different values of hoop stress that can cause accumulated damage to the panel over the time. Therefore, in cumulative damage modeling of aircraft fuselage for the purpose of design and maintenance of aircrafts, it is required to perform several stress simulations for different hoop stress values. Consider a simplified aircraft fuselage panel as shown in Fig. 53. This example constructs a parameterized model that, once trained, can predict the stress and displacement in the panel for different values of hoop stress.

Fig. 53 Geometry and boundary conditions of the simplified aircraft fuselage panel.

Case Setup and Results

The panel material is Aluminium 2024-T3, with \((E, \nu) = (73 \,\text{GPa}, 0.33)\). We train a parameterized model with varying \(\sigma_{hoop} \in (46, 56.5)\). Since the thickness of the panel is very small (i.e., \(2 \, \text{mm}\)), solve the plane stress equations in this example. The python script for this problem can be found at examples/fuselage_panel/panel_solver.py. The plane stress form of the linear elasticity equations in Modulus can be called by using the LinearElasticityPlaneStress class:

    # make list of nodes to unroll graph on
    le = LinearElasticityPlaneStress(lambda_=lambda_, mu=mu)
    elasticity_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("sigma_hoop")],
        output_keys=[
            Key("u"),
            Key("v"),
            Key("sigma_xx"),
            Key("sigma_yy"),
            Key("sigma_xy"),
        ],
        cfg=cfg.arch.fully_connected,
    )
    nodes = le.make_nodes() + [
        elasticity_net.make_node(name="elasticity_network", jit=cfg.jit)
    ]

Fig. 54 shows the Modulus results for panel displacements and stresses. For comparison, the commercial solver results are also presented in Fig. 55. The Modulus and commercial solver results are in close agreement, with a difference of less than 5% in the maximum Von Mises stress.

Fig. 54 Modulus linear elasticity results for the aircraft fuselage panel example with parameterized hoop stress. The results are for \(\sigma_{hoop}\) = 46

Fig. 55 Commercial solver linear elasticity results for the aircraft fuselage panel example with \(\sigma_{hoop}\) = 46

Linear Elasticity in the Variational Form

Linear elasticity equations in the variational form

In addition to the differential form of the linear elasticity equations, Modulus also enables the use of variational form of these equations (Interface Problem by Variational Method). In this section, we derive the linear elasticity equations in the variational form as implemented in Modulus. We will use the non-dimensionalized variables in this derivation. Start by forming the inner product of Equation (142) and a vector test function \(v \in \mathbf{V}\), and integrating over the domain, as follows

(159)\[\int_{\Omega} \hat{\sigma}_{ji,j} v_i d \mathbf{x}+ \int_{\Omega} \hat{f_i} v_i d \mathbf{x} = 0.\]

Using the integration by parts, we have

(160)\[ \int_{\partial \Omega} \hat{T_i} v_i d \mathbf{s} -\int_{\Omega} \hat{\sigma}_{ji} v_{j,i} d \mathbf{x} + \int_{\Omega} \hat{f_i} v_i d \mathbf{x} = 0,\]

where \(T_i\) is the traction. The first term is zero for the traction-free boundaries. Equation (160) is the variational form of the linear elasticity equations that is adopted in Modulus. The term \(\hat{\sigma}_{ji}\) in this equation is computed using the stress-displacement relation in Equation (153).

Problem 3: Plane displacement

This example solves the linear elasticity plane stress equations in the variational form. Consider a square plate that is clamped from one end and is under a displacement boundary condition on the other end, as shown in Fig. 56. The rest of the boundaries are traction-free. The material properties are assumed to be \((E, \nu) = (10 \,\text{MPa}, 0.2)\). This example is adopted from 1.

Case Setup and Results

To solve this problem in the variational form, we enforce the displacement boundary conditions in the differential form, and also generate interior and boundary points to be used in evaluation of the integrals in Equation (160). The python script for this problem can be found at examples/plane_displacement/plane_displacement.py.

Fig. 56 Geometry and boundary conditions of the plane displacement example. This example is adopted from 1.

    # make domain
    domain = Domain()
    bottomBC = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0.0, "v": 0.0},
        batch_size=cfg.batch_size.bottom,
        batch_per_epoch=5000,
        lambda_weighting={"u": 10.0, "v": 10.0},
        criteria=Eq(y, domain_origin[1]),
    )
    domain.add_constraint(bottomBC, "bottomBC_differential")

    topBC = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0.0, "v": 0.1},
        batch_size=cfg.batch_size.top,
        batch_per_epoch=5000,
        lambda_weighting={"u": 10.0, "v": 10.0},
        criteria=Eq(y, domain_origin[1] + domain_dim[1])
        & (x <= domain_origin[0] + domain_dim[0] / 2.0),
    )
    domain.add_constraint(topBC, "topBC_differential")

    # register variational data
    batch_per_epoch = 1
    variational_datasets = {}
    # bottomBC, index : 0
    invar = geo.sample_boundary(
        batch_per_epoch * cfg.batch_size.bottom,
        criteria=Eq(y, domain_origin[1]),
        quasirandom=True,
    )
    invar["area"] *= batch_per_epoch
    variational_datasets["bottom_bc"] = VariationalDataset(
        batch_size=cfg.batch_size.bottom,
        invar=invar,
        outvar_names=["u__x", "u__y", "v__x", "v__y"],
    )

    # topBC, index : 1
    invar = geo.sample_boundary(
        batch_per_epoch * cfg.batch_size.top,
        criteria=Eq(y, domain_origin[1] + domain_dim[1])
        & (x <= domain_origin[0] + domain_dim[0] / 2.0),
        quasirandom=True,
    )
    invar["area"] *= batch_per_epoch
    variational_datasets["top_bc"] = VariationalDataset(
        batch_size=cfg.batch_size.top,
        invar=invar,
        outvar_names=["u__x", "u__y", "v__x", "v__y"],
    )

    # Interior, index : 2
    invar = geo.sample_interior(
        batch_per_epoch * cfg.batch_size.interior,
        bounds={x: bounds_x, y: bounds_y},
        quasirandom=True,
    )
    invar["area"] *= batch_per_epoch
    variational_datasets["interior"] = VariationalDataset(
        batch_size=cfg.batch_size.interior,
        invar=invar,
        outvar_names=["u__x", "u__y", "v__x", "v__y"],
    )

    # make variational constraints
    variational_constraint = VariationalConstraint(
        datasets=variational_datasets,
        nodes=nodes,
        num_workers=1,
        loss=DGLoss(),
    )
    domain.add_constraint(variational_constraint, "variational")

The displacement boundary conditions have been included in the regular PINN loss function. For the variational constraints, we first created VariationalDataset for each bc, by specifying the invar and the required outvar_names. These outvars will then be used while computing the variational loss. A VariationalConstraint is then constructed using the dictionary of variational_datasets and the nodes of the network. For loss, specify the DGLoss() which is a custom loss function that includes the variational loss. The remainder of this subsection, covers how to generate this variational loss (DGLoss()). First, read the neural network solution and the gradients at the interior and boundary points:

class DGLoss(Loss):
    def __init__(self):
        super().__init__()
        test_fn = Test_Function(
            name_ord_dict={
                Legendre_test: [k for k in range(10)],
                Trig_test: [k for k in range(10)],
            },
            box=[
                [domain_origin[0], domain_origin[1]],
                [domain_origin[0] + domain_dim[0], domain_origin[1] + domain_dim[1]],
            ],
            diff_list=["grad"],
        )

        self.v = Vector_Test(test_fn, test_fn, mix=0.02)

    def forward(
        self,
        list_invar,
        list_outvar,
        step: int,
    ):
        torch.cuda.nvtx.range_push("Make_DGLoss")
        torch.cuda.nvtx.range_push("Make_DGLoss_Get_Data")
        # self.v.sample_vector_test()
        # get points on the interior
        x_interior = list_invar[2]["x"]
        y_interior = list_invar[2]["y"]
        area_interior = list_invar[2]["area"]

        # compute solution for the interior
        u_x_interior = list_outvar[2]["u__x"]
        u_y_interior = list_outvar[2]["u__y"]
        v_x_interior = list_outvar[2]["v__x"]
        v_y_interior = list_outvar[2]["v__y"]

        # get points on the boundary
        x_bottom_dir = list_invar[0]["x"]
        y_bottom_dir = list_invar[0]["y"]
        normal_x_bottom_dir = list_invar[0]["normal_x"]
        normal_y_bottom_dir = list_invar[0]["normal_y"]
        area_bottom_dir = list_invar[0]["area"]

        x_top_dir = list_invar[1]["x"]
        y_top_dir = list_invar[1]["y"]
        normal_x_top_dir = list_invar[1]["normal_x"]
        normal_y_top_dir = list_invar[1]["normal_y"]
        area_top_dir = list_invar[1]["area"]

        # compute solution for the boundary
        u_x_bottom_dir = list_outvar[0]["u__x"]
        u_y_bottom_dir = list_outvar[0]["u__y"]
        v_x_bottom_dir = list_outvar[0]["v__x"]
        v_y_bottom_dir = list_outvar[0]["v__y"]

        u_x_top_dir = list_outvar[1]["u__x"]
        u_y_top_dir = list_outvar[1]["u__y"]
        v_x_top_dir = list_outvar[1]["v__x"]
        v_y_top_dir = list_outvar[1]["v__y"]
        torch.cuda.nvtx.range_pop()
        torch.cuda.nvtx.range_push("Make_DGLoss_Test_Function")

In the next step, define the test functions, and compute the test functions and their required gradients at the interior and boundary points:

        # test functions
        vx_x_interior, vy_x_interior = self.v.eval_test("vx", x_interior, y_interior)
        vx_y_interior, vy_y_interior = self.v.eval_test("vy", x_interior, y_interior)
        vx_bottom_dir, vy_bottom_dir = self.v.eval_test("v", x_bottom_dir, y_bottom_dir)
        vx_top_dir, vy_top_dir = self.v.eval_test("v", x_top_dir, y_top_dir)

Here, a set of test functions consisting of Legendre polynomials and trigonometric functions is constructed, and 2% of these test functions are randomly sampled. Compute only the terms that appear in the variational loss in Equation (160). For instance, it is not necessary to compute the derivative of the test functions with respect to input coordinates for the boundary points.

The next step is to compute the stress terms according to the plane stress equations in Equations (157), (158), and also the traction terms:

        torch.cuda.nvtx.range_pop()
        torch.cuda.nvtx.range_push("Make_DGLoss_Computation")
        w_z_interior = -lambda_ / (lambda_ + 2 * mu) * (u_x_interior + v_y_interior)
        sigma_xx_interior = (
            lambda_ * (u_x_interior + v_y_interior + w_z_interior)
            + 2 * mu * u_x_interior
        )
        sigma_yy_interior = (
            lambda_ * (u_x_interior + v_y_interior + w_z_interior)
            + 2 * mu * v_y_interior
        )
        sigma_xy_interior = mu * (u_y_interior + v_x_interior)

        w_z_bottom_dir = (
            -lambda_ / (lambda_ + 2 * mu) * (u_x_bottom_dir + v_y_bottom_dir)
        )
        sigma_xx_bottom_dir = (
            lambda_ * (u_x_bottom_dir + v_y_bottom_dir + w_z_bottom_dir)
            + 2 * mu * u_x_bottom_dir
        )
        sigma_yy_bottom_dir = (
            lambda_ * (u_x_bottom_dir + v_y_bottom_dir + w_z_bottom_dir)
            + 2 * mu * v_y_bottom_dir
        )
        sigma_xy_bottom_dir = mu * (u_y_bottom_dir + v_x_bottom_dir)

        w_z_top_dir = -lambda_ / (lambda_ + 2 * mu) * (u_x_top_dir + v_y_top_dir)
        sigma_xx_top_dir = (
            lambda_ * (u_x_top_dir + v_y_top_dir + w_z_top_dir) + 2 * mu * u_x_top_dir
        )
        sigma_yy_top_dir = (
            lambda_ * (u_x_top_dir + v_y_top_dir + w_z_top_dir) + 2 * mu * v_y_top_dir
        )
        sigma_xy_top_dir = mu * (u_y_top_dir + v_x_top_dir)

        traction_x_bottom_dir = (
            sigma_xx_bottom_dir * normal_x_bottom_dir
            + sigma_xy_bottom_dir * normal_y_bottom_dir
        )
        traction_y_bottom_dir = (
            sigma_xy_bottom_dir * normal_x_bottom_dir
            + sigma_yy_bottom_dir * normal_y_bottom_dir
        )
        traction_x_top_dir = (
            sigma_xx_top_dir * normal_x_top_dir + sigma_xy_top_dir * normal_y_top_dir
        )
        traction_y_top_dir = (
            sigma_xy_top_dir * normal_x_top_dir + sigma_yy_top_dir * normal_y_top_dir
        )

Finally, following the Equation (160), define your variational interior and boundary integral terms, formulate the variational loss, and add that to the overall loss, as follows:

        torch.cuda.nvtx.range_pop()
        torch.cuda.nvtx.range_push("Make_DGLoss_Integral")
        interior_loss = tensor_int(
            area_interior,
            sigma_xx_interior * vx_x_interior
            + sigma_yy_interior * vy_y_interior
            + sigma_xy_interior * (vx_y_interior + vy_x_interior),
        )
        boundary_loss1 = tensor_int(
            area_bottom_dir,
            traction_x_bottom_dir * vx_bottom_dir
            + traction_y_bottom_dir * vy_bottom_dir,
        )
        boundary_loss2 = tensor_int(
            area_top_dir,
            traction_x_top_dir * vx_top_dir + traction_y_top_dir * vy_top_dir,
        )
        torch.cuda.nvtx.range_pop()
        torch.cuda.nvtx.range_push("Make_DGLoss_Register_Loss")
        losses = {
            "variational_plane": torch.abs(
                interior_loss - boundary_loss1 - boundary_loss2
            )
            .pow(2)
            .sum()
        }
        torch.cuda.nvtx.range_pop()
        torch.cuda.nvtx.range_pop()
        return losses

Fig. 57 shows the Modulus results for this plane displacement example. The results are in good agreement with the FEM results reported in 1, verifying the accuracy of the Modulus results in the variational form.

Fig. 57 Modulus results for the plane displacement example.

References

1(1,2,3,4)

Chengping Rao, Hao Sun, and Yang Liu. Physics informed deep learning for computational elastodynamics without labeled data. arXiv preprint arXiv:2006.08472, 2020.