Compute surface area over a given boundary of the mesh


73
views
0
6 weeks ago by
Ovais  
Hi there ,
I am using the following code to compute volume of my mesh. 
from dolfin import*

mesh = Mesh('boxx.xml')

def Make_SubDomains(mesh):
	mf = FacetFunction("size_t", mesh)						
	mf.set_all(0)												
	
	return mf
boundaries = Make_SubDomains(mesh)					           
ds = ds[boundaries] 										  
n_mesh = FacetNormal(mesh)	
X = SpatialCoordinate(mesh)
vol = abs(assemble((-1.0/3.0)*dot(X , n_mesh)*ds(0)))
print(vol)​
My mesh is of 20x7x3 dimensions. I am getting the correct volume through this code. I was wondering how can i compute surface area using a similar code. My end goal is to compute surface area of a prolate ellipsoid. Hence i cant use normal geometry here. This is a toy problem that i have created to communicate the same problem in an easy way.
Thanks.
Community: FEniCS Project

1 Answer


1
6 weeks ago by
Nate  
from dolfin import *
mesh = UnitCubeMesh(1, 1, 1)
facet_function = MeshFunction("size_t", mesh, 2, 0)
sub_boundary = CompiledSubDomain("near(x[0], 0.0) or near(x[1], 0.0) or near(x[2], 0.0)")
sub_boundary.mark(facet_function, 1)
ds = Measure("ds", subdomain_data=facet_function)
surfacearea = assemble(Constant(1.0)*ds(1, domain=mesh))
print(surfacearea)
thanks man! ..grateful.

I have tried this code and it works. I am using it in a nonlinear elasticity problem where the surface area changes with deformation. Any idea what should I add in the above code so that as I add u (deformation) to X (spatial coordinates) , the surface area should also change (in this case increase)
written 6 weeks ago by Ovais  
1
You can do this by explicitly integrating a form derived using Nanson's formula, or by using ALE.move() and integrating 1.0:

from dolfin import *
mesh = UnitCubeMesh(1, 1, 1)
facet_function = MeshFunction("size_t", mesh, 2, 0)
sub_boundary = CompiledSubDomain("near(x[0], 0.0) or near(x[1], 0.0) or near(x[2], 0.0)")
sub_boundary.mark(facet_function, 1)
ds = Measure("ds", subdomain_data=facet_function)
surfacearea = assemble(Constant(1.0)*ds(1, domain=mesh))
print(surfacearea)

# Rotate mesh (which should preserve area)
theta = pi/3.0
X = SpatialCoordinate(mesh)
disp = as_matrix([[cos(theta),-sin(theta),0.0],
                  [sin(theta),cos(theta),0.0],
                  [0.0,0.0,1.0]])*X - X

x = X + disp

def getMetric(phi):
    Dphi = grad(phi)
    return Dphi.T*Dphi
def surfaceJacobian(g,N):
    return sqrt(det(g)*inner(N,inv(g)*N))

g = getMetric(x)
N = FacetNormal(mesh)
da = surfaceJacobian(g,N)
print(assemble(da*ds(1, domain=mesh)))

# Scale rotated mesh by factor of 2.0 (which should multiply area by 4.0)
x = 2.0*(X+disp)

g = getMetric(x)
da = surfaceJacobian(g,N)
print(assemble(da*ds(1, domain=mesh)))

# Alternatively, if the displacement is a Function, you can use ALE.move()
# to hide Nanson's formula:
V = VectorFunctionSpace(mesh,"Lagrange",1)
disp = project(x - X,V)
ALE.move(mesh,disp)
print(assemble(Constant(1.0)*ds(1, domain=mesh)))

# Un-do mesh motion, to revert mesh to reference configuration:
disp.assign(-disp)
ALE.move(mesh,disp)
​
written 6 weeks ago by David Kamensky  
Thanks david , I found this very helpful. In this while I also figured out the following way :-

surface_area = assemble((1.0)*(J)*pow(dot(inv(C)*n_mesh,n_mesh),0.5)*ds(2))​
where J = det(F) , C is the Cauchy green deformation tensor , n_mesh is facet normal  and ds(2) have been marked as the entire outer boundary. Thanks for your interest in my question.
Best,
Ovais
written 6 weeks ago by Ovais  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »