[Feature&Fix] Add sympy_to_func

ppsci/equation/pde/linear_elasticity.py

@@ -33,9 +33,6 @@ class LinearElasticity(base.PDE):
         traction_{x} = n_x \sigma_{xx} + n_y \sigma_{xy} + n_z \sigma_{xz} \\
         traction_{y} = n_y \sigma_{yx} + n_y \sigma_{yy} + n_z \sigma_{yz} \\
         traction_{z} = n_z \sigma_{zx} + n_y \sigma_{zy} + n_z \sigma_{zz} \\
-        navier_{x} = \rho(\dfrac{\partial^2 u}{\partial t^2}) - (\lambda + \mu)(\dfrac{\partial^2 u}{\partial x^2}+\dfrac{\partial^2 v}{\partial y \partial x} + \dfrac{\partial^2 w}{\partial z \partial x}) - \mu(\dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2} + \dfrac{\partial^2 u}{\partial z^2}) \\
-        navier_{y} = \rho(\dfrac{\partial^2 v}{\partial t^2}) - (\lambda + \mu)(\dfrac{\partial^2 v}{\partial x \partial y}+\dfrac{\partial^2 v}{\partial y^2} + \dfrac{\partial^2 w}{\partial z \partial y}) - \mu(\dfrac{\partial^2 v}{\partial x^2} + \dfrac{\partial^2 v}{\partial y^2} + \dfrac{\partial^2 v}{\partial z^2}) \\
-        navier_{z} = \rho(\dfrac{\partial^2 w}{\partial t^2}) - (\lambda + \mu)(\dfrac{\partial^2 w}{\partial x \partial z}+\dfrac{\partial^2 v}{\partial y \partial z} + \dfrac{\partial^2 w}{\partial z^2}) - \mu(\dfrac{\partial^2 w}{\partial x^2} + \dfrac{\partial^2 w}{\partial y^2} + \dfrac{\partial^2 w}{\partial z^2}) \\
     \end{cases}
     $$
 
@@ -66,6 +63,12 @@ def __init__(
         time: bool = False,
     ):
         super().__init__()
+        if lambda_ is None:
+            nu = float(nu)
+            E = float(E)
+            lambda_ = nu * E / ((1 + nu) * (1 - 2 * nu))
+            mu = E / (2 * (1 + nu))
+
         self.E = E
         self.nu = nu
         self.lambda_ = lambda_
