Converting legacy FEniCS code to FEniCSx causes issue

[mohaiminul_islam_Mis]

Hi
I was trying to convert the 3D beam found here FEniCSx.
I am facing 2 issues.

1. I was able to convert the code. But if i change the loading condition to concentrated load FEniCSx gives no solution for complicated mesh (For simple mesh works fine).
2. the .xdmf file written can’t be opened by paraview for Meshtags however opening mesh works fine. And if i try writing the write_function my Jupyter notebook crashes and restarts.

For the First problem, If i use simple mesh everything works (concentrated load or self-weight). If i use a little complicated mesh self-weight problem works fine but concentrated one gives no solution.
I installed FEniCSx with conda and dolfinx version is 0.7.3
My MWE

# static 1D Cantilever Beam with a rectangular cross section
import matplotlib.pyplot as plt
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
                         locate_dofs_topological,Constant,Expression, locate_dofs_geometrical, assemble_scalar)
from dolfinx.io import VTKFile, gmshio, XDMFFile, VTXWriter, VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary,create_interval ,locate_entities, meshtags
from dolfinx import mesh
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, dS, ds, Jacobian, as_vector, sqrt, cross, dot, avg,
diag,VectorElement, split, Measure)
import meshio
import gmsh
## Complicated mesh 
msh = meshio.read("../fenics_nlopt/examples/Truss/gmsh/t3/test.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points = msh.points, cells={"line": msh.cells_dict['line']}))
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

## Simple mesh 
#length = 10.0
#N = 40 # number of elements

#points = np.zeros((N+1, 3))
#points[:, 0] = np.linspace(0, length, N+1)
#cells = np.array([[i, i+1] for i in range(N)])
#meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)

def tangent(domain):
    t = Jacobian(domain)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

t = tangent(domain)

ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
thick = Constant(domain,0.3)
width = thick/3
E = Constant(domain,70e3)
nu = Constant(domain,0.3)
G = E/2/(1+nu)
rho = Constant(domain,2.7e-3)
g = Constant(domain,9.81)

S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = Constant(domain,5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S

Ue = VectorElement("CG", domain.ufl_cell(), 1, dim=3)
W = FunctionSpace(domain, Ue*Ue)

u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

def tgrad(u):
    return dot(grad(u), t)
def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([dot(tgrad(w), t),
                      dot(tgrad(w), a1)-dot(theta, a2),
                      dot(tgrad(w), a2)+dot(theta, a1),
                      dot(tgrad(theta), t),
                      dot(tgrad(theta), a1),
                      dot(tgrad(theta), a2)])
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))

Sig = generalized_stresses(du)
Eps =  generalized_strains(u_)
def left(x):
    return(np.isclose(x[0], 0, 1))

def load(x):
    return np.isclose(x[0], 10, 1) #and np.isclose(x[1], 0, 12)
    
boundaries = [(1, left),
              (2, load)]

    
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])


W0, submap0 = W.sub(0).collapse()
W1, submap1 = W.sub(1).collapse()
left_dofs_x = locate_dofs_geometrical((W.sub(0),W0),left )
left_dofs_r = locate_dofs_geometrical((W.sub(1),W1),left)
#right_dofs = locate_dofs_geometrical(W, right)
zero = Function(W0)
with zero.vector.localForm() as zero_local:
    zero_local.set(0.0)

bcs = [
    dirichletbc(zero, left_dofs_x, W.sub(0)),
    dirichletbc(zero, left_dofs_r, W.sub(1))]
ds = Measure("ds", domain=domain, subdomain_data=facet_tag, metadata = {"quadrature_degree": 4})
dS = Measure("dS", domain=domain, subdomain_data=facet_tag,metadata = {"quadrature_degree": 4})
#dx = Measure("dx",domain=domain, subdomain_data=subdomains, metadata = {"quadrature_degree": 4})
dx_shear = dx(scheme="default",metadata={ "quadrature_degree": 1})

k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
# 
l_form = inner(Constant(domain,-1.), 2*avg(u_[1]))*dS(2)#Works for simple mesh not for complicated one 
#l_form = inner(Constant(domain,-1.), u_[1])*dx# Works for both
u = Function(W)
problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()

assemble_scalar((form(inner(uh,uh)*dx)))

