Error in defining the integral over subdomain

[Akshay]

In order to calculate J-Integral, I want to calculate integration over subdomain. So I defined the subdomain as

But error encountered as

Please help me to define subdomain for integration

Complete Code is as Follows:

from dolfin import *
import numpy as np
import math
parameters[‘form_compiler’][‘cpp_optimize’] = True
parameters[‘form_compiler’][‘optimize’] = True
from mshr import*
set_log_active(False)
import matplotlib.pyplot as plt
from fenics import *
‘’’--------------------------------------------------------------------

Mesh Generation

#--------------------------------------------------------------------’’’

#mesh_size = 400
#mesh = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5),mesh_size,mesh_size)
#mesh = UnitSquareMesh(mesh_size,mesh_size)
mesh = Mesh(“mesh7.xml”)
plot(mesh,title=“Mesh”)
plt.xlabel(“X”)
plt.ylabel(“Y”)
plt.show()

‘’‘Printing Mesh parameters’’’

print(‘calculating the Length Scale Parameter and Total no of Elements’)
print (‘Length Scale Parameter =’,2*mesh.hmin())
print ('Total no of Elements = ', mesh.num_cells())

ndim = mesh.geometry().dim()

‘’’--------------------------------------------------------------------

Defining Functionspaces

#--------------------------------------------------------------------’’’

V = FunctionSpace(mesh,“CG”,1)
W = VectorFunctionSpace(mesh,“CG”,1)
WW = FunctionSpace(mesh,“DG”,0)

p , q = TrialFunction(V) , TestFunction(V)
u , v = TrialFunction(W) , TestFunction(W)

‘’’--------------------------------------------------------------------

Material Properties

#--------------------------------------------------------------------’’’

E = 210e3
nu = 0.31
k1c = 9.6
G_c = 1000*(1-nu*nu)k1c**2/E
#l = 0.015 # Characterstic Length scale parameter
l = 2mesh.hmin()

lmbda , mu = Constant(Enu/((1.0 + nu )(1.0-2.0nu))) , Constant(E/(2(1+nu)))

‘’’--------------------------------------------------------------------

Constitutive Relations

#--------------------------------------------------------------------’’’

def epsilon(u):
return sym(grad(u))

def sigma(u):
return 2muepsilon(u) + lmbda*tr(epsilon(u))*Identity(ndim)

‘’‘Strain Energy Functional and History Field Variable’’’

def psi(u):
return 0.5*(lmbda + mu)(0.5(tr(epsilon(u)) + abs(tr(epsilon(u)))))**2 +
mu*inner(dev(epsilon(u)),dev(epsilon(u)))

def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)

def Crack(x):
return abs(x[1]) < (2*mesh.hmin()) and x[0] <= 0.0

‘’’--------------------------------------------------------------------

Subdividing Boundaries

#--------------------------------------------------------------------’’’
top = CompiledSubDomain(“near(x[1],0.5,1e-8) and on_boundary”)
bottom = CompiledSubDomain(“near(x[1],-0.5,1e-8) and on_boundary”)
left = CompiledSubDomain("(near(x[0],-0.5) or near(x[0],0.5)) and on_boundary")
domain1= CompiledSubDomain(“x[0]*x[0] + x[1]*x[1] - 0.0625 > 0 and x[0]*x[0] + x[1]*x[1] - 0.25 < 0”)

‘’’--------------------------------------------------------------------

Boundaries Conditions

#--------------------------------------------------------------------’’’

load = Expression(“t”,t=0,degree=1)
bcbot = DirichletBC(W,Constant((0,0)),bottom)
bcleft = DirichletBC(W.sub(0),Constant(0),left)
bctop = DirichletBC(W.sub(1),load,top)

bc_u = [bcbot,bctop]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]

‘’’--------------------------------------------------------------------

Numbering the Boundaries

#--------------------------------------------------------------------’’’

Boundaries = MeshFunction(‘size_t’,mesh,ndim-1)
Boundaries.set_all(0)
top.mark(Boundaries,1)
bottom.mark(Boundaries,2)
ds = Measure(‘ds’)(subdomain_data = Boundaries)
n = FacetNormal(mesh)