@@ -76,13 +79,11 @@ def __init__(
 
         # Stress equations
         def stress_disp_xx_compute_func(out):
-            x, y, z, u, v, w = (
+            x, y, u, v = (
                 out["x"],
                 out["y"],
-                out["z"],
                 out["u"],
                 out["v"],
-                out["w"],
             )
             sigma_xx = out["sigma_xx"]
             stress_disp_xx = (
@@ -98,13 +99,11 @@ def stress_disp_xx_compute_func(out):
         self.add_equation("stress_disp_xx", stress_disp_xx_compute_func)
 
         def stress_disp_yy_compute_func(out):
-            x, y, z, u, v, w = (
+            x, y, u, v = (
                 out["x"],
                 out["y"],
-                out["z"],
                 out["u"],
                 out["v"],
-                out["w"],
             )
             sigma_yy = out["sigma_yy"]
             stress_disp_yy = (
@@ -168,7 +167,7 @@ def stress_disp_yz_compute_func(out):
 
         # Equations of equilibrium
         def equilibrium_x_compute_func(out):
-            x, y, z = out["x"], out["y"], out["z"]
+            x, y = out["x"], out["y"]
             sigma_xx, sigma_xy = out["sigma_xx"], out["sigma_xy"]
             equilibrium_x = -jacobian(sigma_xx, x) - jacobian(sigma_xy, y)
             if self.dim == 3:
@@ -182,11 +181,10 @@ def equilibrium_x_compute_func(out):
         self.add_equation("equilibrium_x", equilibrium_x_compute_func)
 
         def equilibrium_y_compute_func(out):
-            x, y, z = out["x"], out["y"], out["z"]
-            sigma_xy, sigma_yy, sigma_yz = (
+            x, y = out["x"], out["y"]
+            sigma_xy, sigma_yy = (
                 out["sigma_xy"],
                 out["sigma_yy"],
-                out["sigma_yz"],
             )
             equilibrium_y = -jacobian(sigma_xy, x) - jacobian(sigma_yy, y)
             if self.dim == 3:
@@ -255,89 +253,22 @@ def traction_y_compute_func(out):
 
         self.add_equation("traction_y", traction_y_compute_func)
 
-        def traction_z_compute_func(out):
-            normal_x, normal_y, normal_z = (
-                out["normal_x"],
-                out["normal_y"],
-                out["normal_z"],
-            )
-            sigma_xz, sigma_yz, sigma_zz = (
-                out["sigma_xz"],
-                out["sigma_yz"],
-                out["sigma_zz"],
-            )
-            traction_z = normal_x * sigma_xz + normal_y * sigma_yz + normal_z * sigma_zz
-            return traction_z
-
-        self.add_equation("traction_z", traction_z_compute_func)
-
-        # Navier equations
-        def navier_x_compute_func(out):
-            x, y, u, v = (
-                out["x"],
-                out["y"],
-                out["u"],
-                out["v"],
-            )
-            duxvywz = jacobian(u, x) + jacobian(v, y)
-            duxxuyyuzz = hessian(u, x) + hessian(u, y)
-            if self.dim == 3:
-                z, w = out["z"], out["w"]
-                duxvywz += jacobian(w, z)
-                duxxuyyuzz += hessian(u, z)
-            navier_x = (
-                -(self.lambda_ + self.mu) * jacobian(duxvywz, x) - self.mu * duxxuyyuzz
-            )
-            if self.time:
-                t = out["t"]
-                navier_x += rho * hessian(u, t)
-            return navier_x
-
-        self.add_equation("navier_x", navier_x_compute_func)
-
-        def navier_y_compute_func(out):
-            x, y, u, v = (
-                out["x"],
-                out["y"],
-                out["u"],
-                out["v"],
-            )
-            duxvywz = jacobian(u, x) + jacobian(v, y)
-            dvxxvyyvzz = hessian(v, x) + hessian(v, y)
-            if self.dim == 3:
-                z, w = out["z"], out["w"]
-                duxvywz += jacobian(w, z)
-                dvxxvyyvzz += hessian(v, z)
-            navier_y = (
-                -(self.lambda_ + self.mu) * jacobian(duxvywz, y) - self.mu * dvxxvyyvzz
-            )
-            if self.time:
-                t = out["t"]
-                navier_y += rho * hessian(v, t)
-            return navier_y
-
-        self.add_equation("navier_y", navier_y_compute_func)
-
         if self.dim == 3:
 
-            def navier_z_compute_func(out):
-                x, y, z, u, v, w = (
-                    out["x"],
-                    out["y"],
-                    out["z"],
-                    out["u"],
-                    out["v"],
-                    out["w"],
+            def traction_z_compute_func(out):
+                normal_x, normal_y, normal_z = (
+                    out["normal_x"],
+                    out["normal_y"],
+                    out["normal_z"],
                 )
-                duxvywz = jacobian(u, x) + jacobian(v, y) + jacobian(w, z)
-                dwxxvyyvzz = hessian(w, x) + hessian(w, y) + hessian(w, z)
-                navier_z = (
-                    -(self.lambda_ + self.mu) * jacobian(duxvywz, z)
-                    - self.mu * dwxxvyyvzz
+                sigma_xz, sigma_yz, sigma_zz = (
+                    out["sigma_xz"],
+                    out["sigma_yz"],
+                    out["sigma_zz"],
                 )
-                if self.time:
-                    t = out["t"]
-                    navier_z += rho * hessian(w, t)
-                return navier_z
+                traction_z = (
+                    normal_x * sigma_xz + normal_y * sigma_yz + normal_z * sigma_zz
+                )
+                return traction_z
 
-            self.add_equation("navier_z", navier_z_compute_func)
+            self.add_equation("traction_z", traction_z_compute_func)