with XDMFFile(domain.comm, "results/tag.xdmf", "w" ,encoding=XDMFFile.Encoding.HDF5) as xdmf:
    xdmf.write_meshtags(facet_tag, domain.geometry)
with XDMFFile(domain.comm, "results/fun.xdmf", "w" ) as xdmf:
    xdmf.write_function(uh)

the test.msh is below

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
6 11 0 0
1 0 0 0 0 
2 10 0 0 0 
3 10 10 0 0 
4 0 10 0 0 
6 20 0 0 0 
10 20 10 0 0 
1 0 0 0 10 0 0 0 2 1 -2 
2 10 0 0 10 10 0 0 2 2 -3 
3 0 10 0 10 10 0 0 2 3 -4 
4 0 0 0 0 10 0 0 2 4 -1 
5 0 0 0 10 10 0 0 2 2 -4 
6 0 0 0 10 10 0 0 2 3 -1 
7 10 0 0 20 0 0 0 2 2 -6 
8 20 0 0 20 10 0 0 2 6 -10 
9 10 10 0 20 10 0 0 2 10 -3 
11 10 0 0 20 10 0 0 2 6 -3 
12 10 0 0 20 10 0 0 2 10 -2 
$EndEntities
$Nodes
17 6 1 6
0 1 0 1
1
0 0 0
0 2 0 1
2
10 0 0
0 3 0 1
3
10 10 0
0 4 0 1
4
0 10 0
0 6 0 1
5
20 0 0
0 10 0 1
6
20 10 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 5 0 0
1 6 0 0
1 7 0 0
1 8 0 0
1 9 0 0
1 11 0 0
1 12 0 0
$EndNodes
$Elements
17 17 1 17
0 1 15 1
1 1 
0 2 15 1
2 2 
0 3 15 1
3 3 
0 4 15 1
4 4 
0 6 15 1
5 5 
0 10 15 1
6 6 
1 1 1 1
7 1 2 
1 2 1 1
8 2 3 
1 3 1 1
9 3 4 
1 4 1 1
10 4 1 
1 5 1 1
11 2 4 
1 6 1 1
12 3 1 
1 7 1 1
13 2 5 
1 8 1 1
14 5 6 
1 9 1 1
15 6 3 
1 11 1 1
16 5 3 
1 12 1 1
17 6 2 
$EndElements

[dokken]

@mohaiminul_islam_Mis

the .xdmf file written can’t be opened by paraview for Meshtags however opening mesh works fine. And if i try writing the write_function my Jupyter notebook crashes and restarts.

For both of these cases you need to store the mesh to the XDMFFile (using '.write_mesh`). This is shown in all demos and tutorials that uses XDMFFile.

[dokken]

@mohaiminul_islam_Mis

zero = Function(W0)
with zero.vector.localForm() as zero_local:
    zero_local.set(0.0)

bcs = [
    dirichletbc(zero, left_dofs_x, W.sub(0)),
    dirichletbc(zero, left_dofs_r, W.sub(1))]

Here the second DirichletBC should use a zero function from W1, not W0.

I would also like to point out that the demo you are linking to solves what Seems to me to be a very different problem (making it hard to do any sensible comparison). Maybe @bleyerj can comment:

[mohaiminul_islam_Mis]

Apologies, i linked the wrong example link.
Here is corrected link

[bleyerj]

The new version of this demo should be available soon. For concentrated loadings, please be sure that you correctly identify and mark inner “facets” which are points in fact.
To avoid any issue with XDMFFile on manifold meshes, use gmsh Python API instead

[mohaiminul_islam_Mis]

Great. looking forward to the new examples.

[mohaiminul_islam_Mis]

Thanks for the reply.
There are couple of problems,
Now i can write the meshtags, But when i open the .xdmf file to paraview, i can’t change the color or i don’t see the change in color for marked points.
As for the write_function My Jupyter notebook still crashes for the below code if i change it in the current problem.

with XDMFFile(domain.comm, "results/fun.xdmf", "w" ) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)

However, if i create a simple cubic domain for test, it works fine.

domain = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
V = FunctionSpace(domain, ("Lagrange", 1))
u = Function(V)

with XDMFFile(domain.comm, "results/test.xdmf", "w") as xdmf:
  xdmf.write_mesh(domain)
  xdmf.write_function(u)

