dolfin team mailing list archive

[luckas <question151902@xxxxxxxxxxxxxxxxxxxxx>]

New question #151902 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/151902

Hi,
I have  a Problem  with the solving of  integral equation in Fenics  , so my question  is ,how to write the Syntax of  this equation : W = 0.5∫▒〖(■(σ11&σ12@σ21&σ22)〗).(■(ε11&ε12@ε21&ε22))dx ,
Where the first σ Matrix is the tensor Sigma (v) line 113 and  the second  ε Matrix is the tensor (grad(v)+grad(v).T)   line 114 in the code.
I will be happy  if anyone will solve it for me.

 
 the code is :



#Unterschiedlich Materialzuweisung
# Last changed: 2010-05-11

from dolfin import *

not_working_in_parallel("mesh smoothing demo")

# Create rectangular mesh
mesh = Rectangle(0.0, 0.0, 0.8, 0.8, 20, 20)
# Define a circular hole
center = Point(0.4, 0.4)
radius = 0.1

class Hole(SubDomain):

    def inside(self, x, on_boundary):
        r = sqrt((x[0] - center[0])**2 + (x[1] - center[0])**2)
        return r < 1.5*radius # slightly larger

    def snap(self, x):
        r = sqrt((x[0] - center[0])**2 + (x[1] - center[1])**2)
        if r < 1.5*radius:
            x[0] = center[0] + (radius / r)*(x[0] - center[0])
            x[1] = center[1] + (radius / r)*(x[1] - center[1])

# Refine and snap mesh
plot(mesh, title="Mesh 0")
num_refinements = 6
for i in range(num_refinements):

    # Mark cells for refinement
    markers = MeshFunction("bool", mesh, mesh.topology().dim())
    markers.set_all(False)
    for cell in cells(mesh):
        if cell.midpoint().distance(center) < 1.1*radius:
            markers[cell.index()] = True
        if cell.midpoint().distance(center) < 0.9*radius:
            markers[cell.index()] = False


    # Refine mesh
    mesh = refine(mesh, markers)



    # Save solution in VTK format
    file = File("andreas.pvd")
    file << mesh


    # Plot mesh
    plot(mesh, title=("Mesh %d" % (i + 1)))

# interactive()

class Material0(SubDomain):
   def inside(self, x, on_boundary):
       return True if x[1] <= 0.8 and x[0] <= 0.8 and sqrt((x[0] - center[0])**2 + (x[1] - center[1])**2) >= 0.1 else False

class Material1(SubDomain):
   def inside(self, x, on_boundary):
       return True if sqrt((x[0] - center[0])**2 + (x[1] - center[1])**2) < 0.1 else False


V = VectorFunctionSpace(mesh, 'CG' , 1)

left , right , top , bottom = compile_subdomains(["(fabs(x[0]) < DOLFIN_EPS) && on_boundary" ,
"(fabs(x[0] - 0.8) < DOLFIN_EPS) && on_boundary" ,
 "(fabs(x[1] - 0.8) < DOLFINE_EPS) && on_boundary" ,
"(fabs(x[1]) < DOLFIN_EPS) && on_boundary" ])

boundary_parts = MeshFunction("uint" ,mesh , mesh.topology().dim()-1)

Q = Function(V)
Q = Expression('-1.0*sin(pi*1.0*x[0]/0.4)*cos(pi*1.0*x[1]/0.4)')

u0 = Constant(0.0)
u1 = Constant(0.4)

bc = [DirichletBC(V.sub(0), u0, left), DirichletBC(V.sub(0), u1, right),
DirichletBC(V.sub(1), u0, bottom)]


v = TestFunction(V)
u = TrialFunction(V)

subdomains = MeshFunction('uint', mesh, 2)
# Mark subdomains with numbers 0 and 1
subdomain0 = Material0()
subdomain0.mark(subdomains, 0)
subdomain1 = Material1()
subdomain1.mark(subdomains, 1)

V0 = FunctionSpace(mesh, 'DG', 0)
G = Function(V0)
nu = Function(V0)
lmbd = Function(V0)

G_values = [48, 30] #GPa
for cell_no in range(len(subdomains.values())):
         subdomain_no = subdomains.values()[cell_no]
         G.vector()[cell_no] = G_values[subdomain_no]

nu_values = [0.34, 0.37]
for cell_no in range(len(subdomains.values())):
         subdomain_no = subdomains.values()[cell_no]
         nu.vector()[cell_no] = nu_values[subdomain_no]
         lmbd.vector()[cell_no] = 2.0*G_values[subdomain_no]*nu_values[subdomain_no]/(1.0-2.0*nu_values[subdomain_no])


#lmbd = 2.0*G*nu / (1.0 - 2.0*nu)

def sigma(v):
        return lmbd*tr(0.5*(grad(v)+grad(v).T))*Identity(v.cell().d) + 2*G*0.5*(grad(v)+grad(v).T)


n = V.cell().n
a = inner(grad(v), sigma(u))*dx
L = inner(v, Q*n)*ds(0)


A = assemble(a, mesh=mesh, exterior_facet_domains=boundary_parts)
b = assemble(L, mesh=mesh, exterior_facet_domains=boundary_parts)
for condition in bc: condition.apply(A, b)
u = Function(V)
solve(A, u.vector(), b)

file = File("Verschiebungen.pvd")
file << u

file = File("G.pvd")
file << G
file = File("nu.pvd")
file << nu
file = File("lmbd.pvd")
file << lmbd


-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.