Integrate along the mark element inside the domain


166
views
2
5 months ago by
Hi, I am trying to integrate function  $\int_{\Gamma}W_sdy-T_i\frac{du_i}{dx}d\Gamma$ΓWsdyTiduidx dΓ    along the marked element.
where:  $W_s=\frac{1}{2}\left[\sigma_{xx}\frac{\partial u_x}{\partial x}+\sigma_{xy}\left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}\right)\frac{\partial u_x}{\partial y}+\sigma_{yy}\frac{\partial u_y}{\partial y}\right]$Ws=12 [σxxuxx +σxy(uxy +uyx )uxy +σyyuyy ] ,  
$T_i\cdot\frac{\partial u_i}{\partial y}=\left[\sigma_{xx}n_1+\sigma_{xy}n_2\right]\frac{\partial u_x}{\partial x}+\left[\sigma_{xy}n_1+\sigma_{yy}n_2\right]\frac{\partial u_y}{\partial y}$Ti·uiy =[σxxn1+σxyn2]uxx +[σxyn1+σyyn2]uyy ,
and  $d\Gamma=\sqrt{\left(\frac{\partial u_x}{\partial x}\right)^2+\left(\frac{\partial u_y}{\partial y}\right)^2}d\eta$dΓ=(uxx )2+(uyy )2dη
I am getting error while wrting the equations. UFLException: Can't add expressions with different free indices.

How can I fix this problem and integrate along the marked boundary
from fenics import *
mesh = UnitSquareMesh(30,30)
plot(mesh, interactive = True)
V = VectorFunctionSpace(mesh, 'CG', 1)
cells = CellFunction('size_t', mesh)
facets = FacetFunction('size_t', mesh)

class ele_circle(SubDomain):
    def inside(self, x, on_boundary):
        tol = 2e-2
        return near((x[0])**2 + (x[1])**2, 0.2**2, tol)

cells.set_all(0)

subdomain0 = ele_circle()
subdomain0.mark(cells,1)

plot(cells, interactive = True)

dA = Measure('ds', domain = mesh, subdomain_data = facets)
dV = Measure('dx', domain = mesh, subdomain_data = cells)

left = CompiledSubDomain('near(x[1],0) && on_boundary')
right = CompiledSubDomain('near(x[1],length) && on_boundary', length = 1.0)

facets.set_all(0)
right.mark(facets,1)

tr = Expression(('A+B*x[1]','0.0'), A = 2.0, B = 0.5)

null = Constant((0.0,0.0))
bc = [DirichletBC(V, null, left)]

u = TrialFunction(V)
del_u = TestFunction(V)

nu = 0.3
E = 210e3
G = E/(2.0*(1+nu))
la = 2.0*G*nu/(1.0-2.0*nu)
mu = G

delta = Identity(2)
i,j = indices(2)

 
eps = as_tensor(1.0/2.0*(u[i].dx(j) + u[j].dx(i)), (i,j))
sigma = as_tensor(la*eps[i,i]*delta[i,j] + 2.0*mu*eps[i,j], (i,j))


a = sigma[j,i]*del_u[i].dx(j)*dV
L = tr[i]*del_u[i]*dA(1)

disp = Function(V)

solve(a==L, disp, bcs = bc)

n = FacetNormal(mesh)

Ws = 0.5*(sigma[0,0]*u[0].dx(0) + sigma[0,1]*( u[0].dx(1) + u[1].dx(0))*u[0].dx(0) + sigma[1,1]*u[1].dx(1) )
Ts = ( sigma[0,0]*n[0] + sigma[0,1]*n[1] )*u[0].dx(0) + ( sigma[0,1]*n[0] + sigma[1,1]*n[1] )*u[1].dx(1)
?
Thanks for the time.

Community: FEniCS Project

1 Answer


1
5 months ago by
Your implementation of Ws and Ts does not in any way represent the formulas you have written down above. The error you get thrown is due to your indices not being conform with the Einstein index convention which is followed in UFL. Especially the part sigma[i,j]*n[i] + sigma[j,j]*n[j] )*u[j].dx(j) in the definition of Ts contains 5 coupled indices. This does not make sense and throws the error you quoted.

Keep in mind that the expression sigma[i,i] is internally expanded to sigma[0,0]+sigma[1,1]+sigma[2,2].

You can manually code the expanded formulas you gave by using sigma[0,0] for  $\sigma_{xx}$σxx  , u[0].dx(1) for  $\frac{\partial u_x}{\partial y}$uxy and so forth.
Thanks for the help. I have manually coded expression for Ws and Ts, Now It is not showing any error. I have updated my code. Now how Can I integrate along the marked element?
written 5 months ago by hirshikesh  
I am trying to integrate along the marked cells but getting error: ArityMismatch: Multiplying expressions with overlapping form argument number 1, argument is v_1.
How can I fix this Error?

Thanks for the time.
from fenics import *
mesh = UnitSquareMesh(30,30)
plot(mesh, interactive = True)
V = VectorFunctionSpace(mesh, 'CG', 1)
cells = CellFunction('size_t', mesh)
facets = FacetFunction('size_t', mesh)

class ele_circle(SubDomain):
    def inside(self, x, on_boundary):
        tol = 2e-2
        return near((x[0])**2 + (x[1])**2, 0.2**2, tol)

cells.set_all(0)

subdomain0 = ele_circle()
subdomain0.mark(cells,1)

plot(cells, interactive = True)

dA = Measure('ds', domain = mesh, subdomain_data = facets)
dV = Measure('dx', domain = mesh, subdomain_data = cells)

left = CompiledSubDomain('near(x[1],0) && on_boundary')
right = CompiledSubDomain('near(x[1],length) && on_boundary', length = 1.0)

facets.set_all(0)
right.mark(facets,1)

tr = Expression(('A+B*x[1]','0.0'), A = 2.0, B = 0.5)

null = Constant((0.0,0.0))
bc = [DirichletBC(V, null, left)]

u = TrialFunction(V)
del_u = TestFunction(V)

nu = 0.3
E = 210e3
G = E/(2.0*(1+nu))
la = 2.0*G*nu/(1.0-2.0*nu)
mu = G

delta = Identity(2)
i,j = indices(2)

 
eps = as_tensor(1.0/2.0*(u[i].dx(j) + u[j].dx(i)), (i,j))
sigma = as_tensor(la*eps[i,i]*delta[i,j] + 2.0*mu*eps[i,j], (i,j))


a = sigma[j,i]*del_u[i].dx(j)*dV
L = tr[i]*del_u[i]*dA(1)

disp = Function(V)

solve(a==L, disp, bcs = bc)


n = FacetNormal(mesh)
Ws = 0.5*(sigma[0,0]*u[0].dx(0) + sigma[0,1]*( u[0].dx(1) + u[1].dx(0))*u[0].dx(0) + sigma[1,1]*u[1].dx(1) )
Ts = ( sigma[0,0]*n[0] + sigma[0,1]*n[1] )*u[0].dx(0) + ( sigma[0,1]*n[0] + sigma[1,1]*n[1] )*u[1].dx(1)

Y = assemble((Ws*dA(1) ))
​
written 4 months ago by hirshikesh  
This is because you don't integrate your calculated displacement field disp but the TrialFunction u.
And are you sure about your expression to calculate the energy? Where does the factor $\frac{\partial u_x}{\partial y}$uxy after the brackets in the second term come from?
written 4 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »