### Compute surface area over a given boundary of the mesh

73

views

0

Hi there ,

I am using the following code to compute volume of my mesh.

Thanks.

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

```
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)
```

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 :-

Best,

Ovais

`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.

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)