Dear Dolfinx Developers,
First and foremost, I would like to express my deepest gratitude and respect for the outstanding work you have accomplished with Dolfinx.
As a researcher in architectural engineering, I am currently working on a research paper utilizing Dolfinx. In our field, we commonly model 1D steel reinforcements as line elements embedded within 3D concrete structures.
To achieve this, I have modeled a 3D concrete volume and 1D steel reinforcements using GMSH, as shown in the code below. However, when I attempt to execute the code, it does not run as expected.
Would it be possible to seek your advice on resolving this issue?
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 8 18:54:54 2025
@author: hoseong
"""
import gmsh
import math
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_mesh
import numpy as np
import ufl
from mpi4py import MPI
import meshio
from dolfinx.io import gmshio
import pyvista as pv
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
########################Mesh creation#######################
gmsh.initialize()
gmsh.model.add('DFG 3D')
#define dimension
length_member = 500
width_member = 100
depth_member = 100
cover_spacing_longitudinal_rebar = 20
embedded_length_longitudinal_rebar = 400
diameter_longitudinal_rebar = 8
#####tag list#####
#volume 1X0
tag_volume_concrete = 130
#Line 1X000
tag_line_longitudinal_rebar = 11000
tag_line_transverse_rebar = 12000
#Point 1X0000
tag_point_count_longitudinal_rebar = 120000
tag_point_count_transverse_rebar = 120000
#volume 2X0
marker_volume_longitudinal_rebar = 210
marker_volume_transverse_rebar = 220
marker_volume_concrete = 230
#Surface 2X00
marker_surface_longitudinal_rebar = 2100
marker_surface_transverse_rebar = 2200
##### model concrete #####
domain_concrete = gmsh.model.occ.addBox(0, 0, 0, length_member, width_member, depth_member, tag=tag_volume_concrete)
##### generate concrete #####
gmsh.model.occ.synchronize()
##### model rebar #####
gmsh.model.occ.addPoint(length_member / 2 - embedded_length_longitudinal_rebar / 2, width_member / 2, cover_spacing_longitudinal_rebar + diameter_longitudinal_rebar / 2, tag=tag_point_count_longitudinal_rebar)
gmsh.model.occ.addPoint(length_member / 2 + embedded_length_longitudinal_rebar / 2, width_member / 2, cover_spacing_longitudinal_rebar + diameter_longitudinal_rebar / 2, tag=tag_point_count_longitudinal_rebar+1)
domain_longitudinal_rebar = gmsh.model.occ.addLine(tag_point_count_longitudinal_rebar, tag_point_count_longitudinal_rebar+1, tag=tag_line_longitudinal_rebar)
gmsh.model.occ.synchronize()
##### embed straight line #####
gmsh.model.mesh.embed(1, [tag_line_longitudinal_rebar], 3, tag_volume_concrete)
gmsh.model.occ.synchronize()
gmsh.model.occ.removeAllDuplicates()
volume_3D = gmsh.model.getEntities(dim=3)
##### allocate marker and name #####
gmsh.model.addPhysicalGroup(3, [tag_volume_concrete], marker_volume_concrete)
gmsh.model.setPhysicalName(3, marker_volume_concrete, 'volume_concrete')
##### allocate marker and name #####
gmsh.model.addPhysicalGroup(1, [tag_line_longitudinal_rebar], marker_volume_longitudinal_rebar)
gmsh.model.setPhysicalName(1, marker_volume_longitudinal_rebar, 'line_transverse_rebar')
resolution = 25
gmsh.model.mesh.setSize(gmsh.model.getEntities(1), resolution) # 1D 선
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01 * resolution)
gmsh.option.setNumber('Mesh.MeshSizeMax', 10 * resolution)
gmsh.option.setNumber("Mesh.MshFileVersion", 4.1) # MSH 파일 버전 (4.1, 2.2 등)
gmsh.option.setNumber("Mesh.Binary", 0) # ASCII 형식 (Binary 사용 시 1)
gmsh.model.occ.synchronize()
#Generate mesh 3 : Cell, 2 : Facet, 1 : Edge
gmsh.model.mesh.generate(3)
#save model
gmsh_model=gmsh.model
#save msh file
gmsh.model.occ.synchronize()
gmsh.fltk.run()
gmsh.write('mesh3D.msh')
gmsh.finalize()
# gmsh.model.occ.synchronize()
# gmsh.fltk.run()
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
"name_to_read": [cell_data]})
return out_mesh
# read msh file
mesh1 = meshio.read("mesh3D.msh")
mesh2=create_mesh(mesh1,"tetra")
# Save 3D mesh
meshio.write("mesh3D.xdmf" , mesh2)
# Save 1D mesh
mesh3=create_mesh(mesh1,"line")
meshio.write("mesh1D.xdmf", mesh3)
with XDMFFile(MPI.COMM_WORLD, "mesh3D.xdmf", "r") as xdmf:
domain_3D = xdmf.read_mesh(name="Grid")
#Save cell tag
cell_marker_1 = xdmf.read_meshtags(domain_3D, name="Grid")
xdmf.close()
with XDMFFile(MPI.COMM_WORLD, "mesh1D.xdmf", "r") as xdmf:
domain_1D = xdmf.read_mesh(name="Grid")
#Save cell tag
# edge_marker_2 = xdmf.read_meshtags(domain_1D, name="Grid")
xdmf.close()
########################solve#######################
# Scaled variable
mu = 1
rho = 1
delta = width_member / length_member
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma
stiffness_axial=1
V = fem.functionspace(domain_3D, ("Lagrange", 1, (3, ))) #tetra element calculate stress and strain at quadrature points
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = domain_3D.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain_3D, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
T = fem.Constant(domain_3D, default_scalar_type((0, 0, 0)))
ds_3D = ufl.Measure("ds", domain=domain_3D)
dx_3D = ufl.Measure("dx", domain=domain_3D)
dx_1D = ufl.Measure("dx", domain=domain_1D)
def epsilon_3D(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma_3D(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon_3D(u)
def tangent_calculation(domain_1D):
tangent_matrix=ufl.Jacobian(domain_1D) # shape is (3,1)
return ufl.as_vector([tangent_matrix[0,0], tangent_matrix[1, 0], tangent_matrix[2, 0]])/ufl.sqrt(ufl.inner(tangent_matrix,tangent_matrix))
def epsilon_1D(delta_displacement_Vector, domain_1D):
tangent_vector=tangent_calculation(domain_1D)
axial_strain=ufl.dot(ufl.dot(ufl.grad(delta_displacement_Vector),
tangent_vector),tangent_vector)
return axial_strain
def sigma_1D(u,domain_1D):
return stiffness_axial* epsilon_1D(u,domain_1D)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain_3D, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma_3D(u), epsilon_3D(v)) * dx_3D
a += sigma_1D(u,domain_1D)*epsilon_1D(v,domain_1D)* dx_1D
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds_3D
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()