[dokken]

@mohaiminul_islam_Mis

There are couple of problems,
Now i can write the meshtags, But when i open the .xdmf file to paraview, i can’t change the color or i don’t see the change in color for marked points.

Use the extract block filter to extract the meshtags from the multi-block system.

Could you please create a minimal example with your other mesh, where you simply load the mesh, create a function space and function (no need to solve any problem) on it and store it to file.
If that doesn’t work, please post that code below.

[mohaiminul_islam_Mis]

Here is a simplified version where I load a simple mesh. define the function spaces and mark the domain.

My Jupyter notebook still crashes when I write function.

import matplotlib.pyplot as plt
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
                         locate_dofs_topological,Constant,Expression, locate_dofs_geometrical, assemble_scalar)
from dolfinx.io import VTKFile, gmshio, XDMFFile, VTXWriter, VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary,create_interval ,locate_entities, meshtags
from dolfinx import mesh
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, dS, ds, Jacobian, as_vector, sqrt, cross, dot, avg,
diag,VectorElement, split, Measure)
import meshio
import gmsh


# Simple mesh 
length = 10.0
N = 40 # number of elements

points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

Ue = VectorElement("CG", domain.ufl_cell(), 1, dim=3)
W = FunctionSpace(domain, Ue*Ue)

def left(x):
    return(np.isclose(x[0], 0, 1))

def load(x):
    return np.isclose(x[0], 20, 1) #and np.isclose(x[1], 0, 12)
    
boundaries = [(1, left),
              (2, load)]

    
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

u = Function(W)


domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "results/tag.xdmf", "w" ) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)
    
with XDMFFile(domain.comm, "results/fun.xdmf", "w" ) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u)

[dokken]

Try the following:

import matplotlib.pyplot as plt
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
                         locate_dofs_topological,Constant,Expression, locate_dofs_geometrical, assemble_scalar)
from dolfinx.io import VTKFile, gmshio, XDMFFile, VTXWriter, VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary,create_interval ,locate_entities, meshtags
from dolfinx import mesh
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, dS, ds, Jacobian, as_vector, sqrt, cross, dot, avg,
diag,VectorElement, split, Measure)
import meshio
import gmsh


# Simple mesh 
length = 10.0
N = 40 # number of elements

points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

Ue = VectorElement("CG", domain.ufl_cell(), 1, dim=3)
W = FunctionSpace(domain, Ue*Ue)

def left(x):
    return(np.isclose(x[0], 0, 1))

def load(x):
    return np.isclose(x[0], 20, 1) #and np.isclose(x[1], 0, 12)
    
boundaries = [(1, left),
              (2, load)]

    
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

u = Function(W)

u0 = u.sub(0).collapse()
u0.name = "u0"
u1 = u.sub(1).collapse()
u1.name = "u1"
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "results/tag.xdmf", "w" ) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)
    
with XDMFFile(domain.comm, "results/fun.xdmf", "w" ) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u0)
    xdmf.write_function(u1)

[mohaiminul_islam_Mis]

Thanks. now it works.
However, i still can’t see colors for the marked points for the facet tags.

[dokken]

Open the file in paraview, choose selecy block filter and extract the second block. Then you need to choose to visualize the colors (where it says solid domain)
See for instance How to write binary data with VTK/XDMFFile? - #9 by dokken

[mohaiminul_islam_Mis]

I am still facing a issue where the triangle gives same result as the legacy dolfin. and the rectangle doesn’t give any solution. I checked the marked points, they seem correct.

MWE

import matplotlib.pyplot as plt
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
                         locate_dofs_topological,Constant,Expression, locate_dofs_geometrical, assemble_scalar)
