In the below MWE, I have a straight beam of length 20 units and height=width=0.3 units. The axis of the beam is oriented along the X-axis. I want to calculate the mean curvature of the top edge of the front surface of the beam (surface facing the screen). The front surface is marked with a physical tag 13 and the top line of the front surface has the physical tag 3. The Gmsh code is given after the MWE.
Although in the MWE, the mean curvature of the front top line will be 0 because it is a straight line, but in my actual computation the geometry is very complicated and also deforms with time.
To calculate the mean curvature, I tried to extract the front surface as a submesh and followed the suggestion in How to compute curvature of a boundary - #2 by dokken. Here is the MWE:
import numpy
import dolfin
from dolfin import *
import meshio
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
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)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
######################################################################################################################
msh = meshio.read("mwe.msh") # Length = 20 microns, H,W = 0.3 microns
######################################################################################################################
tetra_mesh = create_mesh(msh, "tetra", True)
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("surface.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("surface.xdmf") as infile:
infile.read(mvc, "name_to_read")
sf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
def assemble_edge(tag):
mv = MeshView.create(mf, tag)
return assemble(Constant(1.0)*dx(domain=mv)) # Line measure
dx_custom = Measure("ds", domain=mesh, subdomain_data=sf) # Area measure
dv_custom = Measure("dx", domain=mesh, subdomain_data=cf) # Volume measure
print(f"Volume = {assemble(Constant(1.0)*dv_custom(19))}")
print(f"frontbottomline = {assemble_edge(4)}")
print(f"fronttopline = {assemble_edge(3)}")
print(f"frontleftline = {assemble_edge(1)}")
print(f"frontrightline = {assemble_edge(2)}")
print(f"lefttopline = {assemble_edge(5)}")
print(f"leftbottomline = {assemble_edge(6)}")
print(f"righttopline = {assemble_edge(7)}")
print(f"rightbottomline = {assemble_edge(8)}")
print(f"backtopline = {assemble_edge(9)}")
print(f"backbottomline ={assemble_edge(10)}")
print(f"backleftline = {assemble_edge(11)}")
print(f"backrightline = {assemble_edge(12)}")
print(f"frontsurface = {assemble(Constant(1.0)*dx_custom(13))}")
print(f"leftsurface = {assemble(Constant(1.0)*dx_custom(14))}")
print(f"rightsurface = {assemble(Constant(1.0)*dx_custom(15))}")
print(f"topsurface = {assemble(Constant(1.0)*dx_custom(16))}")
print(f"bottomsurface = {assemble(Constant(1.0)*dx_custom(17))}")
print(f"backsurface = {assemble(Constant(1.0)*dx_custom(18))}")
##########################################
# The following code is taken from https://fenicsproject.discourse.group/t/how-to-compute-curvature-of-a-boundary/3193/2
def mesh_to_boundary(v, b_mesh):
"""
Returns a the boundary representation of the CG-1 function v
"""
# Extract the underlying volume and boundary meshes
mesh = v.function_space().mesh()
# We use a Dof->Vertex mapping to create a global
# array with all DOF values ordered by mesh vertices
DofToVert = dolfin.dof_to_vertex_map(v.function_space())
VGlobal = numpy.zeros(v.vector().size())
vec = v.vector().get_local()
for i in range(len(vec)):
Vert = dolfin.MeshEntity(mesh, 0, DofToVert[i])
globalIndex = Vert.global_index()
VGlobal[globalIndex] = vec[i]
VGlobal = SyncSum(VGlobal)
# Use the inverse mapping to se the DOF values of a boundary
# function
surface_space = dolfin.FunctionSpace(b_mesh, "CG", 1)
surface_function = dolfin.Function(surface_space)
mapa = b_mesh.entity_map(0)
DofToVert = dolfin.dof_to_vertex_map(dolfin.FunctionSpace(b_mesh, "CG", 1))
LocValues = surface_function.vector().get_local()
for i in range(len(LocValues)):
VolVert = dolfin.MeshEntity(mesh, 0, mapa[int(DofToVert[i])])
GlobalIndex = VolVert.global_index()
LocValues[i] = VGlobal[GlobalIndex]
surface_function.vector().set_local(LocValues)
surface_function.vector().apply('')
return surface_function
def vector_mesh_to_boundary(func, b_mesh):
v_split = func.split(deepcopy=True)
v_b = []
for v in v_split:
v_b.append(mesh_to_boundary(v, b_mesh))
Vb = dolfin.VectorFunctionSpace(b_mesh, "CG", 1)
vb_out = dolfin.Function(Vb)
scalar_to_vec = dolfin.FunctionAssigner(Vb, [v.function_space() for
v in v_b])
scalar_to_vec.assign(vb_out, v_b)
return vb_out
def SyncSum(vec):
""" Returns sum of vec over all mpi processes.
Each vec vector must have the same dimension for each MPI process """
comm = dolfin.MPI.comm_world
NormalsAllProcs = numpy.zeros(comm.Get_size() * len(vec), dtype=vec.dtype)
comm.Allgather(vec, NormalsAllProcs)
out = numpy.zeros(len(vec))
for j in range(comm.Get_size()):
out += NormalsAllProcs[len(vec) * j:len(vec) * (j + 1)]
return out
def boundary_to_mesh(f, mesh):
b_mesh = f.function_space().mesh()
SpaceV = dolfin.FunctionSpace(mesh, "CG", 1)
SpaceB = dolfin.FunctionSpace(b_mesh, "CG", 1)
F = dolfin.Function(SpaceV)
GValues = numpy.zeros(F.vector().size())
map = b_mesh.entity_map(0) # Vertex map from boundary mesh to parent mesh
d2v = dolfin.dof_to_vertex_map(SpaceB)
v2d = dolfin.vertex_to_dof_map(SpaceV)
dof = SpaceV.dofmap()
imin, imax = dof.ownership_range()
for i in range(f.vector().local_size()):
GVertID = dolfin.Vertex(b_mesh, d2v[i]).index() # Local Vertex ID for given dof on boundary mesh
PVertID = map[GVertID] # Local Vertex ID of parent mesh
PDof = v2d[PVertID] # Dof on parent mesh
value = f.vector()[i] # Value on local processor
GValues[dof.local_to_global_index(PDof)] = value
GValues = SyncSum(GValues)
F.vector().set_local(GValues[imin:imax])
F.vector().apply("")
return F
##################################################################################
fs=SubMesh(mesh,sf,13) # extracting the front surface
boundary_marker_fs = MeshFunction("size_t", fs, fs.topology().dim()-1, 0)
ncells2 = MeshFunction("size_t", fs, fs.topology().dim())
surface_marker_fs = MeshFunction("size_t", fs, fs.topology().dim(), 0)
vmap2 = fs.data().array("parent_vertex_indices", 0)
cmap2 = fs.data().array("parent_cell_indices", fs.topology().dim())
n = 0
for c in cells(fs):
parent_cell = Cell(mesh, cmap2[c.index()])
surface_marker_fs.array()[c.index()] = sf.array()[parent_cell.index()]
for f in facets(parent_cell):
for g in facets(c):
g_vertices = vmap2[g.entities(0)]
if set(f.entities(0)) == set(g_vertices):
boundary_marker_fs.array()[g.index()] = mf.array()[f.index()]
n=n+1
ds_fs = Measure('ds', domain=fs, subdomain_data=boundary_marker_fs) # Boundary marker for submesh fs
dx_fs = Measure('dx', domain=fs, subdomain_data=surface_marker_fs) # Surface marker for submesh fs
##################################################################################
# The following code is taken from https://fenicsproject.discourse.group/t/how-to-compute-curvature-of-a-boundary/3193/2
n = FacetNormal(fs)
V = VectorFunctionSpace(fs, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*ds_fs(3)
l = inner(n, v)*ds_fs(3) # Line measure
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
File("nh.pvd") << nh
bmesh = BoundaryMesh(fs, "exterior")
nb = vector_mesh_to_boundary(nh, bmesh)
Q = FunctionSpace(bmesh, "CG", 1)
p, q = TrialFunction(Q), TestFunction(Q)
a = inner(p,q)*dx_fs(13)
l = inner(div(nb), q)*dx_fs(13)
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
kappab = Function(Q)
solve(A, kappab.vector(), L)
kappa = boundary_to_mesh(kappab, fs)
print(assemble(kappa*ds_fs(3)))
File("kappa.pvd") << kappa
This gives me the error:
boundary_marker_fs.array()[g.index()] = mf.array()[f.index()]
IndexError: index 128652 is out of bounds for axis 0 with size 87177
Here is the Gmsh code for the geometry:
SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.0, 0.0, 1.0};
//+
Point(2) = {0.0, 0.0, 0.3, 1.0};
//+
Point(3) = {20.0, 0.0, 0.3, 1.0};
//+
Point(4) = {20.0, 0.0, 0.0, 1.0};
//+
Point(5) = {0.0, 0.3, 0.0, 1.0};
//+
Point(6) = {20.0, 0.3, 0.0, 1.0};
//+
Point(7) = {20.0, 0.3, 0.3, 1.0};
//+
Point(8) = {0.0, 0.3, 0.3, 1.0};
//+
Line(1) = {8, 7};
//+
Line(2) = {7, 3};
//+
Line(3) = {3, 2};
//+
Line(4) = {2, 8};
//+
Line(5) = {5, 8};
//+
Line(6) = {1, 2};
//+
Line(7) = {1, 5};
//+
Line(8) = {5, 6};
//+
Line(9) = {6, 7};
//+
Line(10) = {1, 4};
//+
Line(11) = {4, 3};
//+
Line(12) = {4, 6};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {6, 4, -5, -7};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {11, -2, -9, -12};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {8, -12, -10, 7};
//+
Plane Surface(4) = {4};
//+
Curve Loop(5) = {11, 3, -6, 10};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {1, -9, -8, 5};
//+
Plane Surface(6) = {6};
//+
Surface Loop(1) = {6, 1, 3, 5, 2, 4};
//+
Volume(1) = {1};
//+
Physical Curve("fronttop", 3) = {1};
//+
Physical Curve("frontbottom", 4) = {3};
//+
Physical Curve("frontleft", 1) = {4};
//+
Physical Curve("frontright", 2) = {2};
//+
Physical Curve("lefttop", 5) = {5};
//+
Physical Curve("leftbottom", 6) = {6};
//+
Physical Curve("righttop", 7) = {9};
//+
Physical Curve("rightbottom", 8) = {11};
//+
Physical Curve("backtop", 9) = {8};
//+
Physical Curve("backbottom", 10) = {10};
//+
Physical Curve("backleft", 11) = {7};
//+
Physical Curve("backright", 12) = {12};
//+
Physical Surface("front", 13) = {1};
//+
Physical Surface("left", 14) = {2};
//+
Physical Surface("right", 15) = {3};
//+
Physical Surface("top", 16) = {6};
//+
Physical Surface("bottom", 17) = {5};
//+
Physical Surface("back", 18) = {4};
//+
Physical Volume("vol", 19) = {1};
I am not sure what’s causing the error and how to fix it and calculate the mean curvature of a line in the geometry.
Any suggestion would be extremely helpful. Thanks in advance!