Hi,
I am trying to solve a hyperelasticity problem on a 2D domain using the virtual displacements form
PK1:F *dx
# Modules Used
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
set_log_active(False)
# Number of Steps for Analytical Solution
n_steps = 21;
lam2 = np.linspace(1.0, 0.5, n_steps)
# Optimization options for the form compiler
parameters["form_compiler"]["representation"] = 'uflacs'
# Generate Mesh and function spaces
mesh = UnitSquareMesh(4,4);
x_corner = np.min(mesh.coordinates()[np.where(mesh.coordinates()[:,1]<1e-3),0])
class North(SubDomain):
def inside(self,x,on_boundary):
return near(x[1],1.0) and on_boundary
class South(SubDomain):
def inside(self,x,on_boundary):
return near(x[1],0.0) and on_boundary
class East(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],1.0) and on_boundary
class West(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],0.0) and on_boundary
def corner(x,on_boundary):
return near(x[0],x_corner) and near(x[1],0.0)
p_order = 2 # Polynomial Order with each element
V = VectorFunctionSpace(mesh,'CG', p_order)
T = TensorFunctionSpace(mesh, 'DG', p_order)
v = TestFunction(V) # Test Function
du = TrialFunction(V) # Incremental Displacement using 2D strain Energy
u= Function(V) # Displacement Control Solution
du2 = TrialFunction(V) # Incremental Displacement using 2D strain Energy
u2= Function(V) # Force Contorl Solution
# Initialize Domains
north=North()
south=South()
east =East()
west =West()
boundaries = FacetFunction("size_t",mesh)
boundaries.set_all(0);
# Label our boundaries
north.mark(boundaries,1);
east.mark(boundaries,2);
south.mark(boundaries,3);
west.mark(boundaries,4);
ds=Measure("ds", domain=mesh,subdomain_data = boundaries)
# Elasticity Parameters
E,nu =10.0, 0.499
C1,D1= Constant(E/(2*(1+nu))), Constant(E/(3*(1-2*nu)))
def strain_energy(u):
# Kinematics
d = len(u)
I = Identity(d)
F = I + grad(u)
C = F.T*F
Ic = tr(C)+1.0
J = det(F)
psi = 0.5*C1*(J**(-2.0/3)*Ic-3)+0.5*D1*(J-1)**2
return psi
def PK1(u):
d = len(u)
I = Identity(d)
F = variable(I + grad(u))
C = F.T*F
Ic=variable(tr(C))+1
J = variable(det(F))
C = variable(C)
psi = 0.5*C1*(J**(-2.0/3)*Ic-3)+0.5*D1*(J-1)**2
#s = diff(psi,F)
P = F*PK2(u)
return P
def PK2(u):
d = len(u)
I = Identity(d)
F = variable(I + grad(u))
C = F.T*F
E = variable(0.5*(C-I))
Ic=variable(tr(C))+1
J = variable(det(F))
C = variable(C)
psi = 0.5*C1*(J**(-2.0/3)*Ic-3)+0.5*D1*(J-1)**2
S = diff(psi,E)
return S
def GreenLagrange(u):
d = len(u)
I = Identity(d)
F = variable(I + grad(u))
C = F.T*F
E = variable(0.5*(C-I))
return E
# Initialize Arrays for plotting
R_d = np.zeros([n_steps,4]) # displacement control reaction forces in the 1,2,3,4,direction
l1_d = np.zeros(n_steps)
psi_d = np.zeros(n_steps)
du_h = 1-lam2 # increment
n =FacetNormal(mesh)
def extract_traction(u,i):
s = PK1(u);
return s[i,0]*n[0]+s[i,1]*n[1]
# Time Step
for i in range(n_steps):
# Dirichlet Boundary Conditions
bc1 = DirichletBC(V.sub(1), Constant(0.0),south) # Roller Bottom
bc2 = DirichletBC(V, Constant((0.0,0.0)),corner, method='pointwise') # pin at a single point
bc3 = DirichletBC(V.sub(1), Constant(-du_h[i]), north) # Displacement Top
# Displacement Boundary Conditions for Disp Control
bcs1= [bc1,bc2,bc3];
# Solve
Res = inner(PK1(u),grad(v))*dx
Jac = derivative(Res,u,du); # Jacobian
solve(Res == 0,u, bcs1, J = Jac)
plot(u, mode='displacement', interactive=True)
Error message:
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /build/dolfin-2bgr9b/dolfin-2017.1.0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
