rename

test/test_rotation_stretch.py → test/rotation_stretch.py

@@ -40,7 +40,7 @@ def test_rotation_stretch(file_path: str | None = None, stretch: float = 0.5):
 
     """
     E = 2e6
-    nu = 0.0
+    nu = 0.3
     mu = E / (2 * (1 + nu))
     lam = E * nu / ((1 + nu) * (1 - 2 * nu))
     parameters = {"mu": mu, "lambda": lam, "density": 1.0}
@@ -299,25 +299,143 @@ def edge_1_1_z(x):
     plt.show()
 
 
-if __name__ == "__main__":
-    test_rotation("output.xdmf")
-    test_rotation_stretch("output.xdmf", stretch=0.01)
+def test_pure_stretch(file_path: str | None = None, stretch: float = 0.5):
+    """
+      ---------
+     /|       /|
+    / |      / |
+    --|------  |
+    | |_____|  |
+    | /     | /
+    |/______|/
 
-    import matplotlib.pyplot as plt
+    """
+    E = 2e6
+    nu = 0.3
+    mu = E / (2 * (1 + nu))
+    lam = E * nu / ((1 + nu) * (1 - 2 * nu))
+    parameters = {"mu": mu, "lambda": lam, "density": 1.0}
+
+    mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 1, 1, 1, df.mesh.CellType.hexahedron)
 
-    stretch = 0.5
-    t = np.linspace(0, 100, 10)
-    circle_left = np.array([unit_circle_left(ti, 100) for ti in t]).reshape(-1, 3)
-    circle_right = np.array([unit_circle_right(ti, 100) for ti in t]).reshape(-1, 3)
-    circle_middle = np.array(
-        [2**0.5 * unit_circle_middle(ti, 100) * (1 + stretch * (ti / 100)) for ti in t]
-    ).reshape(-1, 3)
-    stretched_circle_right = np.array([unit_circle_right(ti, 100) * (1 + stretch * (ti / 100)) for ti in t]).reshape(
-        -1, 3
+    rule = co.helpers.QuadratureRule(
+        quadrature_type=basix.QuadratureType.Default,
+        cell_type=basix.CellType.hexahedron,
+        degree=1,
     )
-    plt.plot(circle_left[:, 0], circle_left[:, 1], label="left")
-    plt.plot(circle_right[:, 0], circle_right[:, 1], label="right")
-    plt.plot(stretched_circle_right[:, 0], stretched_circle_right[:, 1], label="stretched right")
-    plt.plot(circle_middle[:, 0], circle_middle[:, 1], label="middle")
-    plt.legend()
+
+    tdim = mesh.topology.dim
+    fdim = tdim - 1
+    mesh.topology.create_connectivity(fdim, tdim)
+
+    P1 = df.fem.VectorFunctionSpace(mesh, ("Lagrange", 1))
+    u = df.fem.Function(P1)
+    t_end = 100.0
+
+    def edge_0_0_z(x):
+        return np.isclose(x[0], 0.0) & np.isclose(x[1], 0.0)
+
+    def edge_0_1_z(x):
+        return np.isclose(x[0], 0.0) & np.isclose(x[1], 1.0)
+
+    def edge_1_0_z(x):
+        return np.isclose(x[0], 1.0) & np.isclose(x[1], 0.0)
+
+    def edge_1_1_z(x):
+        return np.isclose(x[0], 1.0) & np.isclose(x[1], 1.0)
+
+    dofs_0_0_z = df.fem.locate_dofs_geometrical(P1, edge_0_0_z)
+    dofs_0_1_z = df.fem.locate_dofs_geometrical(P1, edge_0_1_z)
+    dofs_1_0_z = df.fem.locate_dofs_geometrical(P1, edge_1_0_z)
+    dofs_1_1_z = df.fem.locate_dofs_geometrical(P1, edge_1_1_z)
+
+    value_0_0_z = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0]))
+    value_0_1_z = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0]))
+    value_1_0_z = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0]))
+    value_1_1_z = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0]))
+
+    bcs = [
+        df.fem.dirichletbc(value_0_0_z, dofs_0_0_z, P1),
+        df.fem.dirichletbc(value_0_1_z, dofs_0_1_z, P1),
+        df.fem.dirichletbc(value_1_0_z, dofs_1_0_z, P1),
+        df.fem.dirichletbc(value_1_1_z, dofs_1_1_z, P1),
+    ]
+
+    law = co.laws.PyHypoelasticity3D(parameters)
+
+    v_ = ufl.TestFunction(P1)
+    u_ = ufl.TrialFunction(P1)
+
+    mass_form = ufl.inner(u_, v_) * parameters["density"] * ufl.dx
+
+    M = df.fem.petsc.assemble_matrix(df.fem.form(mass_form))
+    M.assemble()
+    ones = df.fem.Function(P1)
+    with ones.vector.localForm() as ones_local:
+        ones_local.set(1.0)
+    M_action = M * ones.vector
+    M_function = df.fem.Function(P1)
+    M_function.vector.array[:] = M_action.array
+
+    M_action.array[:] = 1.0 / M_action.array
+    M_action.ghostUpdate()
+
+    solver = co.cdm.CDM3D(
+        P1,
+        0,
+        None,
+        bcs,
+        M_function,
+        law,
+        rule,
+    )
+    stress = []
+    t = []
+
+    h = 1e-1
+    # n_steps = int(t_end / h)
+    # h = t_end / n_steps
+    v_x = stretch / t_end
+    value_1_0_z.value = np.array([v_x, 0.0, 0.0])
+    value_1_1_z.value = np.array([v_x, 0.0, 0.0])
+    with df.io.XDMFFile(MPI.COMM_WORLD, f"{file_path}output.xdmf", "w") as xdmf:
+        xdmf.write_mesh(mesh)
+    while solver.t < t_end:
+
+        solver.step(h)
+        t.append(solver.t)
+        stress.append(solver.q_fields["mandel_stress"].x.array.copy())
+        # with df.io.XDMFFile(MPI.COMM_WORLD, f"{file_path}output.xdmf", "a") as xdmf:
+        #    xdmf.write_function(solver.fields["u"], solver.t)
+    import matplotlib.pyplot as plt
+
+    stress = np.array(stress).reshape(-1, 6)
+    plt.plot(t, stress[:, 0])
+    plt.plot(t, stress[:, 1])
+    plt.plot(t, stress[:, 2])
     plt.show()
+
+
+if __name__ == "__main__":
+    test_rotation("output.xdmf")
+    test_rotation_stretch("output.xdmf", stretch=0.01)
+    test_pure_stretch("output.xdmf", stretch=0.01)
+
+    # import matplotlib.pyplot as plt
+
+    # stretch = 0.5
+    # t = np.linspace(0, 100, 10)
+    # circle_left = np.array([unit_circle_left(ti, 100) for ti in t]).reshape(-1, 3)
+    # circle_right = np.array([unit_circle_right(ti, 100) for ti in t]).reshape(-1, 3)
+    # circle_middle = np.array(
+    #     [2**0.5 * unit_circle_middle(ti, 100) * (1 + stretch * (ti / 100)) for ti in t]
+    # ).reshape(-1, 3)
+    # stretched_circle_right = np.array([unit_circle_right(ti, 100) * (1 + stretch * (ti / 100)) for ti in t]).reshape(
+    #     -1, 3
+    # )
+    # plt.plot(circle_left[:, 0], circle_left[:, 1], label="left")
+    # plt.plot(circle_right[:, 0], circle_right[:, 1], label="right")
+    # plt.plot(stretched_circle_right[:, 0], stretched_circle_right[:, 1], label="stretched right")
+    # plt.plot(circle_middle[:, 0], circle_middle[:, 1], label="middle")
+    # plt.legend()
+    # plt.show()