boundaries = MeshFunction(‘size_t’,mesh,mesh.topology().dim(),0)
boundaries.set_all(0)
domain1.mark(boundaries,1)
dA = Measure(‘dA’)(subdomain_data = boundaries)

‘’’--------------------------------------------------------------------

Variational Formulation

#--------------------------------------------------------------------’’’

unew , uold = Function(W) , Function(W)
pnew , pold ,Hold = Function(V) , Function(V) , Function(V)
E_du = ((1 - pold)**2)inner(sigma(u),grad(v))dx
E_phi = (G_clinner(grad§,grad(q)) + (((G_c/l) + 2*H(uold,unew,Hold))inner(p,q))
- 2H(uold,unew,Hold)*q)*dx

p_disp = LinearVariationalProblem(lhs(E_du),rhs(E_du),unew,bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi),rhs(E_phi),pnew,bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)

‘’’--------------------------------------------------------------------

Initialization of Iterative Process parameters

#--------------------------------------------------------------------’’’

‘’’--------------------------------------------------------------------

Staggered Scheme

#--------------------------------------------------------------------’’’

t = 0
u_r = 0.007
deltaT = 0.01
tol = 1e-3
conc_f = File ("./ResultsDir/phi.pvd")
fname = open(‘ForcevsDisp.txt’, ‘w’)

uyy_0 =[]
Reaction_0 = []

Staggered scheme

while t<=0.45:
t += deltaT
if t >=0.34:
deltaT = 0.001
load.t=t*u_r
iter = 0
err = 1

while err > tol:
    iter += 1
    solver_disp.solve()
    solver_phi.solve()
    err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
    err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
    err = max(err_u,err_phi)
    
    uold.assign(unew)
    pold.assign(pnew)
    Hold.assign(project(psi(unew), WW))

    if err < tol:
	
        print ('Iterations:', iter, ', Total time', t)

        if round(t*1e4) % 10 == 0:
            conc_f << pnew

            Traction = dot(sigma(unew),n)
            fy = Traction[1]*ds(1)
            fname.write(str(t*u_r) + "\t")
            fname.write(str(assemble(fy)) + "\n")
            uyy_0.append(t*u_r)
            Reaction_0.append(assemble(fy))
            print(assemble(fy))

fname.close()

print (‘Simulation done with no error’)

‘’’--------------------------------------------------------------------

Post Processing

#--------------------------------------------------------------------’’’

ux,uy = split(unew)

plot(ux)
plt.colorbar(plot(ux))
plt.xlabel(‘x’)
plt.ylabel(‘ux’)
plt.show()

plot(uy)
plt.colorbar(plot(uy))
plt.xlabel(‘x’)
plt.ylabel(‘uy’)
plt.show()

#load/displacement applied
plot(unew,title =‘U’ )
plt.colorbar(plot(unew,title = ‘U’))
plt.xlabel(‘x’)
plt.ylabel(‘y’)
plt.show()

stress(yy)

plot(sigma(unew)[1,1],title= "\sigma_{yy} [MPa]")
plt.colorbar(plot(sigma(unew)[1,1]))
plt.show()

#Phase-Field variable
plot(pnew,range_min = 0.,range_max = 1., title = ‘Phase Field Parameter’)
plt.colorbar(plot(pnew,range_min = 0.,range_max = 1))
plt.xlabel(‘x’)
plt.ylabel(‘Phi’)
plt.show()

#Reaction Vs Displacement
plt.plot(uyy_0,Reaction_0)
plt.xlabel(“Displacement”)
plt.ylabel(“Reaction Force”)

‘’’--------------------------------------------------------------------

Writing Excel File

#--------------------------------------------------------------------’’’

from xlsxwriter import *
outWorkbook = Workbook(“out.xlsx”)
outsheet = outWorkbook.add_worksheet()

for i in range(len(uyy_0)):
outsheet.write(i,0,uyy_0[i])

