In [ ]:

    

import getfem as gf
import numpy as np
from scipy.sparse import linalg
from scipy import io


    

In [ ]:

    

h = 2

mo1 = gf.MesherObject('rectangle', [0., 0.], [100., 25.])
mo2 = gf.MesherObject('ball', [25., 12.5], 8.)
mo3 = gf.MesherObject('ball', [50., 12.5], 8.)
mo4 = gf.MesherObject('ball', [75., 12.5], 8.)
mo5 = gf.MesherObject('union', mo2, mo3, mo4)
mo  = gf.MesherObject('set minus', mo1, mo5)

mesh = gf.Mesh('generate', mo1, h, 2)
mesh = gf.Mesh('generate', mo, h, 2)


    

In [ ]:

    

# detect some boundary of the mesh
P = mesh.pts()
cright = (abs(P[0,:] - P[0,:].max()) < 1e-6)
cleft  = (abs(P[0,:] - P[0,:].min()) < 1e-6)
pidright = np.compress(cright,range(0,mesh.nbpts()))
pidleft  = np.compress(cleft,range(0,mesh.nbpts()))
fright = mesh.faces_from_pid(pidright)
fleft  = mesh.faces_from_pid(pidleft)
# create boundary region
NEUMANN_BOUNDARY = 1
DIRICHLET_BOUNDARY = 2
mesh.set_region(NEUMANN_BOUNDARY,fright)
mesh.set_region(DIRICHLET_BOUNDARY,fleft)


    

In [ ]:

    

mfp = gf.MeshFem(mesh, 1)
mfp.set_fem(gf.Fem('FEM_PK(2,1)'))


    

In [ ]:

    

x = mfp.basic_dof_nodes()[0,:]
y = mfp.basic_dof_nodes()[1,:]
tri = mfp.basic_dof_from_cvid()[0].reshape([-1,3])

import pylab as plt
%matplotlib inline

plt.figure(figsize=(10, 10))
plt.gca().set_aspect('equal')
plt.triplot(x,y,tri,color="red")


    

In [ ]:

    

mfd = gf.MeshFem(mesh, 1)
mfd.set_fem(gf.Fem('FEM_PK(2,1)'))


    

In [ ]:

    

mfu = gf.MeshFem(mesh, 2)
mfu.set_fem(gf.Fem('FEM_PK(2,2)'))


    

In [ ]:

    

mim = gf.MeshIm(mesh, gf.Integ('IM_TRIANGLE(3)'))


    

In [ ]:

    

Lambda = np.repeat([1.0],mfd.nbdof())
Mu = np.repeat([1.0],mfd.nbdof())
K = gf.asm_linear_elasticity(mim, mfu, mfd, Lambda, Mu)


    

In [ ]:

    

(H,R) = gf.asm_dirichlet(DIRICHLET_BOUNDARY, mim, mfu, mfd, mfd.eval('[[1, 0], [0, 1]]'), mfd.eval('[0,0]'))
(N,U0) = H.dirichlet_nullspace(R)


    

In [ ]:

    

K.save('mm', 'K.mtx'); K = io.mmread('K.mtx')
N.save('mm', 'N.mtx'); N = io.mmread('N.mtx')


    

In [ ]:

    

Nt = N.transpose()
KK = Nt*K*N


    

In [ ]:

    

F = gf.asm_boundary_source(NEUMANN_BOUNDARY, mim, mfu, mfd, np.repeat([[1.0],[0.0]], mfd.nbdof(), 1))


    

In [ ]:

    

FF = Nt*F


    

In [ ]:

    

UU = linalg.spsolve(KK, FF)


    

In [ ]:

    

U = N*UU
Ux = U.reshape(-1,2)[:,0]
Uy = U.reshape(-1,2)[:,1]


    

In [ ]:

    

x = mfp.basic_dof_nodes()[0,:]
y = mfp.basic_dof_nodes()[1,:]
tri = mfp.basic_dof_from_cvid()[0].reshape([-1,3])

import pylab as plt
%matplotlib inline

plt.figure(figsize=(10, 10))
plt.gca().set_aspect('equal')
plt.triplot(x,y,tri,color="red")

x = mfu.basic_dof_nodes()[0,:]
y = mfu.basic_dof_nodes()[1,:]
plt.quiver(x,y,np.repeat(Ux,2),np.repeat(Uy,2),scale=500,color="blue")


    

In [ ]: