Assemble over boundary faces through Submesh
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).
Although i am able to solve the problem but when i try to compute volume i get following error:-
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))))
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.