from dolfinx.io import VTKFile, gmshio, XDMFFile, VTXWriter, VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary,create_interval ,locate_entities, meshtags
from dolfinx import mesh
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, dS, ds, Jacobian, as_vector, sqrt, cross, dot, avg,
diag,VectorElement, split, Measure)
import meshio
import gmsh
## Complicated mesh 
msh = meshio.read("../fenics_nlopt/examples/Truss/gmsh/t2/x-section.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points = msh.points, cells={"line": msh.cells_dict['line']}))
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

## Simple mesh 
#length = 10.0
#N = 40 # number of elements

#points = np.zeros((N+1, 3))
#points[:, 0] = np.linspace(0, length, N+1)
#cells = np.array([[i, i+1] for i in range(N)])
#meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)

def tangent(domain):
    t = Jacobian(domain)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

t = tangent(domain)

ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
thick = Constant(domain,0.3)
width = thick/3
E = Constant(domain,70e3)
nu = Constant(domain,0.3)
G = E/2/(1+nu)
rho = Constant(domain,2.7e-3)
g = Constant(domain,9.81)

S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = Constant(domain,5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S

Ue = VectorElement("CG", domain.ufl_cell(), 1, dim=3)
W = FunctionSpace(domain, Ue*Ue)

u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

def tgrad(u):
    return dot(grad(u), t)
def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([dot(tgrad(w), t),
                      dot(tgrad(w), a1)-dot(theta, a2),
                      dot(tgrad(w), a2)+dot(theta, a1),
                      dot(tgrad(theta), t),
                      dot(tgrad(theta), a1),
                      dot(tgrad(theta), a2)])
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))

Sig = generalized_stresses(du)
Eps =  generalized_strains(u_)
def left(x):
    return(np.isclose(x[0], 0, 1))

def load(x):
    return np.isclose(x[0], 10, 1) #and np.isclose(x[1], 0, 12)
    
boundaries = [(1, left),
              (2, load)]

    
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])



W0, submap0 = W.sub(0).collapse()
W1, submap1 = W.sub(1).collapse()
left_dofs_x = locate_dofs_geometrical((W.sub(0),W0),left )
left_dofs_r = locate_dofs_geometrical((W.sub(1),W1),left)
#right_dofs = locate_dofs_geometrical(W, right)
zero1 = Function(W0)
with zero1.vector.localForm() as zero1_local:
    zero1_local.set(0.0)
    
zero2 = Function(W1)
with zero2.vector.localForm() as zero2_local:
    zero2_local.set(0.0)

bcs = [
    dirichletbc(zero1, left_dofs_x, W.sub(0)),
    dirichletbc(zero2, left_dofs_r, W.sub(1))]
ds = Measure("ds", domain=domain, subdomain_data=facet_tag, metadata = {"quadrature_degree": 4})
dS = Measure("dS", domain=domain, subdomain_data=facet_tag,metadata = {"quadrature_degree": 4})
#dx = Measure("dx",domain=domain, subdomain_data=subdomains, metadata = {"quadrature_degree": 4})
dx_shear = dx(scheme="default",metadata={ "quadrature_degree": 1})

k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
# 
l_form = inner(Constant(domain,-1.), 2*avg(u_[1]))*dS(2)#Works for simple mesh not for complicated one 
#l_form = inner(Constant(domain,-1.), u_[1])*dx# Works for both
u = Function(W)
problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()

assemble_scalar((form(inner(uh,uh)*dx)))

rectangular mesh rect.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
4 6 0 0
1 0 0 0 0 
2 10 0 0 0 
3 10 10 0 0 
4 0 10 0 0 
1 0 0 0 10 0 0 0 2 1 -2 
2 10 0 0 10 10 0 0 2 2 -3 
3 0 10 0 10 10 0 0 2 3 -4 
4 0 0 0 0 10 0 0 2 1 -4 
5 0 0 0 10 10 0 0 2 1 -3 
6 0 0 0 10 10 0 0 2 2 -4 
$EndEntities
$Nodes
10 4 1 4
0 1 0 1
1
0 0 0
0 2 0 1
2
10 0 0
0 3 0 1
3
10 10 0
0 4 0 1
4
0 10 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 5 0 0
1 6 0 0
$EndNodes
$Elements
10 10 1 10
0 1 15 1
1 1 
0 2 15 1
2 2 
0 3 15 1
3 3 
0 4 15 1
4 4 
1 1 1 1
5 1 2 
1 2 1 1
6 2 3 
1 3 1 1
7 3 4 
1 4 1 1
8 1 4 
1 5 1 1
9 1 3 
1 6 1 1
10 2 4 
$EndElements