for j in range(len(Reaction_0)):
outsheet.write(j,1,Reaction_0[j])
print(“Results written successfully”)

[dokken]

Please encapsulate your code with ``` to for at it properly. However, i do not know Which measure A is supposed to be, as the measures in ufl are defined as:

_integral_types = [
    # === Integration over full topological dimension:
    ("cell", "dx"),                     # Over cells of a mesh
    # ("macro_cell", "dE"),              # Over a group of adjacent cells (TODO: Arbitrary cell group? Not currently used.)

    # === Integration over topological dimension - 1:
    ("exterior_facet", "ds"),           # Over one-sided exterior facets of a mesh
    ("interior_facet", "dS"),           # Over two-sided facets between pairs of adjacent cells of a mesh

    # === Integration over topological dimension 0
    ("vertex", "dP"),                    # Over vertices of a mesh
    # ("vertex", "dV"),                  # TODO: Use this over vertices?
    # ("point", "dP"),                   # TODO: Use this over arbitrary points inside cells?

    # === Integration over custom domains
    ("custom", "dc"),                 # Over custom user-defined domains (run-time quadrature points)
    ("cutcell", "dC"),                # Over a cell with some part cut away (run-time quadrature points)
    ("interface", "dI"),              # Over a facet fragment overlapping with two or more cells (run-time quadrature points)
    ("overlap", "dO"),                # Over a cell fragment overlapping with two or more cells (run-time quadrature points)

[Akshay]

I encounter above error. I modified my code according to infromation provided by you.

‘’’
from dolfin import *
import numpy as np
import math
parameters[‘form_compiler’][‘cpp_optimize’] = True
parameters[‘form_compiler’][‘optimize’] = True
from mshr import*
set_log_active(False)
import matplotlib.pyplot as plt
from fenics import *
‘’’--------------------------------------------------------------------

Mesh Generation

#--------------------------------------------------------------------’’’

mesh_size = 64
#mesh = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5),mesh_size,mesh_size)
mesh = UnitSquareMesh(mesh_size,mesh_size)
#mesh = Mesh(“mesh7.xml”)
plot(mesh,title=“Mesh”)
plt.xlabel(“X”)
plt.ylabel(“Y”)
plt.show()

‘’‘Printing Mesh parameters’’’

print(‘calculating the Length Scale Parameter and Total no of Elements’)
print (‘Length Scale Parameter =’,2*mesh.hmin())
print ('Total no of Elements = ', mesh.num_cells())

ndim = mesh.geometry().dim()

‘’’--------------------------------------------------------------------

Defining Functionspaces

#--------------------------------------------------------------------’’’

V = FunctionSpace(mesh,“CG”,1)
W = VectorFunctionSpace(mesh,“CG”,1)
WW = FunctionSpace(mesh,“DG”,0)

p , q = TrialFunction(V) , TestFunction(V)
u , v = TrialFunction(W) , TestFunction(W)

‘’’--------------------------------------------------------------------

Material Properties

#--------------------------------------------------------------------’’’

E = 210e3
nu = 0.31
k1c = 9.6
G_c = 1000*(1-nu*nu)k1c**2/E
#l = 0.015 # Characterstic Length scale parameter
l = 2mesh.hmin()

lmbda , mu = Constant(Enu/((1.0 + nu )(1.0-2.0nu))) , Constant(E/(2(1+nu)))

‘’’--------------------------------------------------------------------

Constitutive Relations

#--------------------------------------------------------------------’’’

def epsilon(u):
return sym(grad(u))

def sigma(u):
return 2muepsilon(u) + lmbda*tr(epsilon(u))*Identity(ndim)

‘’‘Strain Energy Functional and History Field Variable’’’

def psi(u):
return 0.5*(lmbda + mu)(0.5(tr(epsilon(u)) + abs(tr(epsilon(u)))))**2 +
mu*inner(dev(epsilon(u)),dev(epsilon(u)))

def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)

def Crack(x):
return abs(x[1]) < (2*mesh.hmin()) and x[0] <= 0.0

‘’’--------------------------------------------------------------------

Subdividing Boundaries

#--------------------------------------------------------------------’’’
top = CompiledSubDomain(“near(x[1],0.5,1e-8) and on_boundary”)
bottom = CompiledSubDomain(“near(x[1],-0.5,1e-8) and on_boundary”)
left = CompiledSubDomain("(near(x[0],-0.5) or near(x[0],0.5)) and on_boundary")
domain1= CompiledSubDomain(’(x[0]*x[0] + x[1]*x[1] - 0.0625) >= 0’)

‘’’--------------------------------------------------------------------

Boundaries Conditions

#--------------------------------------------------------------------’’’

load = Expression(“t”,t=0,degree=1)
bcbot = DirichletBC(W,Constant((0,0)),bottom)
bcleft = DirichletBC(W.sub(0),Constant(0),left)
bctop = DirichletBC(W.sub(1),load,top)

bc_u = [bcbot,bctop]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]

‘’’--------------------------------------------------------------------

Numbering the Boundaries

#--------------------------------------------------------------------’’’

Boundaries = MeshFunction(‘size_t’,mesh,ndim-1)
Boundaries.set_all(0)
top.mark(Boundaries,1)
bottom.mark(Boundaries,2)
ds = Measure(‘ds’)(subdomain_data = Boundaries)
n = FacetNormal(mesh)

boundaries = MeshFunction(‘size_t’,mesh,mesh.topology().dim(),0)
boundaries.set_all(0)
domain1.mark(boundaries,1)
dc = Measure(‘dc’)(subdomain_data = boundaries)

‘’’--------------------------------------------------------------------

Variational Formulation

#--------------------------------------------------------------------’’’

unew , uold = Function(W) , Function(W)
pnew , pold ,Hold = Function(V) , Function(V) , Function(V)
E_du = ((1 - pold)**2)inner(sigma(u),grad(v))dx
E_phi = (G_clinner(grad§,grad(q)) + (((G_c/l) + 2*H(uold,unew,Hold))inner(p,q))
- 2H(uold,unew,Hold)*q)*dx

p_disp = LinearVariationalProblem(lhs(E_du),rhs(E_du),unew,bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi),rhs(E_phi),pnew,bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)

‘’’--------------------------------------------------------------------

Initialization of Iterative Process parameters

#--------------------------------------------------------------------’’’

‘’’--------------------------------------------------------------------

Staggered Scheme

#--------------------------------------------------------------------’’’

t = 0
u_r = 0.007
deltaT = 0.01
tol = 1e-3
conc_f = File ("./ResultsDir/phi.pvd")
fname = open(‘ForcevsDisp.txt’, ‘w’)

uyy_0 =[]
Reaction_0 = []

Staggered scheme

while t<=0.45:
t += deltaT
if t >=0.34:
deltaT = 0.001
load.t=t*u_r
iter = 0
err = 1

while err > tol:
    iter += 1
    solver_disp.solve()
    solver_phi.solve()
    err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
    err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
    err = max(err_u,err_phi)
    
    uold.assign(unew)
    pold.assign(pnew)
    Hold.assign(project(psi(unew), WW))

    if err < tol:
	
        print ('Iterations:', iter, ', Total time', t)

        if round(t*1e4) % 10 == 0:
            conc_f << pnew

            Traction = dot(sigma(unew),n)
            fy = Traction[1]*ds(1)
            fname.write(str(t*u_r) + "\t")
            fname.write(str(assemble(fy)) + "\n")
            uyy_0.append(t*u_r)
            Reaction_0.append(assemble(fy))
            print(assemble(fy))

fname.close()

print (‘Simulation done with no error’)
‘’’

[dokken]

You are using the wrong backticks to encapsulate code, Please Edit your post with the correct ones.

[dokken]

@Akshay

boundaries.set_all(0)
domain1.mark(boundaries,1)
dc = Measure(‘dc’)(subdomain_data = boundaries

Also, Please adds description of what you want to obtain with your dc integral, and remove all code the is after the error, to make the example smaller.