dolfin team mailing list archive

[luckas <question151173@xxxxxxxxxxxxxxxxxxxxx>]

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

hi,

I have not made now, was so that I have two body differentiation also
the same Materiall have but now I want to wish that I in any
body have different material, for example in the district would
I have silver with nu = 0.37 and G = 30 GPa and in the other Rectangel
material for such cupfer with nu = 0.34 and G = 48 GPa.

I will send you the code and I just hope that you help me
can.



"This demo program demonstrates various algorithms for mesh smoothing."

__author__ = "Anders Logg (logg@xxxxxxxxx)"
__date__ = "2010-03-02"
__copyright__ = "Copyright (C) 2010 " + __author__
__license__  = "GNU LGPL Version 2.1"

# 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()

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)


nu = 0.3
G = 50 #GPa
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

interactive()


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