triangular mesh tri.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
3 3 0 0
1 0 0 0 0 
2 10 0 0 0 
3 0 10 0 0 
1 0 0 0 10 0 0 0 2 1 -2 
2 0 0 0 10 10 0 0 2 2 -3 
3 0 0 0 0 10 0 0 2 3 -1 
$EndEntities
$Nodes
6 3 1 3
0 1 0 1
1
0 0 0
0 2 0 1
2
10 0 0
0 3 0 1
3
0 10 0
1 1 0 0
1 2 0 0
1 3 0 0
$EndNodes
$Elements
6 6 1 6
0 1 15 1
1 1 
0 2 15 1
2 2 
0 3 15 1
3 3 
1 1 1 1
4 1 2 
1 2 1 1
5 2 3 
1 3 1 1
6 3 1 
$EndElements

[JBo]

Good evening
I wish I was able to answer your question, but I have another one instead.

In the code above you use

Ue = VectorElement("CG", domain.ufl_cell(), 1, dim=3)

it looks like ufl.VectorElement has been replaced with basix.ufl.element)

So I tried the following with your simple beam mesh

# Define the vector element using basix
Ue = basix.ufl.element(
    family=basix.ElementFamily.P,  # Lagrange family
    cell=basix.CellType.interval,  # 1D cell type
    degree=1, 
    shape=(3,))

but I get this error I don’t manage :

---> 13 W = FunctionSpace(domain, Ue*Ue)
unsupported operand type(s) for *: '_BlockedElement' and '_BlockedElement'

Also tried :

Ue = dolfinx.fem.functionspace(domain, ("CG", 2,(gdim,)))
W = FunctionSpace(domain, Ue*Ue)

but got :

----> 2 W = FunctionSpace(domain, Ue*Ue)
TypeError: unsupported operand type(s) for *: 'FunctionSpace' and 'FunctionSpace'

Can anyone help me through this one ?

[francesco-ballarin]

Use Ueprod = basix.ufl.mixed_element([Ue, Uel]), and then build a dolfinx.fem.functionspace on Ueprod. See for instance Mixed formulation for the Poisson equation — DOLFINx 0.8.0 documentation.

[JBo]

Thanks a lot
It works with :

# Define the vector element using basix
Ue = basix.ufl.element(
    family=basix.ElementFamily.P,  # Lagrange family
    cell=basix.CellType.interval,  # 1D cell type
    degree=1,
    shape=(3,))

Ueprod = basix.ufl.mixed_element([Ue, Ue])
W = dolfinx.fem.functionspace(domain, Ueprod)

[JBo]

Now I am trying to print the maximum y-displacement to compare with the predictions of beam theory. I’ve read a few posts, including the Deflection of a membrane - Implementation tutorial, but I haven’t been able to find a code that works.

Could you please give me a hint ?

Here’s what I’ve tried so far:

# Print the maximum y-displacement
# print("Maximum y-displacement =", np.max(u.sub(1).split()[0])
print("Maximum y-displacement =",np.max(-uh.sub(1).x.array))
# print("Maximum y-displacement =",-uh(length,0,0.)[1])

[francesco-ballarin]

u.sub(1).collapse().x.array

may be helpful.

[JBo]

Thanks, indeed it is !

If I understand it right …

# Get the array of displacement values
# u.sub(0) is displacement while u.sub(1) is rotation field
displacement_values = u.sub(0).collapse().x.array       # stacked 11 nodes x,y,z

# Extract y-displacement components
y_displacement = displacement_values[1::3]              # start at 1, steps of 3

# Print y-displacement components
print("y-displacement components:", y_displacement)

then I have another problem !

Max y-displacement = -1238.7051428744687
Beam theory = -0.44091710758377434

I’ll keep on digging then.
My full version of the initial code here below.

import os
arch = os.getenv("ARGS", "real")

try:
    import google.colab  # noqa: F401
except ImportError:
    import ufl
    import dolfinx
else:
    try:
        import ufl
        import dolfinx
    except ImportError:
        if arch != "complex":
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-real.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        else:
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-complex.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        import ufl
        import dolfinx

# Install meshio for mesh handling
!pip install meshio

# Importing necessary libraries and modules
import numpy as np       # Numerical operations
from mpi4py import MPI   # Parallel computing
import basix             # Finite element basis functions
from dolfinx import plot
from dolfinx.fem import (Constant, FunctionSpace, locate_dofs_geometrical,
                         Function, dirichletbc, assemble_scalar, form)
                         # FEM module functions
