Error: Nonmatching dimensions for free indices with same id!

Hello Everyone,

I want to study piezoelectric effect in a 2D beam. The boundary conditions are: 1) the left most bottom corner of the beam is fixed in x and y direction and the rest of the nodes on the bottom edge are restricted only in y-direction 2) load is applied to its top edge 3) potential of 0.0 is given for the bottom edge 4) potential of 0.2 applied to its top edge

The variational formulation is as follows:

where:
fB = body forces
fN = traction forces
qB = volume electric charge given by Div(eletric_displacement)
qN = surface electric charge

My implementation is as follows:

from __future__ import print_function
from dolfin import *

mesh = RectangleMesh(Point(0.0, 0.0), Point(50.0,5.0), 50,10,"crossed")


piezo_Mat=as_matrix([[0.,0.,0.,0.,15,0.], 
                     [0.,0.,0.,15,0.,0.],
                     [30,30,35,0.,0.,0.]]) 

permit_Mat=as_matrix([[0.5,0.,0.], 
                      [0.,0.5,0.],
                      [0.,0.,0.6]]) 

V_mech=VectorElement("CG",mesh.ufl_cell(),2)   
V_elect=FiniteElement("CG",mesh.ufl_cell(),2) 
Mix_ele=MixedElement([V_mech,V_elect])        
Mixed_FS=FunctionSpace(mesh,Mix_ele) 

(u,phi)=TrialFunctions(Mixed_FS)
(v,si) = TestFunctions(Mixed_FS) 

V_ele=FunctionSpace(mesh,"DG",0)

tol=1E-14  
def bottom_fixed_boundary(x,on_boundary):
    return near(x[0],0,tol) and near(x[1],0,tol)
Fixed_left=Constant((0.,0.))    
bc_Fixed=DirichletBC(Mixed_FS.sub(0),Fixed_left,bottom_fixed_boundary,method='pointwise')

def bottom_right_boundary(x,on_boundary):
    return near(x[1],0.0,tol) and x[0]>0.0
Roller_right=Constant(0)
bc_roller=DirichletBC(Mixed_FS.sub(0).sub(1),Roller_right,bottom_right_boundary) 

ground_phi=Constant(0)
bc_bottom_ground=DirichletBC(Mixed_FS.sub(1),ground_phi,bottom_fixed_boundary)

bcs_all=[bc_Fixed,bc_roller,bc_bottom_ground]

nu=0.3  
E0=3000
alpha_val=100
rho0=10
rho_val=[] 

for cell_s in cells(mesh): 
    rho_val.append(rho0) 

class Top(SubDomain):
    tol=1E-14
    def inside(self,x,on_boundary):
        return near(x[1],5.0,tol)
top=Top() 

boundries=MeshFunction('size_t',mesh,1) 
boundries.set_all(0) 
top.mark(boundries,1) 

ds=Measure('ds',domain=mesh,subdomain_data=boundries)

def calculate_E(rho):
    E_updated=Function(V_ele) 
    E_array=E_updated.vector().get_local()
    for i, rho in enumerate (rho): 
        E_array[i]=E0
    E_updated.vector().set_local(E_array)
    return E_updated
E0=calculate_E(rho_val) 

def calculate_alpha(rho_values):
    alpha=Function(V_ele) 
    alpha_array=alpha.vector().get_local()
    for i, rho in enumerate (rho_values): 
        alpha_array[i]=alpha_val 
    alpha.vector().set_local(alpha_array)
    return alpha
alpha0=calculate_alpha(rho_val)     

def calculate_Lame_coefficients(E_val):
    mu_val=(E_val)/(2*(1+nu))  
    lmbda_val=(E_val*nu)/((1+nu)*(1-2*(nu))) 
    return mu_val,lmbda_val
mu, lmbda=calculate_Lame_coefficients(E0) 

def epsilon(u):
    strain_u=0.5*(grad(u)+grad(u).T)
    return strain_u

def strain_voigt(u):
    strain_3_1=as_vector([u[i].dx(i) for i in range(2)] +[u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])
    strain_6_1=as_vector([strain_3_1[0], strain_3_1[1], 0, 0, 0, strain_3_1[2]])
    return strain_6_1

def sigma(u,mu,lmbda):
    eps=epsilon(u)
    stress_u=lmbda*div(u)*Identity(d)+2*mu*eps 
    return stress_u 

def elect_field_E(phi):
    E_field_val=-grad(phi)
    E_field_val_3_1=as_vector([E_field_val[0], E_field_val[1], 0])
    return E_field_val_3_1

def elect_disp_D(u,phi):
    elect_disp=alpha0*((piezo_Mat*strain_voigt(u)) + (permit_Mat*elect_field_E(phi)))      
    return elect_disp

fB=Constant((0.0,0.0))  
fN=Constant(10.0)

qB=elect_disp_D(v,si)
qN=Constant(0.2)

d=u.geometric_dimension()

a = 2*mu*inner(epsilon(u),epsilon(v))*dx + lmbda*dot(tr(epsilon(u)),tr(epsilon(v)))*dx + alpha0*div((permit_Mat*elect_field_E(phi))-(piezo_Mat*strain_voigt(u)))*div(si)*dx   

L=dot(fB,v)*dx + v[1]*fN*ds(1) + div(qB)*si*dx + qN*si*ds(1) 

w=Function(Mixed_FS)   
(u,phi)= split(w) 

solve(a==L,w,bcs_all)

(u,phi)=w.split(True)

print("Displacement: ",u.vector().get_local())
print("Potential: ",phi.vector().get_local())

I am getting an error that says “Nonmatching dimensions for free indices with same ids!

I am not sure if my implementation of the variation formulations is right and what is causing the error. Can I please get some hint or guidance?

Thank you very much.

Please post code that reproduces the error. The above snippet is not sufficient.

Concerning div(qB)*si*dx the div, given that this is a mesh in 2d, expects vectors of length 2 but qB is a 3-vector. You need to resolve this inconsistency.

Thank you very much Miroslav for your kind reply. I tried to work on your suggestions but have one more question. The term qB represents volume electric charge and is given as:

qB = Div(D)

where: D = α ξ ε(u) + α β E(φ)

In this equation the piezoelectric matrix (ξ ) is of size 3X6 and the permittivity matrix (β) is of size 3X3.

For the matrix multiplication to be possible, in the code I have updated the strain tensor (ε(u)) to be of size 6X1 and electric field (E(φ)) i.e. grad(φ) to be a 3X1 matrix. Thus the final result of this term is of size 3X1.

I will be very grateful if you could please provide some more details and help me resolve the issue. Once again, thank you very much for your precious time and valuable guidance. I look forward to hearing from you soon.

If you are sure that for the 2d problem qB is correctly defined as a 3-vector then you need to change what div does. For example, with the following modification your code compiles

L_elect= sum(qB[i].dx(i) for i in range(2))*si*dx + (qN*si)*ds(1)

Note that this is just an illustration to point out the issue with the divergence; the last component of the vector is neglected and the solver produces nonsense.

Thank you very much Miroslav for your valuable guidance. I am working on it.