### Assemble over boundary faces through Submesh

230
views
0
5 months ago by
Hi there,
I have read an earlier post on Assemble over boundary faces through BoundaryMesh (here)  where as I am trying to assemble an integral over 2 x boundary faces using Submesh.
I was computing volume of my mesh after solving a nonlinear elasticity problem for displacement u occurring due to an external force applied over boundaries through Nanson's formula as under:-

using following code:-

1 def Compute_Volume (u):
2    X = SpatialCoordinate ( mesh )
3    vol = assemble (( -1./3.) * dot(X + u, J* inv (F).T* n_mesh )*ds( endo ))
4    return vol *1 e6 # scale volume to ml​

Now i have a different mesh for which i have 2 x surfaces (defined through facet function)  and I want to compute the volume of both surface partitions of my mesh separately after solving a nonlinear elasticity problem for displacement u as under:-(posting almost complete code because there might an issue in some other part of code which may be causing problem in computation of volume).
mesh = Mesh('biv_mesh.xml')
.
.

# Set up mesh boundaries and normal
boundaries = Make_SubDomains(mesh)					           # Identify boundaries as facets and add markes to facets through call of Make_SubDomains function
ds = ds[boundaries]
n_mesh = FacetNormal(mesh)

# Space and functions
U = VectorFunctionSpace(mesh,"Lagrange", 2)
V = VectorFunctionSpace(mesh,"Lagrange", 1)
Q = FunctionSpace(mesh,"Lagrange",1)
W = MixedFunctionSpace([U,V,Q])

uvp = Function(W)
u, v, p = split(uvp)
du, dv, dp = TestFunctions(W)
.
.
.
f_0 = Function(FiberSpace) 	# fiber direction
s_0 = Function(FiberSpace) 	# sheet vector
n_0 = Function(FiberSpace) 	# sheet normal

f = F*f_0
s = F*s_0
n = F*n_0
.
.

Jac = derivative(eq, uvp)
# solve the problem for displacement , velocity and pressure
Solve_System(T_a,p0,p1)

# Finally compute volume
left_submesh = SubMesh(mesh,boundaries,30)
X = SpatialCoordinate(left_submesh)
left_vol = 1e6*abs(assemble((-1.0/3.0)*(dot(X + u, J*inv(F).T*n_mesh)*ds(30))))
Although i am able to solve the problem but when i try to compute volume i get following error:-
fenics@9ded2a1c8504:~/shared/ms22/bivent_soln/Experiment_Snippets\$ python Compute_volume_undeformed.py
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 6.969e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 9.675e-06 (tol = 1.000e-10) r (rel) = 1.388e-08 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 3.058e-12 (tol = 1.000e-10) r (rel) = 4.388e-15 (tol = 1.000e-09)
Newton solver finished in 2 iterations and 2 linear solver iterations.
An Integral without a Domain is now illegal.
Traceback (most recent call last):
File "Compute_volume_undeformed.py", line 184, in <module>
left_vol = 1e6*abs(assemble((-1.0/3.0)*(dot(X + u, J*inv(F).T*n_mesh)*ds(30))))
File "/home/fenics/build/lib/python2.7/site-packages/ufl/measure.py", line 419, in __rmul__
return Form([integral])
File "/home/fenics/build/lib/python2.7/site-packages/ufl/form.py", line 92, in __init__
self._integrals = _sorted_integrals(integrals)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/form.py", line 45, in _sorted_integrals
ufl_assert(d is not None, "An Integral without a Domain is now illegal.")
File "/home/fenics/build/lib/python2.7/site-packages/ufl/assertions.py", line 37, in ufl_assert
if not condition: error(*message)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: An Integral without a Domain is now illegal.​

Can anyone suggest a reason of the error ?.. The only difference between the previous way of computing volume and now is that now i am using SpatialCoordinate of a Submesh instead of complete mesh. Is it completely wrong ?.

Can we use assemble over boundary faces using Submesh instead of mesh ?? and if i want to use assemble over a boundary face , is there some other way ?

Some guidance / help is required, please.
Best regards.

Community: FEniCS Project
1
left_vol = 1e-06*abs(assemble(...*ds(30, domain=mesh)))​