diff --git a/.gitignore b/.gitignore index be120d8..dae3ea7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ .idea/ -CAD_2plates.step \ No newline at end of file +CAD_2plates.step +.venv +__pycache__ \ No newline at end of file diff --git a/src/autobolt/combined.py b/src/autobolt/combined.py index ce6aa57..3652479 100644 --- a/src/autobolt/combined.py +++ b/src/autobolt/combined.py @@ -1,6 +1,50 @@ from .fea_solver import _calculate_fos_from_build123d -from .parametric_cad_solver import create_two_plate_assembly +from .parametric_cad_solver import create_two_plate_assembly, create_two_plate_assembly_components +from .fea_solver_dual import _calculate_fos_from_build123d_dual + + +def calculate_fos_dual( + plate_length_m: float, + plate_width_m: float, + plate_thickness_m: float, + num_holes: int, + hole_radius_m: float, + edge_margin_m: float, + hole_spacing_m: float, + hole_offset_from_bottom_m: float, + plate_gap_mm: float, + elastic_modulus: float, + poissons_ratio: float, + bolt_yield_strength: float, + plate_yield_strength: float, + traction_values: list[tuple], +) -> tuple[float, float]: + + plateA, plateB, bolt_cylinders = create_two_plate_assembly( + plate_length_m, + plate_width_m, + plate_thickness_m, + num_holes, + hole_radius_m, + edge_margin_m, + hole_spacing_m, + hole_offset_from_bottom_m, + plate_gap_mm, + ) + + bolt_fos, plate_fos = _calculate_fos_from_build123d_dual( + plateA, + plateB, + bolt_cylinders, + elastic_modulus, + poissons_ratio, + bolt_yield_strength, + plate_yield_strength, + traction_values, + ) + + return bolt_fos, plate_fos def calculate_fos( plate_length_m: float, diff --git a/src/autobolt/fea_solver_dual.py b/src/autobolt/fea_solver_dual.py new file mode 100644 index 0000000..fc8d668 --- /dev/null +++ b/src/autobolt/fea_solver_dual.py @@ -0,0 +1,224 @@ +import tempfile + +import build123d +import fenics +import gmsh +import meshio + +import numpy as np + +BOLT_TAG = 101 +PLATE_TAG = 102 + +def import_step_and_capture_new_volumes(step_path: str) -> list[int]: + """Import a STEP file and return the list of volume entity tags created by that import.""" + before = set(t for _, t in gmsh.model.getEntities(3)) + gmsh.merge(step_path) + gmsh.model.occ.synchronize() + after = set(t for _, t in gmsh.model.getEntities(3)) + new_tags = sorted(after - before) + return new_tags + +def _run_fenics_simulation( + xdmf_file_mesh: str, + xdmf_file_surface_tags: str, + elastic_modulus: float, + poissons_ratio: float, + bolt_yield_strength: float, + plate_yield_strength: float, + traction_values: list[tuple], + zero_disp_tags: list[int], + traction_tags: list[int] + ): + mesh = fenics.cpp.mesh.Mesh() + with fenics.cpp.io.XDMFFile(xdmf_file_mesh) as infile: + infile.read(mesh) + + mvc_surf = fenics.MeshValueCollection("size_t", mesh, 2) + with fenics.cpp.io.XDMFFile(xdmf_file_surface_tags) as infile: + infile.read(mvc_surf, "name_to_read") + mf = fenics.cpp.mesh.MeshFunctionSizet(mesh, mvc_surf) + + mvc_cell = fenics.MeshValueCollection("size_t", mesh, 3) + with fenics.XDMFFile(xdmf_file_mesh) as infile: + infile.read(mvc_cell, "name_to_read") + cf = fenics.MeshFunction("size_t", mesh, mvc_cell) + + V = fenics.VectorFunctionSpace(mesh, "CG", 2) + + bcs = [ + fenics.DirichletBC(V, fenics.Constant((0.0, 0.0, 0.0)), mf, tag) + for tag in zero_disp_tags + ] + + mu = elastic_modulus / (2.0 * (1.0 + poissons_ratio)) + lmbda = elastic_modulus * poissons_ratio / ( + (1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio) + ) + + def epsilon(u): + return 0.5 * (fenics.grad(u) + fenics.grad(u).T) + + def sigma(u): + return 2.0 * mu * epsilon(u) + lmbda * fenics.tr(epsilon(u)) * fenics.Identity(3) + + u = fenics.TrialFunction(V) + v = fenics.TestFunction(V) + f = fenics.Constant((0.0, 0.0, 0.0)) + + a = fenics.inner(sigma(u), epsilon(v)) * fenics.dx + ds = fenics.Measure("ds", domain=mesh, subdomain_data=mf) + L = fenics.dot(f, v) * fenics.dx + + for tag, traction in zip(traction_tags, traction_values): + L += fenics.dot(fenics.Constant(traction), v) * ds(tag) + + u_sol = fenics.Function(V) + fenics.solve(a == L, u_sol, bcs) + + s_dev = sigma(u_sol) - (1.0 / 3.0) * fenics.tr(sigma(u_sol)) * fenics.Identity(3) + W0 = fenics.FunctionSpace(mesh, "DG", 0) + vm0 = fenics.project(fenics.sqrt(3.0 / 2.0 * fenics.inner(s_dev, s_dev)), W0) + vals = vm0.vector().get_local() + markers = cf.array() + + bolt_vals = vals[markers == BOLT_TAG] + plate_vals = vals[markers == PLATE_TAG] + + max_bolt_vm = float(np.max(bolt_vals)) + max_plate_vm = float(np.max(plate_vals)) + + bolt_fos = float(bolt_yield_strength / max_bolt_vm) + plate_fos = float(plate_yield_strength / max_plate_vm) + + return bolt_fos, plate_fos + +def _calculate_fos_from_build123d_dual( + plateA: build123d.Shape, + plateB: build123d.Shape, + bolt_cylinders: build123d.Shape, + elastic_modulus: float, + poissons_ratio: float, + bolt_yield_strength: float, + plate_yield_strength: float, + traction_values: list[tuple], +): + tempdir = tempfile.TemporaryDirectory() + + xdmf_file_mesh = tempdir.name + "/mesh.xdmf" + xdmf_file_surface_tags = tempdir.name + "surface_tags.xdmf" + mesh_file, zero_disp_tags, traction_tags = generate_tagged_mesh(plateA, plateB, bolt_cylinders, tempdir) + + convert_to_xdmf(mesh_file, "tetra", xdmf_file_mesh) + convert_to_xdmf(mesh_file, "triangle", xdmf_file_surface_tags) + + return _run_fenics_simulation( + xdmf_file_mesh, + xdmf_file_surface_tags, + elastic_modulus, + poissons_ratio, + bolt_yield_strength, + plate_yield_strength, + traction_values, + zero_disp_tags, + traction_tags + ) + +def convert_to_xdmf(mesh_file, cell_type, out_path): + mesh = meshio.read(mesh_file) + + cells = mesh.get_cells_type(cell_type) + cell_data = mesh.get_cell_data("gmsh:physical", cell_type) + out = meshio.Mesh( + points=mesh.points, + cells={cell_type: cells}, + cell_data={"name_to_read": [cell_data]}, + ) + meshio.write(out_path, out) + +def generate_tagged_mesh( + plateA: build123d.Shape, + plateB: build123d.Shape, + bolt_cylinders: build123d.Shape, + tempdir: tempfile.TemporaryDirectory +): + gmsh.initialize() + + step_plateA = tempdir.name + "/plateA.step" + step_plateB = tempdir.name + "/plateB.step" + step_bolts = tempdir.name + "/bolts.step" + + build123d.export_step(plateA, step_plateA) + build123d.export_step(plateB, step_plateB) + build123d.export_step(bolt_cylinders, step_bolts) + + plateA_vols = import_step_and_capture_new_volumes(step_plateA) + plateB_vols = import_step_and_capture_new_volumes(step_plateB) + bolt_vols = import_step_and_capture_new_volumes(step_bolts) + + in_dimtags = [(3, t) for t in (plateA_vols + plateB_vols + bolt_vols)] + _, out_map = gmsh.model.occ.fragment(in_dimtags, []) + gmsh.model.occ.synchronize() + + gmsh.model.occ.removeAllDuplicates() + gmsh.model.occ.synchronize() + + bolt_set_pre = set(bolt_vols) + + bolt_frag_vols: set[int] = set() + plate_frag_vols: set[int] = set() + + for (src_dim, src_tag), frags in zip(in_dimtags, out_map): + frag_vols = [t for dim, t in frags if dim == 3] + if src_tag in bolt_set_pre: + bolt_frag_vols.update(frag_vols) + else: + plate_frag_vols.update(frag_vols) + + gmsh.model.addPhysicalGroup(3, sorted(bolt_frag_vols), tag=BOLT_TAG) + gmsh.model.setPhysicalName(3, BOLT_TAG, "bolt") + gmsh.model.addPhysicalGroup(3, sorted(plate_frag_vols), tag=PLATE_TAG) + gmsh.model.setPhysicalName(3, PLATE_TAG, "plate") + + surfaces = gmsh.model.getEntities(2) + + for (sdim, stag) in surfaces: + gmsh.model.addPhysicalGroup(sdim, [stag], tag=stag) + + surface_info = [] + for dim, s in gmsh.model.getEntities(2): + xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.occ.getBoundingBox(dim, s) + surface_info.append((s, ymin, ymax)) + + global_ymin = min(ymin for _, ymin, _ in surface_info) + global_ymax = max(ymax for _, _, ymax in surface_info) + + def y_face(target_y: float, tol: float = 1e-6) -> int: + for tag, ymin, ymax in surface_info: + if abs(ymax - ymin) < tol and abs(ymin - target_y) < tol: + return tag + raise RuntimeError(f"Could not find planar Y-face at y={target_y}") + + trac_pos_tag = y_face(global_ymin) + trac_neg_tag = y_face(global_ymax) + + zero_disp_tags = [trac_pos_tag] + traction_tags = [trac_neg_tag] + + gmsh.option.setNumber("Mesh.Algorithm3D", 10) + gmsh.option.setNumber("Mesh.Optimize", 1) + gmsh.option.setNumber("Mesh.OptimizeNetgen", 1) + + gmsh.model.mesh.generate(3) + + gmsh.model.mesh.optimize("Netgen") + + gmsh.model.mesh.removeDuplicateNodes() + + msh_file = tempdir.name + "/mesh_merged.msh" + + gmsh.write(msh_file) + + gmsh.finalize() + + return msh_file, zero_disp_tags, traction_tags diff --git a/src/autobolt/parametric_cad_solver.py b/src/autobolt/parametric_cad_solver.py index 8c6e176..64203dc 100644 --- a/src/autobolt/parametric_cad_solver.py +++ b/src/autobolt/parametric_cad_solver.py @@ -1,7 +1,7 @@ import build123d -def create_two_plate_assembly( +def create_two_plate_assembly_components( plate_length_m: float, plate_width_m: float, plate_thickness_m: float, @@ -87,5 +87,31 @@ def create_two_plate_assembly( bolt_cylinders = bp2.part # 4) Assemble and return + return plateA, plateB, bolt_cylinders + + +def create_two_plate_assembly( + plate_length_m: float, + plate_width_m: float, + plate_thickness_m: float, + num_holes: int, + hole_radius_m: float, + edge_margin_m: float, + hole_spacing_m: float, + hole_offset_from_bottom_m: float, + plate_gap_mm: float, +): + plateA, plateB, bolt_cylinders = create_two_plate_assembly_components( + plate_length_m, + plate_width_m, + plate_thickness_m, + num_holes, + hole_radius_m, + edge_margin_m, + hole_spacing_m, + hole_offset_from_bottom_m, + plate_gap_mm, + ) + assembly = build123d.Compound(children=[plateA, plateB, bolt_cylinders]) - return assembly + return assembly \ No newline at end of file diff --git a/src/autobolt/testing.ipynb b/src/autobolt/testing.ipynb new file mode 100644 index 0000000..2cade50 --- /dev/null +++ b/src/autobolt/testing.ipynb @@ -0,0 +1,881 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "1b9d30d5", + "metadata": {}, + "outputs": [], + "source": [ + "import build123d\n", + "import fenics\n", + "import gmsh\n", + "import meshio\n", + "\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ddb3ce0c", + "metadata": {}, + "outputs": [], + "source": [ + "joint_configuration = {\n", + " \"load\": 60000,\n", + " \"desired_safety_factor\": 3.0,\n", + " \"bolt_yield_strength\": 940,\n", + " \"plate_yield_strength\": 250,\n", + " \"preload\": 150000,\n", + " \"pitch\": 1.5,\n", + " \"plate_thickness\": 10,\n", + " \"bolt_elastic_modulus\": 210,\n", + " \"plate_elastic_modulus\": 210,\n", + "\n", + " \"num_bolts\": 4,\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69c9d36c", + "metadata": {}, + "outputs": [], + "source": [ + "plate_width = 0.1 # [m]\n", + "plate_length = 0.2 # [m]\n", + "\n", + "plate_thickness = 10 # [mm]\n", + "plate_elastic_modulus = 210 # [GPa]\n", + "plate_yield_strength = 250 # [MPa]\n", + "\n", + "load = 60000 # [N]\n", + "traction = -load / (plate_thickness / 1000 * plate_length)\n", + "\n", + "num_bolts = 4\n", + "bolt_diameter = 10 # [mm]\n", + "bolt_elastic_modulus = 210 # [GPa]\n", + "bolt_yield_strength = 940 # [MPa]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14d432ae", + "metadata": {}, + "outputs": [], + "source": [ + "plate_thickness_m=plate_thickness / 1000\n", + "num_holes=num_bolts\n", + "elastic_modulus=plate_elastic_modulus * 10**9\n", + "yield_strength=plate_yield_strength * 10**6\n", + "traction_values=[(0, traction, 0)]\n", + "hole_radius_m=bolt_diameter / 2 / 1000\n", + "plate_length_m=plate_length\n", + "plate_width_m=plate_width\n", + "edge_margin_m=plate_length / (2 * num_bolts)\n", + "hole_spacing_m=plate_length / num_bolts\n", + "hole_offset_from_bottom_m=0.020 # [m] vertical position of hole centers (Y from bottom edge)\n", + "plate_gap_mm=0.01 # [mm] gap between the two plates\n", + "poissons_ratio=0.3 # Poisson's ratio for steel" + ] + }, + { + "cell_type": "markdown", + "id": "bde1b7b8", + "metadata": {}, + "source": [ + "# Test code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1592b99f", + "metadata": {}, + "outputs": [], + "source": [ + "import fea_solver_dual\n", + "from parametric_cad_solver import create_two_plate_assembly_components" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5f45e08", + "metadata": {}, + "outputs": [], + "source": [ + "plateA,plateB,boltCylinders = create_two_plate_assembly_components(\n", + " plate_length_m=plate_length_m,\n", + " plate_width_m=plate_width_m,\n", + " plate_thickness_m=plate_thickness_m,\n", + " num_holes=num_holes,\n", + " hole_radius_m=hole_radius_m,\n", + " edge_margin_m=edge_margin_m,\n", + " hole_spacing_m=hole_spacing_m,\n", + " hole_offset_from_bottom_m=hole_offset_from_bottom_m,\n", + " plate_gap_mm=plate_gap_mm\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db13f7f3", + "metadata": {}, + "outputs": [], + "source": [ + "fea_solver_dual._calculate_fos_from_build123d_dual(\n", + " plateA,\n", + " plateB,\n", + " boltCylinders,\n", + " elastic_modulus,\n", + " poissons_ratio,\n", + " bolt_yield_strength * 10**6,\n", + " plate_yield_strength * 10**6,\n", + " traction_values,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "dca3b813", + "metadata": {}, + "source": [ + "# Test old code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46f65838", + "metadata": {}, + "outputs": [], + "source": [ + "import fea_solver\n", + "from parametric_cad_solver import create_two_plate_assembly" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0f2a629", + "metadata": {}, + "outputs": [], + "source": [ + "object = create_two_plate_assembly(\n", + " plate_length_m=plate_length_m,\n", + " plate_width_m=plate_width_m,\n", + " plate_thickness_m=plate_thickness_m,\n", + " num_holes=num_holes,\n", + " hole_radius_m=hole_radius_m,\n", + " edge_margin_m=edge_margin_m,\n", + " hole_spacing_m=hole_spacing_m,\n", + " hole_offset_from_bottom_m=hole_offset_from_bottom_m,\n", + " plate_gap_mm=plate_gap_mm\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b318b5cb", + "metadata": {}, + "outputs": [], + "source": [ + "fea_solver._calculate_fos_from_build123d(\n", + " object,\n", + " elastic_modulus,\n", + " poissons_ratio,\n", + " yield_strength,\n", + " traction_values\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "57d88764", + "metadata": {}, + "source": [ + "# Create assembly step file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38529e0f", + "metadata": {}, + "outputs": [], + "source": [ + "# convert the gap to meters\n", + "gap_m = plate_gap_mm * 1e-4\n", + "\n", + "# 1) Build Plate A (with through‑holes)\n", + "with build123d.BuildPart() as bp:\n", + " build123d.Box(plate_length_m, plate_width_m, plate_thickness_m)\n", + " # center the box around the XY origin\n", + " bp.part = bp.part.move(build123d.Pos(plate_length_m / 2, plate_width_m / 2, 0))\n", + "\n", + " # compute X‑coordinates of hole centers\n", + " x0 = plate_length_m - edge_margin_m\n", + " x_coords = [x0 - i * hole_spacing_m for i in range(num_holes)]\n", + "\n", + " # subtract each hole\n", + " for x in x_coords:\n", + " with build123d.Locations(build123d.Pos(x, hole_offset_from_bottom_m, 0)):\n", + " build123d.Cylinder(\n", + " radius=hole_radius_m,\n", + " height=plate_thickness_m + 1e-3, # +1 mm to ensure clean cut\n", + " mode=build123d.Mode.SUBTRACT,\n", + " )\n", + "plateA = bp.part\n", + "\n", + "# 2) Copy Plate A → shift up by thickness+gap → mirror in Y → re‑position\n", + "plateB = build123d.copy(plateA)\n", + "# shift in Z by plate_thickness + gap\n", + "plateB = plateB.locate(\n", + " build123d.Location((0, 0, plate_thickness_m + gap_m), (0, 0, 1), 0)\n", + ")\n", + "# mirror about the XZ plane (so y→−y)\n", + "plateB = build123d.mirror(plateB, about=build123d.Plane.XZ)\n", + "# then translate back up in Y by 2*hole_center_y to line holes up\n", + "plateB = plateB.locate(\n", + " build123d.Location((0, 2 * hole_offset_from_bottom_m, 0), (0, 0, 1), 0)\n", + ")\n", + "\n", + "# 3) Create the bolt‑like cylinders (conformal, centered in the holes)\n", + "bolt_height_m = 2 * plate_thickness_m + gap_m\n", + "z_offset_m = plate_thickness_m / 2 + gap_m / 2 # start mid‑thickness of Plate A\n", + "\n", + "with build123d.BuildPart() as bp2:\n", + " for x in x_coords:\n", + " with build123d.Locations(\n", + " build123d.Pos(x, hole_offset_from_bottom_m, z_offset_m)\n", + " ):\n", + " build123d.Cylinder(\n", + " radius=hole_radius_m,\n", + " height=bolt_height_m,\n", + " align=(\n", + " build123d.Align.CENTER,\n", + " build123d.Align.CENTER,\n", + " build123d.Align.CENTER,\n", + " ),\n", + " )\n", + "bolt_cylinders = bp2.part\n", + "\n", + "# 4) Assemble and return\n", + "assembly = build123d.Compound(children=[plateA, plateB, bolt_cylinders])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dea990d2", + "metadata": {}, + "outputs": [], + "source": [ + "assembly" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# convert build123d mesh to gsmh mesh to fenics mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f286dc3", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fcffe0c1", + "metadata": {}, + "outputs": [], + "source": [ + "gmsh.initialize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69fece4a", + "metadata": {}, + "outputs": [], + "source": [ + "gmsh.clear()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2c637538", + "metadata": {}, + "outputs": [], + "source": [ + "tempdir = \"temp\"\n", + "\n", + "step_file = tempdir + \"/model.step\"\n", + "build123d.export_step(assembly, step_file)\n", + "\n", + "msh_file = tempdir + \"/mesh_merged.msh\"\n", + "xdmf_file_mesh = tempdir + \"/mesh.xdmf\"\n", + "xdmf_file_surface_tags = tempdir + \"surface_tags.xdmf\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4dcd012a", + "metadata": {}, + "outputs": [], + "source": [ + "def import_step_and_capture_new_volumes(step_path: str) -> list[int]:\n", + " \"\"\"Import a STEP file and return the list of volume entity tags created by that import.\"\"\"\n", + " before = set(t for _, t in gmsh.model.getEntities(3))\n", + " gmsh.merge(step_path)\n", + " gmsh.model.occ.synchronize()\n", + " after = set(t for _, t in gmsh.model.getEntities(3))\n", + " new_tags = sorted(after - before)\n", + " return new_tags" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c8b4b5f", + "metadata": {}, + "outputs": [], + "source": [ + "tempdir = tempfile.TemporaryDirectory()\n", + "\n", + "# --- Export STEP files (no geometry parameter changes) ---\n", + "step_plateA = tempdir.name + \"/plateA.step\"\n", + "step_plateB = tempdir.name + \"/plateB.step\"\n", + "step_bolts = tempdir.name + \"/bolts.step\"\n", + "build123d.export_step(plateA, step_plateA)\n", + "build123d.export_step(plateB, step_plateB)\n", + "build123d.export_step(bolt_cylinders, step_bolts)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f75d5900", + "metadata": {}, + "outputs": [], + "source": [ + "plateA_vols = import_step_and_capture_new_volumes(step_plateA)\n", + "plateB_vols = import_step_and_capture_new_volumes(step_plateB)\n", + "bolt_vols = import_step_and_capture_new_volumes(step_bolts)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db8815a4", + "metadata": {}, + "outputs": [], + "source": [ + "in_dimtags = [(3, t) for t in (plateA_vols + plateB_vols + bolt_vols)]\n", + "_, out_map = gmsh.model.occ.fragment(in_dimtags, [])\n", + "gmsh.model.occ.synchronize()\n", + "\n", + "gmsh.model.occ.removeAllDuplicates()\n", + "gmsh.model.occ.synchronize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc4e0b7c", + "metadata": {}, + "outputs": [], + "source": [ + "bolt_set_pre = set(bolt_vols)\n", + "\n", + "bolt_frag_vols: set[int] = set()\n", + "plate_frag_vols: set[int] = set()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3cff1ba", + "metadata": {}, + "outputs": [], + "source": [ + "for (src_dim, src_tag), frags in zip(in_dimtags, out_map):\n", + " frag_vols = [t for dim, t in frags if dim == 3]\n", + " if src_tag in bolt_set_pre:\n", + " bolt_frag_vols.update(frag_vols)\n", + " else:\n", + " plate_frag_vols.update(frag_vols)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "292fcd0a", + "metadata": {}, + "outputs": [], + "source": [ + "BOLT_TAG = 101\n", + "PLATE_TAG = 102\n", + "gmsh.model.addPhysicalGroup(3, sorted(bolt_frag_vols), tag=BOLT_TAG)\n", + "gmsh.model.setPhysicalName(3, BOLT_TAG, \"bolt\")\n", + "gmsh.model.addPhysicalGroup(3, sorted(plate_frag_vols), tag=PLATE_TAG)\n", + "gmsh.model.setPhysicalName(3, PLATE_TAG, \"plate\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3d636436", + "metadata": {}, + "outputs": [], + "source": [ + "surfaces = gmsh.model.getEntities(2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8da5cfe3", + "metadata": {}, + "outputs": [], + "source": [ + "for (sdim, stag) in surfaces:\n", + " gmsh.model.addPhysicalGroup(sdim, [stag], tag=stag)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bada8290", + "metadata": {}, + "outputs": [], + "source": [ + "surface_info = []\n", + "for dim, s in gmsh.model.getEntities(2):\n", + " xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.occ.getBoundingBox(dim, s)\n", + " surface_info.append((s, ymin, ymax))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a93b4f8", + "metadata": {}, + "outputs": [], + "source": [ + "global_ymin = min(ymin for _, ymin, _ in surface_info)\n", + "global_ymax = max(ymax for _, _, ymax in surface_info)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e76a49db", + "metadata": {}, + "outputs": [], + "source": [ + "def y_face(target_y: float, tol: float = 1e-6) -> int:\n", + " for tag, ymin, ymax in surface_info:\n", + " if abs(ymax - ymin) < tol and abs(ymin - target_y) < tol:\n", + " return tag\n", + " raise RuntimeError(f\"Could not find planar Y-face at y={target_y}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86d891d4", + "metadata": {}, + "outputs": [], + "source": [ + "trac_pos_tag = y_face(global_ymin)\n", + "trac_neg_tag = y_face(global_ymax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74bf4eed", + "metadata": {}, + "outputs": [], + "source": [ + "zero_disp_tags = [trac_pos_tag]\n", + "traction_tags = [trac_neg_tag]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46836935", + "metadata": {}, + "outputs": [], + "source": [ + "gmsh.option.setNumber(\"Mesh.Algorithm3D\", 10) # or 4 if 10 not supported\n", + "gmsh.option.setNumber(\"Mesh.Optimize\", 1)\n", + "gmsh.option.setNumber(\"Mesh.OptimizeNetgen\", 1)\n", + "\n", + "gmsh.model.mesh.generate(3)\n", + "\n", + "gmsh.model.mesh.optimize(\"Netgen\")\n", + "\n", + "gmsh.model.mesh.removeDuplicateNodes()\n", + "gmsh.write(msh_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd688185", + "metadata": {}, + "outputs": [], + "source": [ + "gmsh.finalize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc3be86f", + "metadata": {}, + "outputs": [], + "source": [ + "msh = meshio.read(msh_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15a46b33", + "metadata": {}, + "outputs": [], + "source": [ + "def to_xdmf(mesh, cell_type: str, out_path: str):\n", + " cells = mesh.get_cells_type(cell_type)\n", + " cell_data = mesh.get_cell_data(\"gmsh:physical\", cell_type)\n", + " out = meshio.Mesh(\n", + " points=mesh.points,\n", + " cells={cell_type: cells},\n", + " cell_data={\"name_to_read\": [cell_data]},\n", + " )\n", + " meshio.write(out_path, out)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c4feb0ce", + "metadata": {}, + "outputs": [], + "source": [ + "to_xdmf(msh, \"tetra\", xdmf_file_mesh)\n", + "to_xdmf(msh, \"triangle\", xdmf_file_surface_tags)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8d7c8ba", + "metadata": {}, + "outputs": [], + "source": [ + "mesh = fenics.cpp.mesh.Mesh()\n", + "with fenics.cpp.io.XDMFFile(xdmf_file_mesh) as infile:\n", + " infile.read(mesh)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a06da661", + "metadata": {}, + "outputs": [], + "source": [ + "mvc_surf = fenics.MeshValueCollection(\"size_t\", mesh, 2)\n", + "with fenics.cpp.io.XDMFFile(xdmf_file_surface_tags) as infile:\n", + " infile.read(mvc_surf, \"name_to_read\")\n", + "mf = fenics.cpp.mesh.MeshFunctionSizet(mesh, mvc_surf)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2a655b2", + "metadata": {}, + "outputs": [], + "source": [ + "mvc_cell = fenics.MeshValueCollection(\"size_t\", mesh, 3)\n", + "with fenics.XDMFFile(xdmf_file_mesh) as infile:\n", + " infile.read(mvc_cell, \"name_to_read\")\n", + "cf = fenics.MeshFunction(\"size_t\", mesh, mvc_cell)" + ] + }, + { + "cell_type": "markdown", + "id": "85895ccd", + "metadata": {}, + "source": [ + "# FOS calculation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d6d7346c", + "metadata": {}, + "outputs": [], + "source": [ + "V = fenics.VectorFunctionSpace(mesh, \"CG\", 2)\n", + "\n", + "bcs = [\n", + " fenics.DirichletBC(V, fenics.Constant((0.0, 0.0, 0.0)), mf, tag)\n", + " for tag in zero_disp_tags\n", + "]\n", + "\n", + "mu = elastic_modulus / (2.0 * (1.0 + poissons_ratio))\n", + "lmbda = elastic_modulus * poissons_ratio / (\n", + " (1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84052722", + "metadata": {}, + "outputs": [], + "source": [ + "def epsilon(u):\n", + " return 0.5 * (fenics.grad(u) + fenics.grad(u).T)\n", + "\n", + "def sigma(u):\n", + " return 2.0 * mu * epsilon(u) + lmbda * fenics.tr(epsilon(u)) * fenics.Identity(3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3261aa8", + "metadata": {}, + "outputs": [], + "source": [ + "u = fenics.TrialFunction(V)\n", + "v = fenics.TestFunction(V)\n", + "f = fenics.Constant((0.0, 0.0, 0.0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15d978dc", + "metadata": {}, + "outputs": [], + "source": [ + "a = fenics.inner(sigma(u), epsilon(v)) * fenics.dx" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17f78695", + "metadata": {}, + "outputs": [], + "source": [ + "ds = fenics.Measure(\"ds\", domain=mesh, subdomain_data=mf)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "befe86ae", + "metadata": {}, + "outputs": [], + "source": [ + "L = fenics.dot(f, v) * fenics.dx" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c49657e8", + "metadata": {}, + "outputs": [], + "source": [ + "for tag, traction in zip(traction_tags, traction_values):\n", + " L += fenics.dot(fenics.Constant(traction), v) * ds(tag)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5f17a26", + "metadata": {}, + "outputs": [], + "source": [ + "u_sol = fenics.Function(V)\n", + "fenics.solve(a == L, u_sol, bcs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c40e101", + "metadata": {}, + "outputs": [], + "source": [ + "s_dev = sigma(u_sol) - (1.0 / 3.0) * fenics.tr(sigma(u_sol)) * fenics.Identity(3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c11b29c0", + "metadata": {}, + "outputs": [], + "source": [ + "W0 = fenics.FunctionSpace(mesh, \"DG\", 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "643a4db0", + "metadata": {}, + "outputs": [], + "source": [ + "vm0 = fenics.project(fenics.sqrt(3.0 / 2.0 * fenics.inner(s_dev, s_dev)), W0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d297b3f5", + "metadata": {}, + "outputs": [], + "source": [ + "vals = vm0.vector().get_local()\n", + "markers = cf.array()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7653de2c", + "metadata": {}, + "outputs": [], + "source": [ + "bolt_vals = vals[markers == BOLT_TAG]\n", + "plate_vals = vals[markers == PLATE_TAG]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d02a5dd", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3b66326", + "metadata": {}, + "outputs": [], + "source": [ + "vm0.vector().max()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f349acb1", + "metadata": {}, + "outputs": [], + "source": [ + "max_bolt_vm = float(np.max(bolt_vals))\n", + "max_plate_vm = float(np.max(plate_vals))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "927a32b1", + "metadata": {}, + "outputs": [], + "source": [ + "max_plate_vm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d341c67f", + "metadata": {}, + "outputs": [], + "source": [ + "max_bolt_vm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04417bee", + "metadata": {}, + "outputs": [], + "source": [ + "bolt_fos = float(bolt_yield_strength / max_bolt_vm)\n", + "plate_fos = float(plate_yield_strength / max_plate_vm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45c08753", + "metadata": {}, + "outputs": [], + "source": [ + "print(bolt_fos, plate_fos)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "autobolt", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}