from dolfinx.io import XDMFFile  # I/O operations
from dolfinx.fem.petsc import LinearProblem
                         # Solving linear problems
from dolfinx.mesh import locate_entities, meshtags
                         # Mesh utility functions
from ufl import (Jacobian, as_vector, sqrt, inner, cross, dot, TestFunction,
                TrialFunction, split, diag, grad, Measure, dx, avg)
                         # Defining variational forms
import meshio            # Mesh I/O operations
import matplotlib.pyplot as plt
import pyvista as pv     # Visualization

# Initialize PyVista
pv.start_xvfb()

## Simple mesh
length = 10.0
N = 10 # number of elements

points = np.zeros((N+1, 3))                     # N+1 array of zeros in 3D space
points[:, 0] = np.linspace(0, length, N+1)      # Distributing points along the x-axis from 0 to 'length'
cells = np.array([[i, i+1] for i in range(N)])  # Defining the connectivity of the mesh (line elements)
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))  # Writing the mesh to an XDMF file

# Reading the mesh from the XDMF file
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

# Creating connectivity information for the mesh
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)

# Explicitly compute the entity-to-cell connectivity
domain.topology.create_connectivity(domain.topology.dim, 0)

# mesh test

# Get the geometrical and topological dimensions of the mesh from the 'domain' object.
gdim = domain.geometry.dim    # = 3
tdim = domain.topology.dim    # = 1

# Print the geometrical and topological dimensions of the mesh.
print(f"Geometrical dimension = {gdim}")
print(f"Topological dimension = {tdim}")

# Extract coordinates of vertices
vertices_coordinates = domain.geometry.x

# Create a 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot vertices
ax.scatter(vertices_coordinates[:, 0],
           vertices_coordinates[:, 1],
           vertices_coordinates[:, 2],
           c='b',
           marker='o')

# Set labels
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# Show plot
plt.show()

# Defining a function to compute the tangent vector to the mesh
def tangent(domain):
    t = Jacobian(domain)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

# Calculating the tangent vector for the domain
t = tangent(domain)

# Defining basis vectors orthogonal to the tangent vector
ez = as_vector([0, 0, 1])   # Basis vector in the z-direction
a1 = cross(t, ez)           # First orthogonal vector
a1 /= sqrt(dot(a1, a1))     # Normalizing the first orthogonal vector
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

# Geometric properties
thick_val = 0.3
thick = Constant(domain,thick_val)
width_val = thick_val
width = thick/3

# Material Properties
E_val = 70e3
E = Constant(domain,E_val)      # Young's modulus
nu = Constant(domain,0.3)       # Poisson's ratio
G = E/2/(1+nu)                  # Shear modulus
rho = Constant(domain,2.7e-3)   # Density
g = Constant(domain,9.81)       # Gravitational acceleration

S = thick*width                 # Cross-sectional area
ES = E*S                        # Axial stiffness
EI1 = E*width*thick**3/12       # Bending stiffness about one axis
EI2 = E*width**3*thick/12       # Bending stiffness about the other axis
GJ = G*0.26*thick*width**3      # Torsional stiffness
kappa = Constant(domain,5./6.)  # Shear correction factor
GS1 = kappa*G*S                 # Shear stiffness in one direction
GS2 = kappa*G*S                 # Shear stiffness in the other direction

# Define the vector element using basix
Ue = basix.ufl.element(
    family=basix.ElementFamily.P,  # Lagrange family
    cell=basix.CellType.interval,  # 1D cell type
    degree=1,                      # Polynomial degree
    shape=(3,))                    # Shape of the element (vector with 3 components)

# Creating a mixed element for displacement and rotation
Ueprod = basix.ufl.mixed_element([Ue, Ue])    # Mixed element
W = dolfinx.fem.functionspace(domain, Ueprod) # Creating a function space

# Defining test and trial functions
u_ = TestFunction(W)       # Test function
du = TrialFunction(W)      # Trial function
(w_, theta_) = split(u_)   # Split test function into displacement and rotation
(dw, dtheta) = split(du)   # Split trial function into displacement and rotation

# Defining a function to compute generalized strains
def tgrad(u):
    return dot(grad(u), t)

