Hi, could you please suggest me if I have to change the form of the equation of linear elasticity (Implementation — FEniCSx tutorial) if I want to solve it on curved surfaces, such as, on the surface of a sphere?
Or I just have to supply a mesh for the curved surface and FEniCSx will do the rest.
Thanks.
You need to handle the nullspace appropriately if solving on a sphere. Other than that i think it should «just work». Maybe @jackhale can comment?
You do not need to change the form of equations. You just generate the spherical mesh with a cap and define appropriate physical tags for surfaces/volumes. FEniCSx should handle the rest.
FEniCSx does not know if the surface is curved or not at all, it will assemble the stiffness matrix and load vector using coordinates of the mesh nodes and will solve the system.
Thanks a lot, @dokken and @Ekrem_Ekici for your helpful and prompt suggestions. Also, just to get your comments on what I am trying to solve on a spherical surface, below I copy my code.
If you can kindly tell me whether this makes sense or not, that would be great. Many thanks again for your help!
### Linear elasticity example FEniCSx tutorial ###
######################################################################
### importing packages and defining parameters starts ###
######################################################################
# Scaled variable
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
import pyvista
from dolfinx import mesh, fem, plot, io
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import trimesh
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx,
exp, grad, inner, pi, sin, cos)
#L = 1
#W = 0.2
mu = 1 ### 0.001, 1
#rho = 1
#delta = W / L
#gamma = 0.4 * delta**2
#beta = 1.25
#lambda_ = beta
lambda_ = 1
print("Value of mu: {}, and lambda: {}".format(mu,lambda_))
#g = gamma
######################################################################
### importing packages and defining parameters ends ###
######################################################################
######################################################################
### defining mesh starts ###
######################################################################
# domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
# [20, 6, 6], cell_type=mesh.CellType.hexahedron)
####------------------------------####
gdim = 3
shape = "triangle"
degree = 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
meshUsingTrimesh = trimesh.creation.icosphere(subdivisions=6, radius=1.0)
x = np.array(meshUsingTrimesh.vertices.tolist())
cells = np.array(meshUsingTrimesh.faces.tolist())
mesh1 = mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)
####------------------------------####
#V = fem.VectorFunctionSpace(domain, ("Lagrange", 1))
V = fem.functionspace(mesh1, ("Lagrange", 1, (mesh1.geometry.dim, )))
# Listing space points ### *********************************************************
dof_coordinates = V.tabulate_dof_coordinates()
df = pd.DataFrame(dof_coordinates)
df.to_csv('sphereNodesStokes.csv', index = None, header=False)
######################################################################
### defining mesh ends ###
######################################################################
######################################################################
### boundary conditions starts ###
######################################################################
### clamped BC starts ###
# def clamped_boundary(x):
# return np.isclose(x[0], 0)
#
#
# fdim = domain.topology.dim - 1
# boundary_facets = mesh.locate_entities_boundary(domain, 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)
bc = []
### clamped BC ends ###
### traction free BC starts on other faces ###
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
### traction free BC ends on other faces ###
### defining surface integral over the boundary starts ###
ds = ufl.Measure("ds", domain=domain)
### defining surface integral over the boundary ends ###
######################################################################
### boundary conditions ends ###
######################################################################
######################################################################
### variational formulation starts ###
######################################################################
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = SpatialCoordinate(mesh1)
#f = fem.Constant(domain, default_scalar_type((0, 0, 0.001)))
#f = as_vector((1 * cos(2 * pi * x[0] + 2 * pi * x[1]), 1 * cos(2 * pi * x[0] + 2 * pi * x[1]), 1 * cos(2 * pi * x[0] + 2 * pi * x[1])))
#f = as_vector((x[0], x[1], x[2]))
f = as_vector((-0.1*x[1], 0.1*x[0], 0))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
######################################################################
### variational formulation ends ###
######################################################################
######################################################################
### solving variational problem starts ###
######################################################################
problem = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
######################################################################
### solving variational problem ends ###
######################################################################
usol = uh.x.array
numpoints = int(len(usol)/3)
usol = np.reshape(usol, (numpoints,3))
#print(usol.tolist())
df = pd.DataFrame(usol)
df.to_csv('sphereDispStokes.csv', index = None, header=False)
print("Number of nodes {}".format(len(dof_coordinates)))
print("Number of sol points {}".format(len(usol)))
# ######################################################################
# ### visualization starts ###
# ######################################################################
# # file = "stokesFlow/outdata/planenodes.csv".format()
# file = "./sphereNodesStokes.csv".format()
# df = pd.read_csv(file, delimiter=",", header=None, skiprows=0)
# nodes = df.to_numpy()
# ######
# # file = "stokesFlow/outdata/planedisplacements.csv".format()
# file = "./sphereDispStokes.csv".format()
# df = pd.read_csv(file, delimiter=",", header=None, skiprows=0)
# vals = df.to_numpy()
#
# #################################################
# # %matplotlib qt
# ### for interactive plots
#
# fig = plt.figure(figsize=(8, 8))
# ax = plt.axes(projection='3d')
# # Vector origin location
# gap = 50
# X = nodes[::gap, 0]
# Y = nodes[::gap, 1]
# Z = nodes[::gap, 2]
#
# # Directional vectors
# U = vals[::gap, 0]
# V = vals[::gap, 1]
# W = vals[::gap, 2]
#
# # Creating plot
# # plt.quiver(X, Y, U, V, color='b', units='xy', scale=1)
# ax.quiver(X, Y, Z, U, V, W, color='b')
#
# ax.set_title("Stokes Flow on Sphere Fenicsx")
# ax.set_xlabel('x')
# ax.set_ylabel('y')
# ax.set_zlabel('z')
#
# # Show plot with grid
# plt.grid()
# plt.show()
# ######################################################################
# ### visualization ends ###
# ######################################################################
# ######################################################################
# ### stress computation starts ###
# ######################################################################
# s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
# von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))
# V_von_mises = fem.FunctionSpace(domain, ("DG", 0))
# stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
# stresses = fem.Function(V_von_mises)
# stresses.interpolate(stress_expr)
# warped.cell_data["VonMises"] = stresses.vector.array
# warped.set_active_scalars("VonMises")
# p = pyvista.Plotter()
# p.add_mesh(warped)
# p.show_axes()
# if not pyvista.OFF_SCREEN:
# p.show()
# else:
# stress_figure = p.screenshot(f"stresses.png")
# ######################################################################
# ### stress computation ends ###
# ######################################################################
Plotting the solution looks like this (I plotted the solution vectors at the nodes of the mesh):
I do not know what triMesh is, looks like a mesh generator.
I would prefer to use gmsh to generate dolfinx meshes. You can check this demo to generate spherical mesh. And use gmshio.model_to_mesh for importing it into your main calculation script.
For plotting, save your solution to XDMF file and use ParaView visualize it. The plot you share is not clear to me.
Thanks, @Ekrem_Ekici, for your reply. Very much appreciated.
Yes, triMesh is a mesh generator indeed. Thanks for sharing the links for gmsh, I will have a look.
I have used matplotlib in a separate code for plotting the data (which I have not shared here), from the csv files generated by the FEniCSx code I shared. The plotting part is not executed in my shared code.
Also, just to explain what I solved for and plotted here, I chose a tangential force field for the function “f” in the code. Then obtained the solution for displacements (3D displacements but tangential to the surface of the sphere, I suppose) at each node using FEniCSx, which are saved in the file ‘sphereDispStokes.csv’ (x-, y-, and z-components are in the 1st, 2nd, and 3rd columns respectively). Similarly, the node-coordinates are saved in the file ‘sphereNodesStokes.csv’. Using these two files I plotted the vector field of displacements at the nodes on the sphere (using matplotlib in a separate code, just used for plotting).
Thanks!
