diff --git a/tools/python_utils/circle_wall_input_using_symmetric_mesh.py b/tools/python_utils/circle_wall_input_using_symmetric_mesh.py new file mode 100644 index 000000000..ca3742e45 --- /dev/null +++ b/tools/python_utils/circle_wall_input_using_symmetric_mesh.py @@ -0,0 +1,273 @@ +import os +import sys +import numpy as np + +#sys.path.append('/tools/python_utils') +from gmsh_particles import * +from util import * + +def create_input_file(inp_dir, pp_tag): + """Generates input file for two-particle test""" + + sim_inp_dir = str(inp_dir) + + R1 = 0.001 + mesh_size = R1 / 8. + horizon = 2.2 * mesh_size + particle_dist = 0.9*mesh_size #horizon # surface to surface distance + + ## particle 1 rectangle + Lx, Ly = 4*R1, 3*mesh_size + rect_center = [R1, 0., 0.] + + ## particle 2 circle + cir_center = [R1, rect_center[1] + 0.5*Ly + particle_dist + R1, 0.] + + ## assign free fall velocity to second particle + high_impact = False + free_fall_vel = [0., -0.1, 0.] + if high_impact: + free_fall_vel = [0., -4., 0.] # high impact velocity + + ## time + final_time = 0.01 + num_steps = 50000 + if high_impact: + final_time = 0.0002 + num_steps = 10000 + # final_time = 0.00002 + # num_steps = 2 + num_outputs = 100 + dt_out_n = num_steps / num_outputs + perform_out = True + + ## material + poisson1 = 0.25 + rho1 = 1200. + K1 = 2.16e+7 + E1 = get_E(K1, poisson1) + G1 = get_G(E1, poisson1) + Gc1 = 50. + + poisson2 = 0.25 + rho2 = 1200. + K2 = 2.16e+7 + E2 = get_E(K2, poisson2) + G2 = get_G(E2, poisson2) + Gc2 = 50. + + ## contact + # R_contact = 0.95 * mesh_size + # R_contact = 1.74e-04 + R_contact_factor = 0.95 + + # Kn_V_max = 7.385158e+05 + # Kn = np.power(Kn_V_max, 2) + # compute from bulk modulus + + # from bulk modulus + Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5)) + Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5)) + Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5)) + + beta_n_eps = 0.9 + friction_coeff = 0.5 + damping_active = False + friction_active = False + beta_n_factor = 10. + Kn_factor = 0.1 + + ## gravity + gravity_active = True + gravity = [0., -10., 0.] + + ## neighbor search details + neigh_search_factor = 2. + neigh_search_interval = 1 + neigh_search_criteria = "simple_all" + + ### ---------------------------------------------------------------- ### + # generate mesh and particle location data + ### ---------------------------------------------------------------- ### + plocf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w') + plocf.write("i, x, y, z, r, o\n") + plocf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (0, rect_center[0], rect_center[1], rect_center[2], Lx, 0.)) + plocf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (1, cir_center[0], cir_center[1], cir_center[2], R1, 0.)) + plocf.close() + + zones_mesh_fnames = ["mesh_rect_1", "mesh_cir_2"] + + # generate mesh for particle 1 + rectangle_mesh_symmetric(xc = [0., 0., 0.], Lx = Lx, Ly = Ly, h = mesh_size, filename = zones_mesh_fnames[0] + "_" + str(pp_tag), vtk_out = True, symmetric_mesh = True) + + # generate mesh for particle 2 + circle_mesh_symmetric(xc = [0., 0., 0.], r = R1, h = mesh_size, filename = zones_mesh_fnames[1] + "_" + str(pp_tag), vtk_out = True, symmetric_mesh = True) + + os.system("mkdir -p ../out") + + ### ---------------------------------------------------------------- ### + # generate YAML file + ### ---------------------------------------------------------------- ### + # print('\nGenerating imput file\n') + inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w') + inpf.write("Model:\n") + inpf.write(" Dimension: 2\n") + inpf.write(" Discretization_Type:\n") + inpf.write(" Spatial: finite_difference\n") + inpf.write(" Time: central_difference\n") + inpf.write(" Final_Time: %4.6e\n" % (final_time)) + inpf.write(" Time_Steps: %d\n" % (num_steps)) + + # + # container info + # + inpf.write("Container:\n") + inpf.write(" Geometry:\n") + inpf.write(" Type: rectangle\n") + contain_params = [rect_center[0] - 0.5*Lx, rect_center[1] - 0.5*Ly, 0., rect_center[0] + 0.5*Lx, rect_center[1] + 0.5*Ly + particle_dist + 2*R1, 0.] + inpf.write(" Parameters: " + print_dbl_list(contain_params)) + + # + # zone info + # + inpf.write("Zone:\n") + inpf.write(" Zones: 2\n") + + ## zone 1 (bottom particle) + inpf.write(" Zone_1:\n") + inpf.write(" Is_Wall: false\n") + + ## zone 2 (top particle) + inpf.write(" Zone_2:\n") + inpf.write(" Is_Wall: false\n") + + # + # particle info + # + inpf.write("Particle:\n") + inpf.write(" Zone_1:\n") + inpf.write(" Type: rectangle\n") + inpf.write(" Parameters: " + print_dbl_list([Lx, Ly, rect_center[0], rect_center[1], rect_center[2]])) + inpf.write(" Zone_2:\n") + inpf.write(" Type: circle\n") + p2_geom = [R1, cir_center[0], cir_center[1], cir_center[2]] + inpf.write(" Parameters: " + print_dbl_list(p2_geom)) + + # + # particle generation + # + inpf.write("Particle_Generation:\n") + inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n") + inpf.write(" File_Data_Type: loc_rad_orient\n") + + # + # Mesh info + # + inpf.write("Mesh:\n") + + for i in range(len(zones_mesh_fnames)): + inpf.write(" Zone_%d:\n" % (i+1)) + inpf.write(" File: %s\n" % (zones_mesh_fnames[i] + "_" + str(pp_tag) + ".msh")) + + # Contact info + inpf.write("Contact:\n") + + ## 11 + write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "11", Kn_11) + + ## copy from 11 + copy_contact_zone(inpf, [12, 22], [1, 1]) + + # Neighbor info + inpf.write("Neighbor:\n") + inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria)) + inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor)) + inpf.write(" Search_Interval: %d\n" % (neigh_search_interval)) + + # Material info + inpf.write("Material:\n") + + ## zone 1 + write_material_zone_part(inpf, "1", horizon, rho1, K1, G1, Gc1) + + ## zone 2 + inpf.write(" Zone_2:\n") + inpf.write(" Copy_Material_Data: 1\n") + + # + # Force + # + if gravity_active == True: + inpf.write("Force_BC:\n") + inpf.write(" Gravity: " + print_dbl_list(gravity)) + + # + # IC + # + inpf.write("IC:\n") + inpf.write(" Constant_Velocity:\n") + inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel)) + inpf.write(" Particle_List: [1]\n") + + # + # Displacement + # + inpf.write("Displacement_BC:\n") + inpf.write(" Sets: 1\n") + + inpf.write(" Set_1:\n") + inpf.write(" Particle_List: [0]\n") + inpf.write(" Direction: [1,2]\n") + inpf.write(" Time_Function:\n") + inpf.write(" Type: constant\n") + inpf.write(" Parameters:\n") + inpf.write(" - 0.0\n") + inpf.write(" Spatial_Function:\n") + inpf.write(" Type: constant\n") + inpf.write(" Zero_Displacement: true\n") + + # + # Output info + # + inpf.write("Output:\n") + inpf.write(" Path: ../out/\n") + inpf.write(" Tags:\n") + inpf.write(" - Displacement\n") + inpf.write(" - Velocity\n") + inpf.write(" - Force\n") + inpf.write(" - Force_Density\n") + inpf.write(" - Damage_Z\n") + inpf.write(" - Damage\n") + inpf.write(" - Nodal_Volume\n") + inpf.write(" - Zone_ID\n") + inpf.write(" - Particle_ID\n") + inpf.write(" - Fixity\n") + inpf.write(" - Force_Fixity\n") + inpf.write(" - Contact_Nodes\n") + inpf.write(" - No_Fail_Node\n") + inpf.write(" - Boundary_Node_Flag\n") + inpf.write(" - Theta\n") + inpf.write(" Output_Interval: %d\n" % (dt_out_n)) + inpf.write(" Compress_Type: zlib\n") + inpf.write(" Perform_FE_Out: false\n") + if perform_out: + inpf.write(" Perform_Out: true\n") + else: + inpf.write(" Perform_Out: false\n") + inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n)) + + inpf.write(" Debug: 3\n") + inpf.write(" Tag_PP: %d\n" %(int(pp_tag))) + + # close file + inpf.close() + +##-------------------------------------------------------## +##-------------------------------------------------------## +if __name__ == "__main__": + inp_dir = './' + pp_tag = 0 + if len(sys.argv) > 1: + pp_tag = int(sys.argv[1]) + + create_input_file(inp_dir, pp_tag) \ No newline at end of file diff --git a/tools/python_utils/geom_util.py b/tools/python_utils/geom_util.py new file mode 100644 index 000000000..e0426114f --- /dev/null +++ b/tools/python_utils/geom_util.py @@ -0,0 +1,381 @@ +import os +import sys +import numpy as np +import gmsh +import math + +def rotate(p, theta, axis): + """ + Returns rotation of vector about specified axis by specified angle + + Parameters + ---------- + p: list + Coordinates of vector + theta : float + Angle of rotation + axis : list + Axis of rotation + + Returns + ------- + list + Coordinates of rotated vector + """ + p_np, axis_np = np.array(p), np.array(axis) + ct, st = np.cos(theta), np.sin(theta) + p_dot_n = np.dot(p_np,axis_np) + n_cross_p = np.cross(axis_np, p_np) + + return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p + +def get_ref_rect_points(center, Lx, Ly, add_center = False): + """ + Returns size points on reference rectangle + + Reference rectangle: + + +-------------------+ + | | + | | + | | + +-------------------+ + + Parameters + ---------- + center: list + Coordinates of center of rectangle + Lx : float + Length of rectangle in x-direction + Ly : float + Length of rectangle in y-direction + add_center : bool + True if we include the center to the returned list (first point in the returned list will be the center if the flag is true) + + Returns + ------- + list + Coordinates of points + """ + + p1 = [center[0] - Lx/2., center[1] - Ly/2., center[2]] + p2 = [center[0] + Lx/2., center[1] - Ly/2., center[2]] + p3 = [center[0] + Lx/2., center[1] + Ly/2., center[2]] + p4 = [center[0] - Lx/2., center[1] + Ly/2., center[2]] + + return [p1, p2, p3, p4] if not add_center else [center, p1, p2, p3, p4] + +def get_ref_triangle_points(center, radius, add_center = False): + """ + Returns size points on reference triangle + + Reference triangle: + + + v2 + | \ + | \ + | \ + | o + v1 + | / + | / + | / + + v3 + + Radius: o-v1 + Axis vector: o-v1 + Rotation axis: [0, 0, 1] + + Parameters + ---------- + center: list + Coordinates of center of triangle + radius : float + Radius of triangle + add_center : bool + True if we include the center to the returned list (first point in the returned list will be the center if the flag is true) + + Returns + ------- + list + Coordinates of points + """ + axis = [1., 0., 0.] + rotate_axis = [0., 0., 1.] + + points = [] + if add_center: + points.append(center) + + for i in range(3): + xi = rotate(axis, i*2*np.pi/3., rotate_axis) + points.append([center[i] + radius * xi[i] for i in range(3)]) + + return points + +def get_ref_hex_points(center, radius, add_center = False): + """ + Returns size points on reference hexagon + + Reference hexagon: + + v3 v2 + + + + + + + o + + v4 x v1 + + + + + v5 v6 + + Radius: x-v1 + Axis vector: x-v1 + Rotation axis: [0, 0, 1] + + Parameters + ---------- + center: list + Coordinates of center of hexagon + radius : float + Radius of hexagon + add_center : bool + True if we include the center to the returned list (first point in the returned list will be the center if the flag is true) + + Returns + ------- + list + Coordinates of points + """ + axis = [1., 0., 0.] + rotate_axis = [0., 0., 1.] + points = [] + if add_center: + points.append(center) + + for i in range(6): + xi = rotate(axis, i*np.pi/3., rotate_axis) + points.append([center[i] + radius * xi[i] for i in range(3)]) + + return points + +def get_ref_drum_points(center, radius, width, add_center = False): + """ + Returns size points on reference concave polygon (we refer to it as drum) + + Reference concave polygon: + + v3 v2 + + -------------------------------- + + \ / + \ / + + o + + /v4 x v1 \ + / \ + + -------------------------------- + + v5 v6 + + Radius: x-v2 + Width: v1-v4 + Axis vector: x-v1 + Rotation axis: [0, 0, 1] + + Parameters + ---------- + center: list + Coordinates of center of polygon + radius : float + Radius of polygon (distance between center x and vertex v2) + width : float + Width of neck, i.e. distance between vertex v2 and v4 + add_center : bool + True if we include the center to the returned list (first point in the returned list will be the center if the flag is true) + + Returns + ------- + list + Coordinates of points + """ + axis = [1., 0., 0.] + rotate_axis = [0., 0., 1.] + + x1 = rotate(axis, np.pi/3., rotate_axis) + x2 = rotate(axis, -np.pi/3., rotate_axis) + + points = [] + if add_center: + points.append(center) + + # v1 + points.append([center[i] + width*0.5*axis[i] for i in range(3)]) + # v2 + points.append([center[i] + radius*x1[i] for i in range(3)]) + # v3 + points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)]) + # v4 + points.append([center[i] - width*0.5*axis[i] for i in range(3)]) + # v5 + v6 = [center[i] + radius*x2[i] for i in range(3)] + points.append([v6[i] - radius*axis[i] for i in range(3)]) + # v6 + points.append(v6) + + return points + + +def does_intersect(p, particles, padding): + """ + Returns true if particle p intersects or is near enough to existing particles + + Parameters + ---------- + p : list + Coordinates of center and radius of particle [x,y,z,r] + particles : list + List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius. + padding: float + Minimum distance between circle boundaries such that if two circles + + Returns + ------- + bool + True if particle intersects or is near enough to one of the particle in the list + """ + if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p)) + if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles)) + if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding)) + + for q in particles: + if len(q) < 4: raise Exception('q = {} in particles must have atleast 4 elements'.format(q)) + pq = np.array([p[i] - q[i] for i in range(3)]) + if np.linalg.norm(pq) <= p[3] + q[3] + padding: + return True + return False + +def does_intersect_rect(p, particles, padding, rect, is_3d = False): + """ + Returns true if particle p is sufficiently close or outside the rectangle (in 2d) or cuboid (in 3d) + + Parameters + ---------- + p : list + Coordinates of center and radius of particle [x,y,z,r] + particles : list + List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius. + padding: float + Minimum distance between circle boundaries such that if two circles + rect: list + Coordinates of left-bottom and right-top corner points of rectangle (2d) or cuboid (3d). E.g. [x1 y1, z1, x2, y2, z2] + is_3d: bool + True if we are dealing with cuboid + + Returns + ------- + bool + True if particle intersects or is near enough to the rectangle + """ + if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p)) + if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles)) + if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding)) + if len(rect) < 6: raise Exception('rect = {} must have 6 elements'.format(rect)) + + pr = [p[0] - p[3], p[1] - p[3], p[2], p[0] + p[3], p[1] + p[3], p[2]] + if is_3d: + pr[2] -= p[3] + pr[5] += p[3] + + if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding: + if is_3d: + if pr[2] < rect[2] + padding or pr[5] > rect[5] - padding: + return True + else: + return True + return False + + +def get_gmsh_entities(): + m = {} + for e in gmsh.model.getEntities(): + bnd = gmsh.model.getBoundary([e]) + nod = gmsh.model.mesh.getNodes(e[0], e[1]) + ele = gmsh.model.mesh.getElements(e[0], e[1]) + m[e] = (bnd, nod, ele) + + return m + +# transform the mesh and create new discrete entities to store it +# source - https://gitlab.onelab.info/gmsh/gmsh/-/blob/master/examples/api/mirror_mesh.py +def gmsh_transform(m, offset_entity, offset_node, offset_element, tx, ty, tz): + for e in sorted(m): + gmsh.model.addDiscreteEntity( + e[0], e[1] + offset_entity, + [(abs(b[1]) + offset_entity) * int(math.copysign(1, b[1])) for b in m[e][0]]) + coord = [] + for i in range(0, len(m[e][1][1]), 3): + x = m[e][1][1][i] * tx + y = m[e][1][1][i + 1] * ty + z = m[e][1][1][i + 2] * tz + coord.append(x) + coord.append(y) + coord.append(z) + gmsh.model.mesh.addNodes(e[0], e[1] + offset_entity, + [n + offset_node for n in m[e][1][0]], coord) + gmsh.model.mesh.addElements(e[0], e[1] + offset_entity, m[e][2][0], + [[t + offset_element for t in typ] for typ in m[e][2][1]], + [[n + offset_node for n in typ] for typ in m[e][2][2]]) + if (tx * ty * tz) < 0: # reverse the orientation + gmsh.model.mesh.reverse([(e[0], e[1] + offset_entity)]) + +def gmsh_transform_general(m, offset_entity, offset_node, offset_element, angle, axis=[0, 0, 1]): + """Rotate mesh by angle (in radians) around specified axis + + Parameters: + m - mesh data dictionary + offset_entity - offset for entity numbers + offset_node - offset for node numbers + offset_element - offset for element numbers + angle - rotation angle in radians + axis - rotation axis vector (default is z-axis [0,0,1]) + """ + + for e in sorted(m): + gmsh.model.addDiscreteEntity( + e[0], e[1] + offset_entity, + [(abs(b[1]) + offset_entity) * int(math.copysign(1, b[1])) for b in m[e][0]]) + coord = [] + # Apply rotation to each point + for i in range(0, len(m[e][1][1]), 3): + xi = np.array([m[e][1][1][i], m[e][1][1][i + 1], m[e][1][1][i + 2]]) + + xi_rot = rotate(xi, angle, axis) + + coord.append(xi_rot[0]) + coord.append(xi_rot[1]) + coord.append(xi_rot[2]) + + gmsh.model.mesh.addNodes(e[0], e[1] + offset_entity, + [n + offset_node for n in m[e][1][0]], coord) + gmsh.model.mesh.addElements(e[0], e[1] + offset_entity, m[e][2][0], + [[t + offset_element for t in typ] for typ in m[e][2][1]], + [[n + offset_node for n in typ] for typ in m[e][2][2]]) + +def gmsh_translate(xc): + """Translate mesh by vector xc + + Parameters: + xc - translation vector [x, y, z] + """ + # Get the nodes + nodeTags, coord, paramCoord = gmsh.model.mesh.getNodes() + numNodes = len(nodeTags) + + # Translation vector + translation_vector = np.array([xc[0], xc[1], xc[2]]) + + # Update each node individually + for i in range(numNodes): + nodeTag = nodeTags[i] + old_coord = coord[3*i:3*(i+1)] + new_coord = old_coord + translation_vector + # Handle paramCoord properly - use empty list if paramCoord is None or empty + param = paramCoord[3*i:3*(i+1)] if paramCoord is not None and len(paramCoord) > 0 else [] + gmsh.model.mesh.setNode(nodeTag, new_coord, param) \ No newline at end of file diff --git a/tools/python_utils/gmsh_particles.py b/tools/python_utils/gmsh_particles.py new file mode 100644 index 000000000..7cc0ffe69 --- /dev/null +++ b/tools/python_utils/gmsh_particles.py @@ -0,0 +1,1450 @@ +import os +import sys +import numpy as np +import gmsh +import math + +from util import * +from geom_util import * + + +def circle_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a circular mesh, either symmetric or full. + xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh + r - radius of the circle + h - mesh size + symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full circle + """ + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # we first assume circle center is at origin and then translate the symmetric mesh to the desired location + xc_mesh = [0., 0., 0.] + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) + p2 = gmsh.model.geo.addPoint(xc_mesh[0] + r, xc_mesh[1], xc_mesh[2], h) + p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + r, xc_mesh[2], h) + + l1 = gmsh.model.geo.addCircleArc(p2, p1, p3) + l2 = gmsh.model.geo.addLine(p1, p2) + l3 = gmsh.model.geo.addLine(p3, p1) + + c1 = gmsh.model.geo.addCurveLoop([l2, l1, l3]) + + p1 = gmsh.model.geo.addPlaneSurface([c1]) + + #gmsh.model.occ.addBox(0,0,0, 1,0.5,0.5) + gmsh.model.geo.synchronize() + # gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.1) + # gmsh.model.mesh.setSize([(0, 2)], 0.01) + gmsh.model.mesh.generate(3) + # gmsh.write('mesh.vtk') + + # get the mesh data + m = get_gmsh_entities() + + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis + + # remove the duplicate nodes that will have been created on the internal + # boundaries + gmsh.model.mesh.removeDuplicateNodes() + + # translate the mesh to the specified center coordinates + gmsh_translate(xc) + + else: + # Create circle using built-in gmsh circle + c = gmsh.model.occ.addCircle(xc[0], xc[1], xc[2], r) + + # Create curve loop from circle + cl = gmsh.model.occ.addCurveLoop([c]) + + # Create surface from curve loop + s = gmsh.model.occ.addPlaneSurface([cl]) + + # Add center point + p = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2], h) + + # Synchronize + gmsh.model.occ.synchronize() + + # Embed center point in surface + gmsh.model.mesh.embed(0, [p], 2, s) + + # generate mesh + gmsh.model.mesh.generate(3) + + # write + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def ellipse_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create an elliptical mesh, either symmetric or full. + xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh + rx - radius in x direction + ry - radius in y direction + h - mesh size + symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full ellipse + """ + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # we first assume ellipse center is at origin and then translate the symmetric mesh to the desired location + xc_mesh = [0., 0., 0.] + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center + p2 = gmsh.model.geo.addPoint(xc_mesh[0] + rx, xc_mesh[1], xc_mesh[2], h) # Right + p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + ry, xc_mesh[2], h) # Top + + # Create elliptical arc using 3 points + l1 = gmsh.model.geo.addEllipseArc(p2, p1, p2, p3) # Right to top + l2 = gmsh.model.geo.addLine(p1, p2) # Center to right + l3 = gmsh.model.geo.addLine(p3, p1) # Top to center + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop([l2, l1, l3]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + # get the mesh data + m = get_gmsh_entities() + + # Mirror the mesh to create full ellipse + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis + + # Remove duplicate nodes on internal boundaries + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center coordinates + gmsh_translate(xc) + + else: + # Create points for full ellipse + p1 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h) # Center + p2 = gmsh.model.geo.addPoint(xc[0] + rx, xc[1], xc[2], h) # Right + p3 = gmsh.model.geo.addPoint(xc[0], xc[1] + ry, xc[2], h) # Top + p4 = gmsh.model.geo.addPoint(xc[0] - rx, xc[1], xc[2], h) # Left + p5 = gmsh.model.geo.addPoint(xc[0], xc[1] - ry, xc[2], h) # Bottom + + # Create elliptical arcs + l1 = gmsh.model.geo.addEllipseArc(p2, p1, p2, p3) # Right to top + l2 = gmsh.model.geo.addEllipseArc(p3, p1, p3, p4) # Top to left + l3 = gmsh.model.geo.addEllipseArc(p4, p1, p4, p5) # Left to bottom + l4 = gmsh.model.geo.addEllipseArc(p5, p1, p5, p2) # Bottom to right + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # synchronize + gmsh.model.geo.synchronize() + + # embed p0 (dim 0) in surface 1 (dim 2) + gmsh.model.mesh.embed(0, [p1], 2, s1) + + # generate mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def sphere_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a spherical mesh, either symmetric or full. + xc - center point coordinates [x, y, z] + r - radius of the sphere + h - mesh size + symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full sphere + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/8th of sphere at origin, then mirror it + xc_mesh = [0., 0., 0.] + + # Create first octant of sphere using OpenCASCADE + sphere = gmsh.model.occ.addSphere(xc_mesh[0], xc_mesh[1], xc_mesh[2], r) + + # cubes to cut the sphere + cut_left = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2]-r, r, 2*r, 2*r) + cut_back = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2]-r, 2*r, r, 2*r) + cut_bottom = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2]-r, 2*r, 2*r, r) + + # Cut sphere to get first octant using 3 planes + gmsh.model.occ.cut([(3,sphere)], [(3,cut_left), (3,cut_back), (3,cut_bottom)]) + + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # gmsh.fltk.run() + + # Get mesh data for mirroring + m = get_gmsh_entities() + + # Mirror mesh in x direction + gmsh_transform(m, 1000, 10000, 100000, -1, 1, 1) + # Mirror in y direction + gmsh_transform(m, 2000, 20000, 200000, 1, -1, 1) + # Mirror in z direction + gmsh_transform(m, 3000, 30000, 300000, 1, 1, -1) + # Mirror in xy plane + gmsh_transform(m, 4000, 40000, 400000, -1, -1, 1) + # Mirror in yz plane + gmsh_transform(m, 5000, 50000, 500000, 1, -1, -1) + # Mirror in xz plane + gmsh_transform(m, 6000, 60000, 600000, -1, 1, -1) + # Mirror in xyz + gmsh_transform(m, 7000, 70000, 700000, -1, -1, -1) + + # Translate to specified center if not at origin + gmsh_translate(xc) + + else: + # Create points for full sphere + # Create full sphere using OpenCASCADE + sphere = gmsh.model.occ.addSphere(xc[0], xc[1], xc[2], r) + + # Synchronize the OpenCASCADE geometry + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Synchronize + gmsh.model.geo.synchronize() + + # Embed center point in volume + p1 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h) + gmsh.model.mesh.embed(0, [p1], 3, sphere) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def rectangle_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a rectangular mesh, either symmetric or full. + xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh + Lx - length in x direction + Ly - length in y direction + h - mesh size + symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full rectangle + """ + # Initialize gmsh + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # center at origin + xc_mesh = [0., 0., 0.] + + # Create points for 1/4 rectangle + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) + p2 = gmsh.model.geo.addPoint(xc_mesh[0]+0.5*Lx, xc_mesh[1], xc_mesh[2], h) + p3 = gmsh.model.geo.addPoint(xc_mesh[0]+0.5*Lx, xc_mesh[1]+0.5*Ly, xc_mesh[2], h) + p4 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1]+0.5*Ly, xc_mesh[2], h) + + # Create lines + l1 = gmsh.model.geo.addLine(p1, p2) + l2 = gmsh.model.geo.addLine(p2, p3) + l3 = gmsh.model.geo.addLine(p3, p4) + l4 = gmsh.model.geo.addLine(p4, p1) + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # Synchronize and generate mesh + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + # Get the mesh data for mirroring + m = get_gmsh_entities() + + # Mirror the mesh in all quadrants + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # translate the mesh to the specified center coordinates + gmsh_translate(xc) + + else: + # Create points for full rectangle + p0 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h) # Center + p1 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom left + p2 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom right + p3 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top right + p4 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top left + + # Create lines + l1 = gmsh.model.geo.addLine(p1, p2) # Bottom + l2 = gmsh.model.geo.addLine(p2, p3) # Right + l3 = gmsh.model.geo.addLine(p3, p4) # Top + l4 = gmsh.model.geo.addLine(p4, p1) # Left + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # synchronize + gmsh.model.geo.synchronize() + + # embed p0 (dim 0) in surface 1 (dim 2) + gmsh.model.mesh.embed(0, [p0], 2, s1) + + # generate mesh + gmsh.model.mesh.generate(3) + + # write + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def hexagon_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a hexagon mesh, either symmetric or full. + xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh + r - radius of the hexagon + h - mesh size + symmetric_mesh - if True, creates 1/6 mesh and mirrors it. If False, creates full hexagon + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/6 hexagon mesh at origin first, then rotate and translate + xc_mesh = [0., 0., 0.] + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center (origin) + p2 = gmsh.model.geo.addPoint(xc_mesh[0] + r, xc_mesh[1], xc_mesh[2], h) # Point on x-axis (p1) + p3 = gmsh.model.geo.addPoint(xc_mesh[0] + r*np.cos(np.pi/3), xc_mesh[1] + r*np.sin(np.pi/3), xc_mesh[2], h) # p2 + + # Create lines for 1/6th of hexagon + l1 = gmsh.model.geo.addLine(p1, p2) # Center to p1 + l2 = gmsh.model.geo.addLine(p2, p3) # p1 to p2 + l3 = gmsh.model.geo.addLine(p3, p1) # p2 to center + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + # Get mesh data for rotation + m = get_gmsh_entities() + + # Rotate the mesh 5 times by 60 degrees + for i in range(5): + angle = (i + 1) * np.pi/3 + gmsh_transform_general(m, 1000*(i+1), 1000000*(i+1), 1000000*(i+1), angle) + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center + gmsh_translate(xc) + + else: + # For non-symmetric mesh, create full hexagon directly + points = [] + for i in range(6): + angle = i * np.pi/3 + px = xc[0] + r * np.cos(angle) + py = xc[1] + r * np.sin(angle) + points.append(gmsh.model.geo.addPoint(px, py, xc[2], h)) + + points.append(gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h)) + + # Create lines connecting points + lines = [] + for i in range(5): + lines.append(gmsh.model.geo.addLine(points[i], points[i+1])) + lines.append(gmsh.model.geo.addLine(points[5], points[0])) + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop(lines) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # synchronize + gmsh.model.geo.synchronize() + + # embed p1 (dim 0) in surface 1 (dim 2) + gmsh.model.mesh.embed(0, [points[-1]], 2, s1) + + # generate mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def drum2d_mesh_symmetric(xc = [0., 0., 0.], r = 1., width = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a drum2d mesh, either symmetric or full. + xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh + r - radius of the drum2d + width - width of the drum2d + h - mesh size + symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full drum2d + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/4 drum2d mesh at origin first, then rotate and translate + xc_mesh = [0., 0., 0.] + x_drum = get_ref_drum_points(xc_mesh, r, width) + v1, v2, v3 = x_drum[0], x_drum[1], x_drum[2] + v23 = [0.5*(v2[0] + v3[0]), 0.5*(v2[1] + v3[1]), 0.5*(v2[2] + v3[2])] + + # 1/4th mesh is a polygon with 4 points (o, v1, v2, (v2+v3)/2) + # Create points for 1/4th mesh + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center point + p2 = gmsh.model.geo.addPoint(v1[0], v1[1], v1[2], h) + p3 = gmsh.model.geo.addPoint(v2[0], v2[1], v2[2], h) + p4 = gmsh.model.geo.addPoint(v23[0], v23[1], v23[2], h) + + # Create lines + lines = [] + points = [p1, p2, p3, p4, p1] # Add p1 again at end to complete loop + for i in range(4): + lines.append(gmsh.model.geo.addLine(points[i], points[i+1])) + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop(lines) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + # Get mesh data for mirroring + m = get_gmsh_entities() + + # Mirror mesh to create full drum2d + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis + + # Remove duplicate nodes on internal boundaries + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center + gmsh_translate(xc) + + else: + # For non-symmetric mesh, create full drum2d directly + x_drum = get_ref_drum_points(xc, r, width) + # add points + points = [] + points.append(gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h)) + for i in range(len(x_drum)): + points.append(gmsh.model.geo.addPoint(x_drum[i][0], x_drum[i][1], x_drum[i][2], h)) + + # Create lines connecting vertices + lines = [] + points.append(points[1]) # to close the loop + for i in range(len(x_drum)): + lines.append(gmsh.model.geo.addLine(points[i+1], points[i+2])) + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop(lines) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # synchronize + gmsh.model.geo.synchronize() + + # embed p1 (dim 0) in surface 1 (dim 2) + gmsh.model.mesh.embed(0, [points[0]], 2, s1) + + # generate mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def polygon_mesh_symmetric(points, theta, xc=[0., 0., 0.], h=0.1, filename='mesh', vtk_out=False, symmetric_mesh = True): + """ + Create a symmetric mesh by rotating a polygon segment n times, where n = 360/theta. + + Parameters: + points - List of [x,y,z] coordinates defining the polygon segment. First point must be [0,0,0] + and first two edges must form angle theta at origin + theta - Angle in radians between first two edges at origin. Must divide 2pi evenly. + xc - Center point for final translation + h - Mesh size + filename - Output filename + vtk_out - Whether to output VTK file + """ + # Validate inputs + if not points or len(points) < 3: + raise ValueError("Need at least 3 points to define a polygon segment") + + if not np.allclose(points[0], [0., 0., 0.]): + raise ValueError("First point must be at origin [0,0,0]") + + # Check if theta divides 360 evenly + n = int(round(2*np.pi/theta)) + if not np.isclose(2*np.pi, n * theta): + raise ValueError(f"Angle {theta} degrees must divide 360 evenly") + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + # Create points for the polygon segment + gmsh_points = [] + for pt in points: + gmsh_points.append(gmsh.model.geo.addPoint(pt[0], pt[1], pt[2], h)) + + # Create lines connecting points + lines = [] + for i in range(len(gmsh_points)-1): + lines.append(gmsh.model.geo.addLine(gmsh_points[i], gmsh_points[i+1])) + # Close the polygon segment + lines.append(gmsh.model.geo.addLine(gmsh_points[-1], gmsh_points[0])) + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop(lines) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # Synchronize and generate initial mesh + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + # Get mesh data for rotation + m = get_gmsh_entities() + + # Rotate n-1 times by theta + for i in range(1, n): + angle = i * theta # Convert to radians + gmsh_transform_general(m, 1000*i, 1000000*i, 1000000*i, angle) + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to final position if needed + gmsh_translate(xc) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def cylindrical2d_wall_mesh(center=[0., 0., 0.], outer_radius=1.0, inner_radius=0.8, + bar_width=0.2, bar_length=0.3, h=0.1, filename='mesh', vtk_out=False): + """ + Create a cylindrical wall mesh with bars. + + Parameters: + center - Center point coordinates [x, y, z] + outer_radius - Outer radius of the wall + inner_radius - Inner radius of the wall + bar_width - Width of the bars + bar_length - Length of the bars + h - Mesh size + filename - Output filename + vtk_out - Whether to output VTK file + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + # Create points + # Center point + p1 = gmsh.model.geo.addPoint(center[0], center[1], center[2], h) + + # Outer circle points + p2 = gmsh.model.geo.addPoint(center[0] + outer_radius, center[1], center[2], h) + p3 = gmsh.model.geo.addPoint(center[0], center[1] + outer_radius, center[2], h) + p4 = gmsh.model.geo.addPoint(center[0] - outer_radius, center[1], center[2], h) + p5 = gmsh.model.geo.addPoint(center[0], center[1] - outer_radius, center[2], h) + + # Inner points with bar + p6 = gmsh.model.geo.addPoint(center[0] + inner_radius, center[1] + 0.5*bar_width, center[2], h) + p7 = gmsh.model.geo.addPoint(center[0] + inner_radius, center[1] - 0.5*bar_width, center[2], h) + p8 = gmsh.model.geo.addPoint(center[0], center[1] + inner_radius, center[2], h) + p9 = gmsh.model.geo.addPoint(center[0] - inner_radius, center[1], center[2], h) + p10 = gmsh.model.geo.addPoint(center[0], center[1] - inner_radius, center[2], h) + + # Bar end points + p11 = gmsh.model.geo.addPoint(center[0] + inner_radius - bar_length, center[1] + 0.5*bar_width, center[2], h) + p12 = gmsh.model.geo.addPoint(center[0] + inner_radius - bar_length, center[1] - 0.5*bar_width, center[2], h) + + # Create circular arcs for outer circle + c1 = gmsh.model.geo.addCircleArc(p2, p1, p3) + c2 = gmsh.model.geo.addCircleArc(p3, p1, p4) + c3 = gmsh.model.geo.addCircleArc(p4, p1, p5) + c4 = gmsh.model.geo.addCircleArc(p5, p1, p2) + + # Create circular arcs for inner circle + c5 = gmsh.model.geo.addCircleArc(p6, p1, p8) + c6 = gmsh.model.geo.addCircleArc(p8, p1, p9) + c7 = gmsh.model.geo.addCircleArc(p9, p1, p10) + c8 = gmsh.model.geo.addCircleArc(p10, p1, p7) + + # Create lines for the bar + l1 = gmsh.model.geo.addLine(p6, p11) + l2 = gmsh.model.geo.addLine(p11, p12) + l3 = gmsh.model.geo.addLine(p12, p7) + + # Create curve loops + outer_loop = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4]) + inner_loop = gmsh.model.geo.addCurveLoop([c5, c6, c7, c8, -l3, -l2, -l1]) + + # Create surface with hole + s1 = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop]) + + # Synchronize and generate mesh + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(2) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def triangle_mesh_symmetric(xc=[0., 0., 0.], r=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True): + """ + Create a triangle mesh, either symmetric or full. + + Parameters: + xc - center point coordinates [x, y, z] + r - radius of the circumscribed circle + h - mesh size + filename - output filename + vtk_out - whether to output VTK file + symmetric_mesh - if True, creates 1/3 mesh and rotates it twice. If False, creates full triangle + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/3 triangle mesh at origin first, then rotate and translate + xc_mesh = [0., 0., 0.] + # Get reference triangle points + x_tri = get_ref_triangle_points(xc_mesh, r) + v1, v2, v3 = x_tri[0], x_tri[1], x_tri[2] + + # Create points for 1/3rd mesh (triangle from center to first two vertices) + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center + p2 = gmsh.model.geo.addPoint(v1[0], v1[1], v1[2], h) # First vertex + p3 = gmsh.model.geo.addPoint(v2[0], v2[1], v2[2], h) # Second vertex + + # Create lines + l1 = gmsh.model.geo.addLine(p1, p2) # Center to first vertex + l2 = gmsh.model.geo.addLine(p2, p3) # First to second vertex + l3 = gmsh.model.geo.addLine(p3, p1) # Second vertex back to center + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + # Get mesh data for rotation + m = get_gmsh_entities() + + # Rotate twice by 120 degrees to complete the triangle + for i in range(2): + angle = (i + 1) * 2*np.pi/3 # 120 degrees in radians + gmsh_transform_general(m, 1000*(i+1), 1000000*(i+1), 1000000*(i+1), angle) + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center + gmsh_translate(xc) + + else: + # For non-symmetric mesh, create full triangle directly + x_tri = get_ref_triangle_points(xc, r) + + # Create points for the full triangle + points = [] + points.append(gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h)) # Center point + for v in x_tri: + points.append(gmsh.model.geo.addPoint(v[0], v[1], v[2], h)) + + # Create lines connecting vertices + lines = [] + for i in range(2): + lines.append(gmsh.model.geo.addLine(points[i+1], points[i+2])) + lines.append(gmsh.model.geo.addLine(points[3], points[1])) + + # Create curve loop and surface + c1 = gmsh.model.geo.addCurveLoop(lines) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + # Synchronize + gmsh.model.geo.synchronize() + + # Embed center point in surface + gmsh.model.mesh.embed(0, [points[0]], 2, s1) + + # Generate mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def cuboid_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 1., Lz = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a cuboid mesh, either symmetric or full. + xc - center point for symmetric mesh, or bottom-left-back corner for non-symmetric mesh + Lx - length in x direction + Ly - length in y direction + Lz - length in z direction + h - mesh size + symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full cuboid + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/8th of cuboid at origin, then mirror it + xc_mesh = [0., 0., 0.] + + # Create box for first octant + box = gmsh.model.occ.addBox(xc_mesh[0], xc_mesh[1], xc_mesh[2], + 0.5*Lx, 0.5*Ly, 0.5*Lz) + + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Get mesh data for mirroring + m = get_gmsh_entities() + + # Mirror the mesh to create full cuboid + # First create the other octants in the positive z hemisphere + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across yz plane + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across xz plane + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across z axis + + # Then mirror the entire upper hemisphere to create lower hemisphere + gmsh_transform(m, 4000, 4000000, 4000000, 1, 1, -1) # Mirror across xy plane + gmsh_transform(m, 5000, 5000000, 5000000, -1, 1, -1) # Mirror across yz plane + gmsh_transform(m, 6000, 6000000, 6000000, 1, -1, -1) # Mirror across xz plane + gmsh_transform(m, 7000, 7000000, 7000000, -1, -1, -1) # Mirror across z axis + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center coordinates + gmsh_translate(xc) + + else: + # Create full cuboid directly using OpenCASCADE + # Create box centered at xc + box = gmsh.model.occ.addBox(xc[0] - 0.5*Lx, xc[1] - 0.5*Ly, xc[2] - 0.5*Lz, + Lx, Ly, Lz) + + # Add center point + p1 = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2], h) + + # Synchronize + gmsh.model.occ.synchronize() + + # Embed center point in volume + gmsh.model.mesh.embed(0, [p1], 3, box) + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def ellipsoid_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, rz = 0.3, h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create an ellipsoid mesh, either symmetric or full. + xc - center point coordinates [x, y, z] + rx - radius in x direction + ry - radius in y direction + rz - radius in z direction + h - mesh size + symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full ellipsoid + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/8th of ellipsoid at origin, then mirror it + xc_mesh = [0., 0., 0.] + + # Create sphere of radius 1 + sphere = gmsh.model.occ.addSphere(xc_mesh[0], xc_mesh[1], xc_mesh[2], 1.0) + + # Scale it to create ellipsoid + gmsh.model.occ.dilate([(3, sphere)], xc_mesh[0], xc_mesh[1], xc_mesh[2], rx, ry, rz) + + # Create cutting boxes for first octant + cut_left = gmsh.model.occ.addBox(xc_mesh[0]-rx, xc_mesh[1]-ry, xc_mesh[2]-rz, rx, 2*ry, 2*rz) + cut_back = gmsh.model.occ.addBox(xc_mesh[0]-rx, xc_mesh[1]-ry, xc_mesh[2]-rz, 2*rx, ry, 2*rz) + cut_bottom = gmsh.model.occ.addBox(xc_mesh[0]-rx, xc_mesh[1]-ry, xc_mesh[2]-rz, 2*rx, 2*ry, rz) + + # Cut sphere to get first octant + gmsh.model.occ.cut([(3,sphere)], [(3,cut_left), (3,cut_back), (3,cut_bottom)]) + + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Get mesh data for mirroring + m = get_gmsh_entities() + + # Mirror the mesh to create full ellipsoid + # First create the other octants in the positive z hemisphere + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across yz plane + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across xz plane + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across z axis + + # Then mirror the entire upper hemisphere to create lower hemisphere + gmsh_transform(m, 4000, 4000000, 4000000, 1, 1, -1) # Mirror across xy plane + gmsh_transform(m, 5000, 5000000, 5000000, -1, 1, -1) # Mirror across yz plane + gmsh_transform(m, 6000, 6000000, 6000000, 1, -1, -1) # Mirror across xz plane + gmsh_transform(m, 7000, 7000000, 7000000, -1, -1, -1) # Mirror across z axis + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center coordinates + gmsh_translate(xc) + + else: + # Create full ellipsoid directly using OpenCASCADE + # First create sphere of radius 1 + sphere = gmsh.model.occ.addSphere(xc[0], xc[1], xc[2], 1.0) + + # Scale it to create ellipsoid + gmsh.model.occ.dilate([(3, sphere)], xc[0], xc[1], xc[2], rx, ry, rz) + + # Synchronize + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def cylinder_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 1., mesh_size = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create a cylinder mesh, either symmetric or full. + xc - center point coordinates [x, y, z] (center of bottom face) + r - radius of cylinder + h - height of cylinder + mesh_size - mesh element size + symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full cylinder + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/8th of cylinder at origin, then mirror it + xc_mesh = [0., 0., 0.] + + # Create full cylinder first + cyl = gmsh.model.occ.addCylinder(xc_mesh[0], xc_mesh[1], xc_mesh[2], + 0, 0, h, # height vector in z direction + r) + + # Create cutting boxes for first octant + cut_left = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2], r, 2*r, h) + cut_back = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2], 2*r, r, h) + + # Cut cylinder to get first quarter + gmsh.model.occ.cut([(3,cyl)], [(3,cut_left), (3,cut_back)]) + + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Get mesh data for mirroring + m = get_gmsh_entities() + + # Mirror the mesh to create full cylinder + # First create the other quarters + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across yz plane + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across xz plane + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across z axis + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center coordinates + gmsh_translate(xc) + + else: + # Create full cylinder directly using OpenCASCADE + cyl = gmsh.model.occ.addCylinder(xc[0], xc[1], xc[2], # base center + 0, 0, h, # height vector in z direction + r) # radius + + # Add center points at top and bottom faces + p1 = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2], mesh_size) # bottom center + p2 = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2] + h, mesh_size) # top center + + # Synchronize + gmsh.model.occ.synchronize() + + # Embed center points in respective faces + # First get all surfaces + surfaces = gmsh.model.getEntities(2) + # Find top and bottom surfaces (they are parallel to xy-plane) + for s in surfaces: + com = gmsh.model.occ.getCenterOfMass(s[0], s[1]) + if abs(com[2] - xc[2]) < 1e-6: # bottom face + gmsh.model.mesh.embed(0, [p1], 2, s[1]) + elif abs(com[2] - (xc[2] + h)) < 1e-6: # top face + gmsh.model.mesh.embed(0, [p2], 2, s[1]) + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def annulus_circle_mesh_symmetric(xc = [0., 0., 0.], r_outer = 1., r_inner = 0.5, h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True): + """ + Create an annulus (ring) mesh, either symmetric or full. + xc - center point coordinates [x, y, z] + r_outer - radius of outer circle + r_inner - radius of inner circle + h - mesh size + symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full annulus + """ + if r_inner >= r_outer: + raise ValueError("Inner radius must be smaller than outer radius") + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + xc_mesh = [0., 0., 0.] + + # add points for the quadrant of annulus circle + p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) + p2 = gmsh.model.geo.addPoint(xc_mesh[0] + r_outer, xc_mesh[1], xc_mesh[2], h) + p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + r_outer, xc_mesh[2], h) + p4 = gmsh.model.geo.addPoint(xc_mesh[0] + r_inner, xc_mesh[1], xc_mesh[2], h) + p5 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + r_inner, xc_mesh[2], h) + + l1 = gmsh.model.geo.addCircleArc(p2, p1, p3) + l2 = gmsh.model.geo.addCircleArc(p4, p1, p5) + l3 = gmsh.model.geo.addLine(p4, p2) + l4 = gmsh.model.geo.addLine(p5, p3) + + c1 = gmsh.model.geo.addCurveLoop([l3, l1, -l4, -l2]) + s1 = gmsh.model.geo.addPlaneSurface([c1]) + + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(3) + + m = get_gmsh_entities() + + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis + + # remove the duplicate nodes that will have been created on the internal + # boundaries + gmsh.model.mesh.removeDuplicateNodes() + + # translate the mesh to the specified center coordinates + gmsh_translate(xc) + + else: + + # Outer circle points + p1 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h) + p2 = gmsh.model.geo.addPoint(xc[0] + r_outer, xc[1], xc[2], h) + p3 = gmsh.model.geo.addPoint(xc[0], xc[1] + r_outer, xc[2], h) + p4 = gmsh.model.geo.addPoint(xc[0] - r_outer, xc[1], xc[2], h) + p5 = gmsh.model.geo.addPoint(xc[0], xc[1] - r_outer, xc[2], h) + + # Inner circle points + p6 = gmsh.model.geo.addPoint(xc[0] + r_inner, xc[1], xc[2], h) + p7 = gmsh.model.geo.addPoint(xc[0], xc[1] + r_inner, xc[2], h) + p8 = gmsh.model.geo.addPoint(xc[0] - r_inner, xc[1], xc[2], h) + p9 = gmsh.model.geo.addPoint(xc[0], xc[1] - r_inner, xc[2], h) + + # Create circular arcs for outer circle + c1 = gmsh.model.geo.addCircleArc(p2, p1, p3) + c2 = gmsh.model.geo.addCircleArc(p3, p1, p4) + c3 = gmsh.model.geo.addCircleArc(p4, p1, p5) + c4 = gmsh.model.geo.addCircleArc(p5, p1, p2) + + # Create circular arcs for inner circle + c5 = gmsh.model.geo.addCircleArc(p6, p1, p7) + c6 = gmsh.model.geo.addCircleArc(p7, p1, p8) + c7 = gmsh.model.geo.addCircleArc(p8, p1, p9) + c8 = gmsh.model.geo.addCircleArc(p9, p1, p6) + + # Create curve loops + outer_loop = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4]) + inner_loop = gmsh.model.geo.addCurveLoop([c5, c6, c7, c8]) + + #gmsh.fltk.run() + + # surface + s1 = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop]) + + # designate surface as physical surface + gmsh.model.addPhysicalGroup(2, [s1], 1) + + # Synchronize and generate mesh + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(2) + + # Write output files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def annulus_rectangle_mesh(xc=[0., 0., 0.], Lx=1., Ly=1., hole_Lx=0.3, hole_Ly=0.3, h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True): + """ + Create a rectangular mesh with a rectangular hole in the center (annulus rectangle). + + Parameters + ---------- + xc : list + Center coordinates [x, y, z] + Lx : float + Length of outer rectangle in x direction + Ly : float + Length of outer rectangle in y direction + hole_Lx : float + Length of inner hole in x direction + hole_Ly : float + Length of inner hole in y direction + h : float + Mesh size + filename : str + Output filename without extension + vtk_out : bool + If True, also writes a VTK file + symmetric_mesh : bool + If True, creates 1/4 mesh and mirrors it. If False, creates full mesh + """ + # Initialize gmsh + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + if symmetric_mesh: + # Create 1/4 mesh at origin first, then mirror and translate + xc_mesh = [0., 0., 0.] + + # Create points for outer rectangle (1/4) + p1 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1], xc_mesh[2], h) + p2 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1] + 0.5*Ly, xc_mesh[2], h) + p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + 0.5*Ly, xc_mesh[2], h) + p4 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1], xc_mesh[2], h) + p5 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1] + 0.5*hole_Ly, xc_mesh[2], h) + p6 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + 0.5*hole_Ly, xc_mesh[2], h) + + # Create lines + l1 = gmsh.model.geo.addLine(p1, p2) + l2 = gmsh.model.geo.addLine(p2, p3) + l3 = gmsh.model.geo.addLine(p3, p6) + l4 = gmsh.model.geo.addLine(p6, p5) + l5 = gmsh.model.geo.addLine(p5, p4) + l6 = gmsh.model.geo.addLine(p4, p1) + + # Create curve loops + cl = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5, l6]) + + # Create plane surface with hole + s = gmsh.model.geo.addPlaneSurface([cl]) + + # Synchronize and generate mesh + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(2) + + # Get mesh data for mirroring + m = get_gmsh_entities() + + # Mirror the mesh in all quadrants + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across y-axis + gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across x-axis + gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across origin + + # Remove duplicate nodes + gmsh.model.mesh.removeDuplicateNodes() + + # Translate to specified center coordinates + gmsh_translate(xc) + + else: + # Create points for outer rectangle + p1 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom left + p2 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom right + p3 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top right + p4 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top left + + # Create points for inner hole + p5 = gmsh.model.geo.addPoint(xc[0] - 0.5*hole_Lx, xc[1] - 0.5*hole_Ly, xc[2], h) # Bottom left + p6 = gmsh.model.geo.addPoint(xc[0] + 0.5*hole_Lx, xc[1] - 0.5*hole_Ly, xc[2], h) # Bottom right + p7 = gmsh.model.geo.addPoint(xc[0] + 0.5*hole_Lx, xc[1] + 0.5*hole_Ly, xc[2], h) # Top right + p8 = gmsh.model.geo.addPoint(xc[0] - 0.5*hole_Lx, xc[1] + 0.5*hole_Ly, xc[2], h) # Top left + + # Create lines for outer rectangle + l1 = gmsh.model.geo.addLine(p1, p2) # Bottom + l2 = gmsh.model.geo.addLine(p2, p3) # Right + l3 = gmsh.model.geo.addLine(p3, p4) # Top + l4 = gmsh.model.geo.addLine(p4, p1) # Left + + # Create lines for inner hole + l5 = gmsh.model.geo.addLine(p5, p6) # Bottom + l6 = gmsh.model.geo.addLine(p6, p7) # Right + l7 = gmsh.model.geo.addLine(p7, p8) # Top + l8 = gmsh.model.geo.addLine(p8, p5) # Left + + # Create curve loops + outer_loop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4]) + inner_loop = gmsh.model.geo.addCurveLoop([l5, l6, l7, l8]) + + # Create plane surface with hole + s1 = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop]) + + # Synchronize and generate mesh + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(2) + + # Write mesh files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + # Finalize gmsh + gmsh.finalize() + +def open_rectangle_mesh(xc = [0., 0., 0.], Lx = 1., Ly = 1., hole_Lx = 0.5, hole_Ly = 0.5, h = 0.1, filename = 'mesh', vtk_out = False): + """ + Create a rectangular mesh with a rectangular hole in the center (annulus rectangle). + + Parameters + ---------- + xc : list + Center coordinates [x, y, z] + Lx : float + Length of outer rectangle in x direction + Ly : float + Length of outer rectangle in y direction + hole_Lx : float + Length of inner hole in x direction + hole_Ly : float + Length of inner hole in y direction + h : float + Mesh size + filename : str + Output filename without extension + vtk_out : bool + If True, also writes a VTK file + """ + # Initialize gmsh + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + # Annulus rectangle with top removed + # xc is the center of the cocentric rectangles + # we mesh the right half of the annulus rectangle + xc_mesh = [0., 0., 0.] + points = [] + points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1] - 0.5*Ly, xc_mesh[2], h)) + points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1] + 0.5*Ly, xc_mesh[2], h)) + points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1] + 0.5*Ly, xc_mesh[2], h)) + points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1] - 0.5*hole_Ly, xc_mesh[2], h)) + points.append(gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] - 0.5*hole_Ly, xc_mesh[2], h)) + points.append(gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] - 0.5*Ly, xc_mesh[2], h)) + + lines = [] + for i in range(len(points) - 1): + lines.append(gmsh.model.geo.addLine(points[i], points[i + 1])) + lines.append(gmsh.model.geo.addLine(points[-1], points[0])) + cl = gmsh.model.geo.addCurveLoop(lines) + s = gmsh.model.geo.addPlaneSurface([cl]) + gmsh.model.geo.synchronize() + gmsh.model.mesh.generate(2) + + m = get_gmsh_entities() + + # mirror the mesh across the x-axis + gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) + gmsh.model.mesh.removeDuplicateNodes() + gmsh_translate(xc) + + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +def open_pipe_mesh(xc=[0., 0., 0.], axis=[0., 0., 1.], length=2., outer_radius=1., wall_thickness=0.1, h=0.1, filename='mesh', vtk_out=False): + """ + Create a 3D pipe mesh with specified axis, closed bottom and open top. + The pipe is created by: + 1. Creating an outer cylinder + 2. Subtracting an inner cylinder to create walls + 3. Subtracting the top surface to create the opening + + Parameters + ---------- + xc : list + Center coordinates [x, y, z] of the base center + axis : list + Axis vector defining pipe orientation [ax, ay, az] + length : float + Length of the pipe along axis + outer_radius : float + Outer radius of the pipe + wall_thickness : float + Thickness of the pipe wall and bottom + h : float + Mesh size + filename : str + Output filename without extension + vtk_out : bool + If True, also writes a VTK file + """ + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + # Normalize axis vector + axis_norm = np.sqrt(axis[0]**2 + axis[1]**2 + axis[2]**2) + if axis_norm == 0: + raise ValueError("Axis vector cannot be zero") + axis = [x/axis_norm for x in axis] + + # Create outer cylinder + l = length + wall_thickness + outer_cylinder = gmsh.model.occ.addCylinder(xc[0], xc[1], xc[2], + axis[0]*l, axis[1]*l, axis[2]*l, + outer_radius) + + # Create inner cylinder (to be subtracted) + inner_cylinder = gmsh.model.occ.addCylinder(xc[0]+wall_thickness*axis[0], xc[1]+wall_thickness*axis[1], xc[2]+wall_thickness*axis[2], + axis[0]*l, axis[1]*l, axis[2]*l, + outer_radius - wall_thickness) + + # Create a box slightly larger than the cylinder at the top to cut off the top surface + # First create a box aligned with coordinate axes + dx = 2 * outer_radius + box = gmsh.model.occ.addBox(xc[0] - dx, xc[1] - dx, xc[2] + l - wall_thickness, + 2*dx, 2*dx, 2*wall_thickness) + + # If axis is not aligned with z, rotate the box to align with axis + if not (abs(axis[0]) < 1e-10 and abs(axis[1]) < 1e-10): + # Calculate rotation angle and axis + z_axis = [0, 0, 1] + rotation_axis = np.cross(z_axis, axis) + rotation_angle = np.arccos(np.dot(z_axis, axis)) + + # Rotate the box around center point + gmsh.model.occ.rotate([(3, box)], + xc[0], xc[1], xc[2], + rotation_axis[0], rotation_axis[1], rotation_axis[2], + rotation_angle) + + # First subtract inner cylinder from outer cylinder + pipe = gmsh.model.occ.cut([(3, outer_cylinder)], [(3, inner_cylinder), (3, box)]) + + # Synchronize before meshing + gmsh.model.occ.synchronize() + + # Set mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h) + + # Generate 3D mesh + gmsh.model.mesh.generate(3) + + # Write mesh files + gmsh.write(filename + '.msh') + if vtk_out: + gmsh.write(filename + '.vtk') + + gmsh.finalize() + +##-------------------------------------------------------## +##-------------------------------------------------------## +if __name__ == "__main__": + inp_dir = './' + + test_meshes = ['circle', 'ellipse', 'sphere', 'cuboid', 'ellipsoid', 'rectangle', 'hexagon', 'drum2d', 'triangle', 'polygon', 'cylindrical2d_wall', 'cylinder', 'annulus_circle', 'annulus_rectangle', 'open_rectangle', 'open_pipe'] + test_meshes = ['open_pipe'] + symm_flag = 0 # 0 - both, 1 - symmetric, 2 - non-symmetric + + for mesh in test_meshes: + if mesh == 'circle': + if symm_flag == 0 or symm_flag == 1: + circle_mesh_symmetric(xc = [-3., -3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + circle_mesh_symmetric(xc = [-3., -3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'ellipse': + if symm_flag == 0 or symm_flag == 1: + ellipse_mesh_symmetric(xc = [-1., -1., 0.], rx = 1., ry = 0.5, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + ellipse_mesh_symmetric(xc = [-1., -1., 0.], rx = 1., ry = 0.5, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'sphere': + if symm_flag == 0 or symm_flag == 1: + sphere_mesh_symmetric(xc = [-5., -5., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + sphere_mesh_symmetric(xc = [-5., -5., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'rectangle': + if symm_flag == 0 or symm_flag == 1: + rectangle_mesh_symmetric(xc = [1., 1., 0.], Lx = 1., Ly = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + rectangle_mesh_symmetric(xc = [1., 1., 0.], Lx = 1., Ly = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'hexagon': + if symm_flag == 0 or symm_flag == 1: + hexagon_mesh_symmetric(xc = [3., 3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + hexagon_mesh_symmetric(xc = [3., 3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'drum2d': + if symm_flag == 0 or symm_flag == 1: + drum2d_mesh_symmetric(xc = [5., 5., 0.], r = 1., width = 0.5, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + drum2d_mesh_symmetric(xc = [5., 5., 0.], r = 1., width = 0.5, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'triangle': + if symm_flag == 0 or symm_flag == 1: + triangle_mesh_symmetric(xc = [7., 7., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + triangle_mesh_symmetric(xc = [7., 7., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'polygon': + theta = np.pi/6. + R, a = 1., 0.25 + v1 = [0., 0., 0.] + v2 = [R*np.cos(0.5*theta), -R*np.sin(0.5*theta), 0.] + v4 = [R*np.cos(0.5*theta), R*np.sin(0.5*theta), 0.] + v3 = [R + a, 0., 0.] + if symm_flag == 0 or symm_flag == 1: + polygon_mesh_symmetric(points = [v1, v2, v3, v4], theta = theta, xc = [9., 9., 9.], h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + polygon_mesh_symmetric(points = [v1, v2, v3, v4], theta = theta, xc = [9., 9., 9.], h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'cylindrical2d_wall': + cylindrical2d_wall_mesh( + center=[0., 0., 0.], + outer_radius=1.0, + inner_radius=0.8, + bar_width=0.2, + bar_length=0.3, + h=0.05, + filename='./' + mesh + '_sym', + vtk_out=True + ) + elif mesh == 'cuboid': + if symm_flag == 0 or symm_flag == 1: + cuboid_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 0.5, Lz = 0.3, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + cuboid_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 0.5, Lz = 0.3, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'ellipsoid': + if symm_flag == 0 or symm_flag == 1: + ellipsoid_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, rz = 0.3, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + ellipsoid_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, rz = 0.3, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'cylinder': + if symm_flag == 0 or symm_flag == 1: + cylinder_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 2., mesh_size = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + cylinder_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 2., mesh_size = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'annulus_circle': + if symm_flag == 0 or symm_flag == 1: + annulus_circle_mesh_symmetric(xc = [0., 0., 0.], r_outer = 1., r_inner = 0.5, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + annulus_circle_mesh_symmetric(xc = [0., 0., 0.], r_outer = 1., r_inner = 0.5, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'annulus_rectangle': + if symm_flag == 0 or symm_flag == 1: + annulus_rectangle_mesh(xc = [0., 0., 0.], Lx = 1., Ly = 0.8, hole_Lx = 0.2, hole_Ly = 0.4, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True) + if symm_flag == 0 or symm_flag == 2: + annulus_rectangle_mesh(xc = [0., 0., 0.], Lx = 1., Ly = 0.8, hole_Lx = 0.2, hole_Ly = 0.4, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False) + elif mesh == 'open_rectangle': + open_rectangle_mesh(xc = [0., 0., 0.], Lx = 1.2, Ly = 0.8, hole_Lx = 1., hole_Ly = 0.6, h = 0.1, filename = './' + mesh, vtk_out = True) + elif mesh == 'open_pipe': + open_pipe_mesh(xc=[0., 0., 0.], axis=[1., 1., 0.], length=2., outer_radius=0.5, + wall_thickness=0.2, h=0.1, filename='./' + mesh, vtk_out=True) + else: + print(f"Mesh {mesh} not found") + \ No newline at end of file diff --git a/tools/python_utils/util.py b/tools/python_utils/util.py index a099ada3a..55c4ee9d9 100644 --- a/tools/python_utils/util.py +++ b/tools/python_utils/util.py @@ -27,7 +27,7 @@ def print_dbl_list(arg, prefix = ""): return a def print_int_list(arg, prefix = ""): - a = prefix + '[' + print_list(arg, '%4.6e') + ']\n' + a = prefix + '[' + print_list(arg, '%d') + ']\n' return a def print_point_gmsh(arg, n): @@ -118,394 +118,35 @@ def get_max(l): l = np.array(l) return l.max() -def rotate(p, theta, axis): - """ - Returns rotation of vector about specified axis by specified angle - - Parameters - ---------- - p: list - Coordinates of vector - theta : float - Angle of rotation - axis : list - Axis of rotation - - Returns - ------- - list - Coordinates of rotated vector - """ - p_np, axis_np = np.array(p), np.array(axis) - ct, st = np.cos(theta), np.sin(theta) - p_dot_n = np.dot(p_np,axis_np) - n_cross_p = np.cross(axis_np, p_np) - - return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p - -def get_ref_hex_points(center, radius, add_center = False): - """ - Returns size points on reference hexagon - - Reference hexagon: - - v3 v2 - + + - - - + o + - v4 x v1 - - + + - v5 v6 - - Axis is a vector from x to v1 and radius is distance between x and v1. - - Parameters - ---------- - center: list - Coordinates of center of hexagon - radius : float - Radius of hexagon - add_center : bool - True if we include the center to the returned list (first point in the returned list will be the center if the flag is true) - - Returns - ------- - list - Coordinates of points - """ - axis = [1., 0., 0.] - rotate_axis = [0., 0., 1.] - points = [] - if add_center: - points.append(center) - - for i in range(6): - xi = rotate(axis, i*np.pi/3., rotate_axis) - points.append([center[i] + radius * xi[i] for i in range(3)]) - - return points - -def get_ref_drum_points(center, radius, width, add_center = False): - """ - Returns size points on reference concave polygon (we refer to it as drum) - Reference concave polygon: - - v3 v2 - + -------------------------------- + - \ / - \ / - + o + - /v4 x v1 \ - / \ - + -------------------------------- + - v5 v6 +def write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, zone_string, Kn): + inpf.write(" Zone_%s:\n" % (zone_string)) + inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor)) - Axis is a vector from x to v1, radius is distance between x and v2, and width of neck is the distance between v2 and v4. - - Parameters - ---------- - center: list - Coordinates of center of polygon - radius : float - Radius of polygon (distance between center x and vertex v2) - width : float - Width of neck, i.e. distance between vertex v2 and v4 - add_center : bool - True if we include the center to the returned list (first point in the returned list will be the center if the flag is true) - - Returns - ------- - list - Coordinates of points - """ - axis = [1., 0., 0.] - rotate_axis = [0., 0., 1.] - - x1 = rotate(axis, np.pi/3., rotate_axis) - x2 = rotate(axis, -np.pi/3., rotate_axis) - - points = [] - if add_center: - points.append(center) - - # v1 - points.append([center[i] + width*0.5*axis[i] for i in range(3)]) - # v2 - points.append([center[i] + radius*x1[i] for i in range(3)]) - # v3 - points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)]) - # v4 - points.append([center[i] - width*0.5*axis[i] for i in range(3)]) - # v5 - v6 = [center[i] + radius*x2[i] for i in range(3)] - points.append([v6[i] - radius*axis[i] for i in range(3)]) - # v6 - points.append(v6) - - return points - - -def does_intersect(p, particles, padding): - """ - Returns true if particle p intersects or is near enough to existing particles - - Parameters - ---------- - p : list - Coordinates of center and radius of particle [x,y,z,r] - particles : list - List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius. - padding: float - Minimum distance between circle boundaries such that if two circles - - Returns - ------- - bool - True if particle intersects or is near enough to one of the particle in the list - """ - if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p)) - if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles)) - if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding)) - - for q in particles: - if len(q) < 4: raise Exception('q = {} in particles must have atleast 4 elements'.format(q)) - pq = np.array([p[i] - q[i] for i in range(3)]) - if np.linalg.norm(pq) <= p[3] + q[3] + padding: - return True - return False - -def does_intersect_rect(p, particles, padding, rect, is_3d = False): - """ - Returns true if particle p is sufficiently close or outside the rectangle (in 2d) or cuboid (in 3d) - - Parameters - ---------- - p : list - Coordinates of center and radius of particle [x,y,z,r] - particles : list - List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius. - padding: float - Minimum distance between circle boundaries such that if two circles - rect: list - Coordinates of left-bottom and right-top corner points of rectangle (2d) or cuboid (3d). E.g. [x1 y1, z1, x2, y2, z2] - is_3d: bool - True if we are dealing with cuboid - - Returns - ------- - bool - True if particle intersects or is near enough to the rectangle - """ - if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p)) - if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles)) - if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding)) - if len(rect) < 6: raise Exception('rect = {} must have 6 elements'.format(rect)) - - pr = [p[0] - p[3], p[1] - p[3], p[2], p[0] + p[3], p[1] + p[3], p[2]] - if is_3d: - pr[2] -= p[3] - pr[5] += p[3] - - if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding: - if is_3d: - if pr[2] < rect[2] + padding or pr[5] > rect[5] - padding: - return True - else: - return True - return False - -def generate_circle_gmsh_input(filename, center, radius, mesh_size, pp_tag = None): - """ - Creates .geo file for discretization of circle using gmsh - - Parameters - ---------- - filename : str - Filename of .geo file to be created - center : list - Coordinates of center of circle - radius: float - Radius of circle - mesh_size: float - Mesh size - pp_tag: str, optional - Postfix .geo file with this tag - """ - pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else '' - geof = open(filename + pp_tag_str + '.geo','w') - gmsh_file_hdr(geof) - - ## points - points = [] - points.append([center[0], center[1], center[2]]) - points.append([center[0] + radius, center[1], center[2]]) - points.append([center[0] - radius, center[1], center[2]]) - points.append([center[0], center[1] + radius, center[2]]) - points.append([center[0], center[1] - radius, center[2]]) - for i in range(len(points)): - geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1)) - - ## circular arc - geof.write(print_cir_gmsh([2,1,4], 1)) - geof.write(print_cir_gmsh([4,1,3], 2)) - geof.write(print_cir_gmsh([3,1,5], 3)) - geof.write(print_cir_gmsh([5,1,2], 4)) - - ## line loop - geof.write(print_lineloop_gmsh([2, 3, 4, 1], 1)) - - ## plane surface - geof.write("Plane Surface(1) = {1};\n") - - ## add center point to plane surface - geof.write("Point{1} In Surface {1};") - - ## close file - geof.close() - -def generate_rectangle_gmsh_input(filename, rectangle, mesh_size, pp_tag = None): - """ - Creates .geo file for discretization of rectangle using gmsh - - Parameters - ---------- - filename : str - Filename of .geo file to be created - rectangle : list - Coordinates of left-bottom and right-top corner points of rectangle - mesh_size: float - Mesh size - pp_tag: str, optional - Postfix .geo file with this tag - """ - pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else '' - geof = open(filename + pp_tag_str + '.geo','w') - gmsh_file_hdr(geof) - - ## points - points = [] - points.append([rectangle[0], rectangle[1], rectangle[2]]) - points.append([rectangle[3], rectangle[1], rectangle[2]]) - points.append([rectangle[3], rectangle[4], rectangle[2]]) - points.append([rectangle[0], rectangle[4], rectangle[2]]) - for i in range(len(points)): - geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1)) - - ## lines - for i in range(4): - if i < 3: - geof.write(print_line_gmsh([i+1,i+2], i+1)) - else: - geof.write(print_line_gmsh([i+1,1], i+1)) - - ## line loop - geof.write(print_lineloop_gmsh([i+1 for i in range(4)], 1)) - - ## plane surface - geof.write("Plane Surface(1) = {1};\n") - - ## add center point to plane surface - geof.write("Point{1} In Surface {1};") - - ## physical surface - tag = '"' + "a" + '"' - geof.write("Physical Surface(%s) = {1};\n" % (tag)) - - ## close file - geof.close() - -def generate_hexagon_gmsh_input(filename, center, radius, mesh_size, pp_tag = None): - """ - Creates .geo file for discretization of hexagon using gmsh - - Parameters - ---------- - filename : str - Filename of .geo file to be created - center : list - Coordinates of center of hexagon - radius: float - Radius - mesh_size: float - Mesh size - pp_tag: str, optional - Postfix .geo file with this tag - """ - pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else '' - geof = open(filename + pp_tag_str + '.geo','w') - gmsh_file_hdr(geof) - - ## points - points = get_ref_hex_points(center, radius, True) - for i in range(len(points)): - geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1)) - - ## lines - for i in range(6): - if i < 5: - geof.write(print_line_gmsh([i+2,i+3], i+1)) - else: - geof.write(print_line_gmsh([i+2,2], i+1)) - - ## line loop - geof.write(print_lineloop_gmsh([i+1 for i in range(6)], 1)) - - ## plane surface - geof.write("Plane Surface(1) = {1};\n") - - ## add center point to plane surface - geof.write("Point{1} In Surface {1};") - - ## close file - geof.close() - - -def generate_drum_gmsh_input(filename, center, radius, width, mesh_size, pp_tag = None): - """ - Creates .geo file for discretization of drum (concave polygon, see get_ref_drum_points) using gmsh - - Parameters - ---------- - filename : str - Filename of .geo file to be created - center : list - Coordinates of center - radius: float - Radius - width: float - Neck width (see get_ref_drum_points()) - mesh_size: float - Mesh size - pp_tag: str, optional - Postfix .geo file with this tag - """ - pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else '' - geof = open(filename + pp_tag_str + '.geo','w') - gmsh_file_hdr(geof) - - ## points - points = get_ref_drum_points(center, radius, width, True) - for i in range(len(points)): - geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1)) - - ## lines - for i in range(6): - if i < 5: - geof.write(print_line_gmsh([i+2,i+3], i+1)) - else: - geof.write(print_line_gmsh([i+2,2], i+1)) - - ## line loop - geof.write(print_lineloop_gmsh([i+1 for i in range(6)], 1)) - - ## plane surface - geof.write("Plane Surface(1) = {1};\n") - - ## add center point to plane surface - geof.write("Point{1} In Surface {1};") - - ## close file - geof.close() + if damping_active == False: + inpf.write(" Damping_On: false\n") + if friction_active == False: + inpf.write(" Friction_On: false\n") + + inpf.write(" Kn: %4.6e\n" % (Kn)) + inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps)) + inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff)) + inpf.write(" Kn_Factor: %4.6e\n" % (Kn_factor)) + inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor)) + +def write_material_zone_part(inpf, zone_string, horizon, rho, K, G, Gc): + inpf.write(" Zone_%s:\n" % (zone_string)) + inpf.write(" Type: PDState\n") + inpf.write(" Horizon: %4.6e\n" % (horizon)) + inpf.write(" Density: %4.6e\n" % (rho)) + inpf.write(" Compute_From_Classical: true\n") + inpf.write(" K: %4.6e\n" % (K)) + inpf.write(" G: %4.6e\n" % (G)) + inpf.write(" Gc: %4.6e\n" % (Gc)) + inpf.write(" Influence_Function:\n") + inpf.write(" Type: 1\n") + +def copy_contact_zone(inpf, zone_numbers, zone_copy_from): + for s in zone_numbers: + inpf.write(" Zone_%d:\n" % (s)) + inpf.write(" Copy_Contact_Data: " + print_int_list(zone_copy_from)) \ No newline at end of file