Looking more closely at this post, one thing strikes me:
Does the 1D structure align with the edges of the 3D space?
If not, I don’t see how one can expect an accurate interpolation from 1D to 3D, as all fields would have to be extrapolated to match the interpolation points of the 3D structure.
In your code, it seems like they are matching, and they give me a result I would expect with
Code:
import numpy as np
import gmsh
import meshio
# Initialize Gmsh
gmsh.initialize()
# Create a new model
gmsh.model.add("model3D")
gmsh.option.setNumber("Geometry.OCCBooleanSimplify", 0)
mesh_name = "pMM2x"
# geom.characteristic_length_max = mesh_size
# Define parameters
r = 0.05 # Small radius
R = 0.25 # Big radius
H = 1.0 # Height
a = 4 * H / np.sqrt(3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 2 * r)
def add_link(ϕ, rotate_axis):
# Define the three points
points = [
[0.0, 0.0, 0.0],
[r, 0.0, 0.0], # Point 1
[R, 0.0, H / 2.0], # Point 2
[r, 0.0, H], # Point 3
[0.0, 0.0, H], # Point 4
]
# Define the axis of rotation (e.g., the z-axis)
axis = [0.0, 0.0, 1.0]
point_on_axis = [0.0, 0.0, 0.0]
# Add points to the model
geom_points = [gmsh.model.occ.addPoint(*p) for p in points]
sym_axis = gmsh.model.occ.addLine(geom_points[-1], geom_points[0])
# print(sym_axis)
# Create lines connecting the points
lines = []
for i in range(len(geom_points) - 1):
lines.append(gmsh.model.occ.addLine(geom_points[i], geom_points[i + 1]))
lines.append(gmsh.model.occ.addLine(geom_points[-1], geom_points[0]))
# Create a wire (curve loop) from the lines
wire = gmsh.model.occ.addCurveLoop(lines)
# Create a surface from the wire
surface = gmsh.model.occ.addPlaneSurface([wire])
# Rotate the surface around the specified axis to create a 3D body
θ = 2 * np.pi / 3
s = surface
bodies = []
for i in range(3):
# Revolve the surface
result = gmsh.model.occ.revolve(
[(2, s)], # Surface to revolve
point_on_axis[0],
point_on_axis[1],
point_on_axis[2], # Point on axis
axis[0],
axis[1],
axis[2], # Axis of rotation
θ, # Angle of rotation
)
# Extract the resulting entities
top = result[0][1] # The new surface created by the revolution
body = result[1][1] # The 3D volume created by the revolution
bodies.append(body)
s = top # Update the surface for the next revolution
# Add spheres at the ends
s1 = gmsh.model.occ.addSphere(points[0][0], points[0][1], points[0][2], r)
s2 = gmsh.model.occ.addSphere(points[-1][0], points[-1][1], points[-1][2], r)
# Combine all bodies into one
b = gmsh.model.occ.fuse([(3, s1), (3, s2)], [(3, body) for body in bodies])
# Rotate the combined body by the specified angle and axis
gmsh.model.occ.rotate(
b[0], # The combined body
point_on_axis[0],
point_on_axis[1],
point_on_axis[2], # Point on axis
rotate_axis[0],
rotate_axis[1],
rotate_axis[2], # Axis of rotation
ϕ, # Angle of rotation
)
gmsh.model.occ.rotate(
[(1, sym_axis)], # The combined body
point_on_axis[0],
point_on_axis[1],
point_on_axis[2], # Point on axis
rotate_axis[0],
rotate_axis[1],
rotate_axis[2], # Axis of rotation
ϕ, # Angle of rotation
)
return b[0], sym_axis # Return the tag of the resulting body
def create_metaatom(translate):
normal_planes = [
[1.0, -1.0, 0.0],
[1.0, -1.0, 0.0],
[1.0, 1.0, 0.0],
[1.0, 1.0, 0.0],
]
ϕ = [
np.arccos(1 / np.sqrt(3)),
2.0 * np.pi - np.arccos(1 / np.sqrt(3)),
np.arccos(-1 / np.sqrt(3)),
2.0 * np.pi - np.arccos(-1 / np.sqrt(3)),
]
links = []
sk = []
for i, n in enumerate(normal_planes):
link, ax = add_link(ϕ[i], n)
# print(ax)
gmsh.model.occ.translate([(1, ax)], translate[0], translate[1], translate[2])
links.append(link[0][:])
sk.append(ax)
# Combine all links into one meta-atom
meta_atom = gmsh.model.occ.fuse([links[0]], links[1:])
# Translate the meta-atom
gmsh.model.occ.translate(meta_atom[0], translate[0], translate[1], translate[2])
return meta_atom[0][:], sk # Return the tag of the resulting meta-atom
def create_SC_site(center):
translations = [
[0.0, 0.0, 0.0],
[a / 2.0, a / 2.0, 0.0],
[0.0, a / 2.0, a / 2.0],
[a / 2.0, 0.0, a / 2.0],
]
unit_cell = []
skeleton = []
for t in translations:
l, sk = create_metaatom(t)
for link in sk:
gmsh.model.occ.translate([(1, link)], center[0], center[1], center[2])
# gmsh.model.occ.translate(sk, center[0], center[1], center[2])
unit_cell.append(l[0][:])
skeleton += sk
# skeleton.append(sk[0][:])
# Combine all meta-atoms into one unit cell
u_cell = gmsh.model.occ.fuse([unit_cell[0]], unit_cell[1:])
# cell_skeleton = gmsh.model.occ.fuse([skeleton[0]], skeleton[1:])
# Translate the unit cell
gmsh.model.occ.translate(u_cell[0], center[0], center[1], center[2])
# gmsh.model.occ.translate(cell_skeleton[0], center[0], center[1], center[2])
return u_cell[0][:], skeleton # Return the tag of the resulting unit cell
# # Create two unit cells
N = 2
MM = []
MM_skeleton = []
for i in range(N):
# for j in range(N):
# for k in range(N):
l, bone = create_SC_site([i * a, 0.0, 0.0])
MM.append(l[0][:])
MM_skeleton += bone
# #l1 = create_SC_site([a, 0.0, 0.0])
# print(MM)
union = gmsh.model.occ.fuse([MM[0]], MM[1:])
# skeleton = gmsh.model.occ.fuse([MM_skeleton[0]], MM_skeleton[1:])
print(MM_skeleton)
# Synchronize the model
gmsh.model.occ.synchronize()
marker = 1
gmsh.model.addPhysicalGroup(3, [union[0][:][0][1]], marker)
gmsh.model.addPhysicalGroup(1, MM_skeleton, marker)
# Save the geometry to a file
gmsh.write(mesh_name + ".geo_unrolled")
# Generate the mesh (optional)
gmsh.model.mesh.generate()
# Save the mesh to a file (optional)
gmsh.write(mesh_name + ".msh")
# Finalize Gmsh
gmsh.finalize()
msh = meshio.read(mesh_name + ".msh")
# Transform mesh to xdmf file. No boundary conditions are necessary so no surface meshes are imported
body_element_tag = "tetra"
line_data = msh.cell_data_dict["gmsh:physical"][body_element_tag]
meshio.write(
mesh_name + ".xdmf",
meshio.Mesh(
points=msh.points,
cells={body_element_tag: msh.cells_dict[body_element_tag]},
cell_data={"bnd_marker": [line_data]},
),
)
surface_element_tag = "line"
line_data = msh.cell_data_dict["gmsh:physical"][surface_element_tag]
meshio.write(
mesh_name + "_skeleton.xdmf",
meshio.Mesh(
points=msh.points,
cells={surface_element_tag: msh.cells_dict[surface_element_tag]},
cell_data={"bnd_marker": [line_data]},
),
)
import dolfinx
import numpy as np
from mpi4py import MPI
from basix.ufl import element
mesh_name = "pMM2x"
# Import Meshes from xdmf file
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, mesh_name + ".xdmf", "r") as xdmf_infile:
body = xdmf_infile.read_mesh(name="Grid")
with dolfinx.io.XDMFFile(
MPI.COMM_WORLD, mesh_name + "_skeleton.xdmf", "r"
) as xdmf_infile:
skeleton = xdmf_infile.read_mesh(name="Grid")
# Define relevant vector elements and spaces
el0 = element("Lagrange", skeleton.basix_cell(), 1)
V0 = dolfinx.fem.functionspace(skeleton, el0)
el1 = element("Lagrange", body.basix_cell(), 1)
V1 = dolfinx.fem.functionspace(body, el1)
padding = 1e-14
t = dolfinx.fem.Function(V0, name="Tangent")
t.interpolate(lambda x: x[0] + 2 * x[1] - x[2])
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "skeleton_tangent.xdmf", "w") as xdmf_outfile:
xdmf_outfile.write_mesh(skeleton)
t.name = "Tangent"
xdmf_outfile.write_function(t)
fine_mesh_cell_map = body.topology.index_map(body.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
interpolation_data = dolfinx.fem.create_interpolation_data(
V1, V0, cells, padding=padding
)
n = dolfinx.fem.Function(V1, name="Director")
n.interpolate_nonmatching(t, cells, interpolation_data)
n.x.scatter_forward()
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "body_tangent.xdmf", "w") as xdmf_outfile:
xdmf_outfile.write_mesh(body)
n.name = "Tangent"
xdmf_outfile.write_function(n)
If you turn up the padding, you get something different (padding = 1)
as one then extrapolates data from the closest part of the skeleton for each cell of the mesh.