def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([dot(tgrad(w), t),
                      dot(tgrad(w), a1)-dot(theta, a2),
                      dot(tgrad(w), a2)+dot(theta, a1),
                      dot(tgrad(theta), t),
                      dot(tgrad(theta), a1),
                      dot(tgrad(theta), a2)])

# Defining a function to compute generalized stresses
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
               generalized_strains(u))

# Calculating stresses and strains
Sig = generalized_stresses(du)   # Generalized stresses
Eps =  generalized_strains(u_)   # Generalized strains

# Locate left boundary
def left(x):
    return(np.isclose(x[0], 0, 1))

# Locate load application
def load(x):
    return np.isclose(x[0], length, 1)

# Defining boundaries and marking facets
boundaries = [(1, left),              # Marker 1 for left boundary
              (2, load)]              # Marker 2 for load application point

facet_indices, facet_markers = [], [] # Lists to store facet indices and markers
fdim = domain.topology.dim - 1        # Facet dimension

# Locating and marking facets for each boundary
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)      # Locate facets
    facet_indices.append(facets)                         # Add facet indices
    facet_markers.append(np.full_like(facets, marker))   # Add facet markers

facet_indices = np.hstack(facet_indices).astype(np.int32)# Combine facet indices
facet_markers = np.hstack(facet_markers).astype(np.int32)# Combine facet markers

sorted_facets = np.argsort(facet_indices)                # Sort facets

facet_tag = meshtags(domain,
                     fdim,
                     facet_indices[sorted_facets],
                     facet_markers[sorted_facets])       # Create mesh tags

# Collapsing subspaces of the function space
W0, submap0 = W.sub(0).collapse()   # Collapsing displacement field subspace
W1, submap1 = W.sub(1).collapse()   # Collapsing rotation field subspace

left_dofs_x = locate_dofs_geometrical((W.sub(0),W0),left )
left_dofs_r = locate_dofs_geometrical((W.sub(1),W1),left)

# Create zero function for the first subspace (W0)
zero_0 = Function(W0)
with zero_0.vector.localForm() as zero_local:
    zero_local.set(0.0)

# Create zero function for the second subspace (W1)
zero_1 = Function(W1)
with zero_1.vector.localForm() as zero_local:
    zero_local.set(0.0)

# Define boundary conditions
bcs = [
    dirichletbc(zero_0, left_dofs_x, W.sub(0)),
    dirichletbc(zero_1, left_dofs_r, W.sub(1))]

# Define a measure for surface integrals
ds = Measure("ds",
             domain=domain,
             subdomain_data=facet_tag,
             metadata = {"quadrature_degree": 4})

# Define a measure for surface integrals in 3D
dS = Measure("dS",
             domain=domain,
             subdomain_data=facet_tag,metadata = {"quadrature_degree": 4})

# Define a measure for volume integrals, used for shear stress
dx_shear = dx(scheme="default",metadata={ "quadrature_degree": 1})

# Define the bilinear form (stiffness matrix)
k_form = (sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]])
             +(Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear)

P = -1.0   # load intensity

# Define the linear form (load vector)
l_form = inner(Constant(domain,P),
               2*avg(u_[1]))*dS(2)          # u_[1] means y , dS(2) means mark 2

u = Function(W)                             # Define the solution function

# Define linear problem
problem = LinearProblem(k_form,
                        l_form,
                        u=u,
                        bcs=bcs)
uh=problem.solve()                          # Solve the linear problem

# Get the array of displacement values
displacement_values = u.sub(0).collapse().x.array       # stacked 11 nodes x,y,z

# Extract y-displacement components
y_displacement = displacement_values[1::3]              # start at 1, steps of 3

# Print the maximum y-displacement
print("Max y-displacement =",min(y_displacement)) # outputs > ? = [0. 0. 0. 0. 0. -46.28571429 0. 0. 

# Calculate maximum y-displacement using beam theory
w_max = (P * length**3) / (48 * E_val*width_val*thick_val**3/12)
print("Beam theory =", w_max)

[francesco-ballarin]

@JBo

# Get the array of displacement values
# u.sub(0) is displacement while u.sub(1) is rotation field
displacement_values = u.sub(0).collapse().x.array       # stacked 11 nodes x,y,z

In that case you want u.sub(0).sub(1).collapse().x.array, where sub(0) gets the displacement field, and sub(1) gets its y component.