Linear Elastic Analysis of a Beam with 1D Line Element Embedded in a 3D Tetra Elements

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()

3D-1D embedded grids are currently not supported in DOLFINx.
I’ve made two pull requests that would be a starting point for this at:

1 Like

Thank you, Dokken.
Once this feature is implemented, I believe users like me will be more actively drawn to Dolfinx.
I truly appreciate your efforts in addressing various inquiries. Thank you for your hard work!

Dear Dokken,

As shown in the attached code above, if the node positions of the 3D mesh and 1D mesh are the same, would it be possible to perform the following analysis in the current Dolfinx?

  • Establish a connection between the nodes of the 3D mesh and the 1D mesh where the line is located.
  • Impose constraints so that the nodes share the same displacement.

If there is any related tutorial, I would greatly appreciate it as a reference.

Thank you so much for your invaluable help!

This should be possible.
There are some examples of this at DOLFINX_MPC:

It would need some adaptation, as I assume your 1D lines are embedded in the interior of the mesh.

Will the whole skeleton made by 1D lines have the same displacement?

1 Like

Thank you for taking an interest in this issue!

In response to your inquiry, I have prepared a drawing as shown below.

The black lines represent the 3D mesh, and the red lines represent the 1D mesh.
The 1D mesh shares the same coordinates with the 3D mesh, and the nodes at these shared coordinates have the same displacements. However, not all nodes within the 1D mesh necessarily have identical displacements.

Embedding the mesh would significantly reduce computational cost, but if that is not feasible, implementing the constraints based on this concept would still be incredibly useful. (In Abaqus or LS-DYNA, this concept is referred to as “tie constraints.” I have emailed you an excerpt from a related document in case you might find it helpful.)

It seems that similar inquiries have been raised by others, likely researchers or engineers from the structural field like myself. I deeply appreciate your assistance, and once this feature is completed and my research is finalized, I will actively promote Dolfinx at conferences.
(The automatic differentiation feature is absolutely remarkable. Automating complex tensor differentiation and simplifying coding is tremendously beneficial for researchers.)

Dr. Dokken,

Thank you for your dedication to advancing both academic and industrial fields.

Attached below is the material regarding the tie constraints mentioned in my previous inquiry.

I sincerely appreciate your invaluable assistance.

그림2.jpg

Wishing you a prosperous and happy New Year!

Best regards,

2025년 1월 15일 (수) 오전 1:59, Jorgen Dokken via FEniCS Project <notifications@fenicsproject1.discoursemail.com>님이 작성:

(Attachment 그림2.jpg is missing)