Assemble over boundary faces through Submesh


230
views
0
5 months ago by
Ovais  
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)
.
.
.
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
1
what about using

left_vol = 1e-06*abs(assemble(...*ds(30, domain=mesh)))​
?
written 5 months ago by Hernán Mella  
let me try this one. Thanks.
written 5 months ago by Ovais  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »