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

Hoseong_Jeong · January 14, 2025, 8:19am

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

dokken · January 14, 2025, 8:58am

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:

Hoseong_Jeong · January 14, 2025, 9:15am

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!

Hoseong_Jeong · January 14, 2025, 12:12pm

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!

dokken · January 14, 2025, 4:49pm

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

github.com/jorgensd/dolfinx_mpc

Rigidity constraint

opened 03:33PM - 11 Dec 21 UTC

closed 02:23PM - 15 Feb 22 UTC

cmaurini

Thanks a lot for the work on this library, which is really useful. I am tryi…ng to use dolfinx_mpc to impose a boundary to have the same displacement (something like [this](https://fenicsproject.discourse.group/t/rigid-plate-boundary-condition/2056 )) In particular for a linear elasticity problem with a vector-valued displacement field I want only to have "rigidity" on one of the components. I managed to obtain the results with the following method. However, this is very slow and surely not the good way of doing it. Could you suggest me how to proceed? ``` def create_rigid_constraint(mpc, dofs, subspace=None): dof_coordinates = mpc.V.tabulate_dof_coordinates()[dofs, :] for i in range(len(dof_coordinates) - 1): mpc.create_general_constraint( {dof_coordinates[i + 1, :].tobytes(): {dof_coordinates[0, :].tobytes(): 1}}, subspace, subspace, ) right_facets = locate_entities_boundary(mesh, 1, right_boundary) rigid_dofs = fem.locate_dofs_topological(V, 1, right_facets) mpc = dolfinx_mpc.MultiPointConstraint(V) create_rigid_constraint(mpc, rigid_dofs, 0) ``` Besides, when considering periodic boundary conditions could you consider adding an option to impose it only on one of the components of a vector-valued field. This would be also a possible solution to the problem above, because "rigidity" can easily be written with a particular periodicity map. I think that this "rigidity" constraint is very useful in many applications as shown also by the several questions on the forums about that. Maybe it could be added between the built-in methods off `MultiPointConstraint`? For completeness, I give below also a full code for a simple example to play with ``` import dolfinx.fem as fem import dolfinx_mpc.utils import numpy as np from dolfinx.common import Timer, list_timings, TimingType from dolfinx.generation import UnitSquareMesh from dolfinx.io import XDMFFile from dolfinx.mesh import locate_entities_boundary from mpi4py import MPI from petsc4py import PETSc from ufl import (Identity, SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx, grad, inner, sym, tr, ) mesh = UnitSquareMesh(MPI.COMM_WORLD, 50, 50) V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) # Generate Dirichlet BC on lower boundary (Fixed) u_bc = fem.Function(V) with u_bc.vector.localForm() as u_local: u_local.set(0.0) def left_boundary(x): return np.isclose(x[0], np.finfo(float).eps) left_facets = locate_entities_boundary(mesh, 1, left_boundary) topological_dofs = fem.locate_dofs_topological(V, 1, left_facets) bc = fem.DirichletBC(u_bc, topological_dofs) bcs = [bc] # Define variational problem u = TrialFunction(V) v = TestFunction(V) # Elasticity parameters E = PETSc.ScalarType(1.0) nu = 0.3 mu = fem.Constant(mesh, E / (2.0 * (1.0 + nu))) lmbda = fem.Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) # Stress computation def sigma(v): return 2.0 * mu * sym(grad(v)) + lmbda * tr(sym(grad(v))) * Identity(len(v)) x = SpatialCoordinate(mesh) # Define variational problem u = TrialFunction(V) v = TestFunction(V) a_form = inner(sigma(u), grad(v)) * dx L_form = inner(as_vector((10 * (x[1] - 0.5) ** 2, 0)), v) * dx # Create MPC for imposing u[1] constant on the bottom face def right_boundary(x): return np.isclose(x[0], 1.0) right_facets = locate_entities_boundary(mesh, 1, right_boundary) rigid_dofs = fem.locate_dofs_topological(V, 1, right_facets) def create_rigid_constraint(mpc, dofs, subspace=None): dof_coordinates = mpc.V.tabulate_dof_coordinates()[dofs, :] for i in range(len(dof_coordinates) - 1): mpc.create_general_constraint( {dof_coordinates[i + 1, :].tobytes(): {dof_coordinates[0, :].tobytes(): 1}}, subspace, subspace, ) with Timer("~RIGID: Initialize MPC"): mpc = dolfinx_mpc.MultiPointConstraint(V) create_rigid_constraint(mpc, rigid_dofs, 0) mpc.finalize() problem = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs) solver = problem.solver u_h = problem.solve() u_h.name = "u" outfile = XDMFFile(MPI.COMM_WORLD, "rigid_boundary.xdmf", "w") outfile.write_mesh(mesh) outfile.write_function(u_h) list_timings(MPI.COMM_WORLD, [TimingType.wall]) ```

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?

Hoseong_Jeong · January 15, 2025, 1:32am

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

Hoseong_Jeong · January 15, 2025, 1:38am

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)

Topic		Replies	Views
Unexpected Behaviour with .sub in fenicsx dolfinx	2	391	September 30, 2022
How to define 1d element embedded in 3d space? General	0	29	October 31, 2024
Creating a 1D beam (full cylinder) clamped-free dolfinx	10	448	June 9, 2024
Error with 1d mesh immersed in 2d with dolfinx dolfinx	3	350	June 22, 2023
Unable to properly locate dofs dolfinx	3	523	March 30, 2021

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

Related topics