Absolutely! Here’s a step-by-step Python code to **mesh half a rectangular plate with a circular hole** using 4-node quadrilateral elements, with user input for `a` (width) and `r` (hole radius). We'll set up the mesh, boundary conditions, and outline FEM setup (though full FEM assembly and solve is a big task, we'll give a strong foundation). Libraries like `numpy` and `matplotlib` will be used for mesh visualization. **Assumptions:** - Plate dimensions: Length \( L = 2a \), Width \( a \) (since a/r >= 5, plate is much bigger than hole) - Only *right half* is modeled (due to symmetry) - Circular hole is at plate center - UDL on top edge, vertical symmetry on left edge, both sides (left and right) fixed horizontally - Top and bottom edges free to move vertically --- ### 1. **Imports and Inputs** ```python import numpy as np import matplotlib.pyplot as plt # User inputs a = float(input("Enter plate width a (m): ")) r = float(input("Enter hole radius r (m), with a/r >= 5: ")) assert a/r >= 5, "a/r must be >= 5" E = 200e9 # Young's modulus (Pa) v = 0.3 # Poisson's ratio L = 2 * a # Plate length # Mesh parameters nx = 40 # Number of elements along length ny = 20 # Number of elements along width ``` --- ### 2. **Mesh Generation (Half Plate with Hole)** We'll mesh only the right half (`x` from 0 to L/2). ```python def generate_mesh(a, L, r, nx, ny): # Generate grid points x = np.linspace(0, L/2, nx+1) y = np.linspace(-a/2, a/2, ny+1) nodes = [] node_id_map = {} node_id = 0 for j in range(ny+1): for i in range(nx+1): xp, yp = x[i], y[j] # Only include nodes outside the hole if (xp - L/2)**2 + yp**2 >= r**2: nodes.append([xp, yp]) node_id_map[(i, j)] = node_id node_id += 1 else: node_id_map[(i, j)] = None nodes = np.array(nodes) # Generate elements elements = [] for j in range(ny): for i in range(nx): n1 = node_id_map[(i, j)] n2 = node_id_map[(i+1, j)] n3 = node_id_map[(i+1, j+1)] n4 = node_id_map[(i, j+1)] # Only add elements if all nodes exist (outside hole) if None not in (n1, n2, n3, n4): elements.append([n1, n2, n3, n4]) return nodes, elements nodes, elements = generate_mesh(a, L, r, nx, ny) ``` --- ### 3. **Visualization** ```python def plot_mesh(nodes, elements): plt.figure(figsize=(10, 5)) for elem in elements: xy = nodes[elem + [elem[0]]] # close the loop plt.plot(xy[:, 0], xy[:, 1], 'b-') plt.plot(nodes[:, 0], nodes[:, 1], 'ro', markersize=2) plt.gca().set_aspect('equal') plt.title("Mesh of Half Plate with Central Hole") plt.xlabel("x [m]") plt.ylabel("y [m]") plt.show() plot_mesh(nodes, elements) ``` --- ### 4. **Boundary Conditions Outline** - **Left edge (x=0):** symmetry, so \( u_x = 0 \) - **Right edge (x=L/2):** fixed horizontally (\( u_x = 0 \)) - **Top/bottom (y=±a/2):** free in x, can move in y - **UDL** applied on **top edge** (\( y = a/2 \)) ```python # Find node indices for boundary conditions tol = 1e-8 left_nodes = [i for i, (x, y) in enumerate(nodes) if abs(x) < tol] right_nodes = [i for i, (x, y) in enumerate(nodes) if abs(x - L/2) < tol] top_nodes = [i for i, (x, y) in enumerate(nodes) if abs(y - a/2) < tol] bottom_nodes = [i for i, (x, y) in enumerate(nodes) if abs(y + a/2) < tol] print(f"Nodes on left edge (u_x=0): {left_nodes}") print(f"Nodes on right edge (u_x=0): {right_nodes}") print(f"Nodes on top edge (UDL): {top_nodes}") print(f"Nodes on bottom edge: {bottom_nodes}") ``` --- ### 5. **FEM System Setup (Outline)** A full FEM implementation is lengthy. Here’s the outline: ```python # Placeholder for assembling the FEM system num_nodes = nodes.shape[0] dofs = 2 * num_nodes # 2 DOF per node (u,v) K = np.zeros((dofs, dofs)) F = np.zeros(dofs) # 1. Loop over elements, calculate element stiffness matrix (see standard 4-node quad formulation) # 2. Assemble global stiffness matrix K and force vector F # 3. Apply boundary conditions (u_x = 0 at left/right, UDL at top edge nodes) # 4. Solve for displacements: U = np.linalg.solve(K, F) ``` --- ## **Full Script (Mesh + Visualization + BCs):** ```python import numpy as np import matplotlib.pyplot as plt # User inputs a = float(input("Enter plate width a (m): ")) r = float(input("Enter hole radius r (m), with a/r >= 5: ")) assert a/r >= 5, "a/r must be >= 5" E = 200e9 # Young's modulus (Pa) v = 0.3 # Poisson's ratio L = 2 * a # Plate length # Mesh parameters nx = 40 # Number of elements along length ny = 20 # Number of elements along width def generate_mesh(a, L, r, nx, ny): x = np.linspace(0, L/2, nx+1) y = np.linspace(-a/2, a/2, ny+1) nodes = [] node_id_map = {} node_id = 0 for j in range(ny+1): for i in range(nx+1): xp, yp = x[i], y[j] if (xp - L/2)**2 + yp**2 >= r**2: nodes.append([xp, yp]) node_id_map[(i, j)] = node_id node_id += 1 else: node_id_map[(i, j)] = None nodes = np.array(nodes) elements = [] for j in range(ny): for i in range(nx): n1 = node_id_map[(i, j)] n2 = node_id_map[(i+1, j)] n3 = node_id_map[(i+1, j+1)] n4 = node_id_map[(i, j+1)] if None not in (n1, n2, n3, n4): elements.append([n1, n2, n3, n4]) return nodes, elements nodes, elements = generate_mesh(a, L, r, nx, ny) def plot_mesh(nodes, elements): plt.figure(figsize=(10, 5)) for elem in elements: xy = nodes[elem + [elem[0]]] # close the loop plt.plot(xy[:, 0], xy[:, 1], 'b-') plt.plot(nodes[:, 0], nodes[:, 1], 'ro', markersize=2) plt.gca().set_aspect('equal') plt.title("Mesh of Half Plate with Central Hole") plt.xlabel("x [m]") plt.ylabel("y [m]") plt.show() plot_mesh(nodes, elements) # Boundary conditions tol = 1e-8 left_nodes = [i for i, (x, y) in enumerate(nodes) if abs(x) < tol] right_nodes = [i for i, (x, y) in enumerate(nodes) if abs(x - L/2) < tol] top_nodes = [i for i, (x, y) in enumerate(nodes) if abs(y - a/2) < tol] bottom_nodes = [i for i, (x, y) in enumerate(nodes) if abs(y + a/2) < tol] print(f"Nodes on left edge (u_x=0): {left_nodes}") print