### Assemble over boundary faces through Submesh

230

views

0

Hi there,

I have read an earlier post on Assemble over boundary faces through

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

Now i have a different mesh for which i have 2 x surfaces (defined through facet function) and

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)
.
.
.
FiberSpace = VectorFunctionSpace(mesh,'Quadrature',4)
f_0 = Function(FiberSpace) # fiber direction
s_0 = Function(FiberSpace) # sheet vector
n_0 = Function(FiberSpace) # sheet normal
Load('fiber_fd.h5','f_0',f_0)
Load('fiber_fd.h5','s_0',s_0)
Load('fiber_fd.h5','n_0',n_0)
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

Please login to add an answer/comment or follow this question.

?