Hey there,
The following code generates a line mesh from (0,0) to (1,0) and solves a simple FEM-problem that can be interpreted as a simple bar of length, cross-section and Youngs Modulus = 1 where the point (0,0) is fixed and a force of magnitude 0.1 pulls on the right-hand side (at (1,0)) in positive x-direction.
import ufl
from dolfinx import io, fem, mesh
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
import numpy as np
gmsh.initialize()
meshSize = 0.05
point1 = gmsh.model.geo.add_point(0.,0.,0.,meshSize)
point2 = gmsh.model.geo.add_point(1.,0.,0.,meshSize)
edge1 = gmsh.model.geo.add_line(point1,point2)
gmsh.model.geo.synchronize()
edges = [edge1]
model1 = gmsh.model.add_physical_group(1,edges)
gmsh.model.mesh.generate(1)
domain = io.ufl_mesh_from_gmsh(1,2)
geometry_data = io.extract_gmsh_geometry(gmsh.model)
topology_data = io.extract_gmsh_topology_and_markers(gmsh.model)
cells = topology_data[1]["topology"]
nodes = geometry_data[:,:2]
struct_mesh = mesh.create_mesh(MPI.COMM_WORLD, cells, nodes, domain)
Vu = fem.VectorFunctionSpace(struct_mesh,("CG",1))
du = ufl.TrialFunction(Vu)
u_ = ufl.TestFunction(Vu)
# bilinear form
aM = ufl.inner(ufl.grad(du),ufl.grad(u_))*ufl.dx
# "fixing" node at (0,0)
b_dofs_l = fem.locate_dofs_geometrical(Vu, lambda x : (np.isclose(x[0],0.)) & (np.isclose(x[1],0.)))
bc_left = fem.dirichletbc(fem.Constant(struct_mesh,np.array([0.,0.])),b_dofs_l,Vu)
# fixing y-dof at (1,0)
b_dofs_r = fem.locate_dofs_geometrical(Vu, lambda x : (np.isclose(x[0],1.) & (np.isclose(x[1],0.))))
bc_right = fem.dirichletbc(fem.Constant(struct_mesh,0.),b_dofs_r,Vu.sub(1))
bcs = [bc_left]#, bc_right]
# assemble lhs
A = fem.petsc.assemble_matrix(fem.form(aM),bcs=bcs)
A.assemble()
# create solver
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setOperators(A)
# create empty rhs-vector
LM = ufl.inner(u_,fem.Constant(struct_mesh,np.array([0.,0.])))*ufl.dx
b = fem.petsc.create_vector(fem.form(LM))
# applying a force to the node at the very right, at (1,0)
xForce = 0.1
b.setValue(b_dofs_r[0]*2,xForce) # accessing x-component of node b_dofs_r
fem.petsc.set_bc(b,bcs=bcs)
b.assemble()
# create a function to write the solution into
u = fem.Function(Vu)
solver.solve(b,u.vector)
If I run the code like this, I get the expected result: The value of u increases linearly from 0. to .1 - exactly like the (trivial) analytical solution.
If I change the line bcs = [bc_left]#, bc_right] to bcs = [bc_left, bc_right], I expect that the boundary condition bc_right is applied too. This BC would constrain the y-dof at the end node ((1,0)). In this case this shouldn’t have an effect on the solution, but sets my engineering mind (which is concerned with static/kinematic determinacy considerations) to peace.
However, the result is not as expected at all. The u-function is constant (value 0.) between (0,0) and (0.5,0) and then increases linearly to 0.05 up to the point (1,0). It seems like instead of constraining the y-dof of (1,0), bc_right constrained the x-dof of (0.5,0).
How can this be explained?
PS: Another question: I noticed that b_dofs_r is 20 in this case. However, node 20 in the array nodes is at (0.95,0) - while node 1 (the second one) is at (1.,0.). Are the nodes reshuffled when the mesh is imported from gmsh to fenicsx?