I have modified the code in bits so that it is clear to discuss where the problem is arising -
Making geometry
import gmsh
import sys
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True,
"eliminate_zeros": True,
"precompute_basis_const": True,
"precompute_ip_const": True}
un = "trial"
b = 1
lc = 1e-1 * 10
scale = 1
def createGeometryAndMesh(b):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
P1 = gmsh.model.occ.addPoint(-b , 0, 0, lc)
P2 = gmsh.model.occ.addPoint(b, scale *0, 0, lc)
P3 = gmsh.model.occ.addPoint(-b, b , 0, lc)
P4 = gmsh.model.occ.addPoint(b , b , 0, lc)
L1 = gmsh.model.occ.addLine(P1, P2)
L2 = gmsh.model.occ.addLine(P2, P4)
L3 = gmsh.model.occ.addLine(P4, P3)
L4 = gmsh.model.occ.addLine(P3, P1)
Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1)
# Synchronize the CAD model
gmsh.model.occ.synchronize()
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [Surface1])
gmsh.model.setPhysicalName(2, pg1, "beam")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
def checkForEvent():
action = gmsh.onelab.getString("ONELAB/Action")
if len(action) and action[0] == "check":
gmsh.onelab.setString("ONELAB/Action", [""])
createGeometryAndMesh()
gmsh.graphics.draw()
return True
# Refine the mesh
gmsh.model.mesh.refine()
gmsh.write("geo.msh");
gmsh.finalize()
createGeometryAndMesh(b)
Importing mesh using meshio
import meshio
msh = meshio.read("geo.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data
triangle_mesh = meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read": [triangle_data]})
meshio.write("mesh.xdmf",triangle_mesh)
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
Defining boundary conditions
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
b_c = Boundary()
b_c.mark(boundary_markers, 3)
File("msh3.pvd") << boundary_markers
# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
V = VectorFunctionSpace(mesh, "CG", 2)
W_tensor = TensorFunctionSpace(mesh, "DG", 0)
def boundary_bot(x, on_boundary):
tol = 1E-7
return on_boundary and near(x[1], 0, tol)
bc_bot = DirichletBC(V, [0, 0, 0], boundary_bot)
def boundary_top(x, on_boundary):
tol = 1E-7
return on_boundary and near(x[1], b, tol)
disp = Constant([0.0, 0.0, 0.0])
bc_top = DirichletBC(V, disp, boundary_top)
bcs = [bc_bot, bc_top]
Defining function spaces and kinematics
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
u0 = Function(V) # Initial displacement
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)
dss = ds(subdomain_data=boundary_subdomains)
def kinematics(u):
# Kinematics
d = u.geometric_dimension() # Space dimension
I = Identity(d) # Identity tensor
F = variable(I + grad(u)) # Deformation gradient
C = F.T * F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
return Ic, J, F, C
Strain and energy values preprocessing
from scipy import interpolate
# Strain data points and energy values
strain_values = [0.0, 0.6, 0.19, 0.25]
energy_values = [0.0, 0.03, 0.005, 0.047]
if any(strain_values[i] >= strain_values[i+1] for i in range(len(strain_values)-1)):
# If there are repeated or non-unique values, preprocess the data
unique_strain_values, unique_indices = np.unique(strain_values, return_index=True)
unique_energy_values = [energy_values[i] for i in unique_indices]
strain_values = list(unique_strain_values)
energy_values = list(unique_energy_values)
# Calculate the spline for strain-energy interpolation
spline = CubicSpline(strain_values, energy_values)
# Solve variational problem
store_u = np.linspace(0, -0.4, 50 , endpoint=True)
normal = Constant((0, 1, 0))
# Define the stress tensor function space
V_tensor = TensorFunctionSpace(mesh, "DG", 0)
Iterative loop for setting up solver and solving
force = []
for i in store_u:
disp.assign(Constant([0.0, i, 0.0]))
Ic, J, F, C = kinematics(u)
# Define the custom strain energy density function
def psi(F, u, i):
strain = i
energy = Constant(spline(strain))
stress = diff(energy, F)
return energy, stress
psi, stress = psi(F, u, i)
# Total potential energy
Pi= psi * dx
# Compute first variation of Pi (directional derivative about u in the direction of v)
dPi = derivative(Pi, u, v)
# Compute second variation of Pi (directional derivative about du in the direction of v)
ddPi= derivative(dPi, u, du)
# Define the nonlinear variational problem with Jacobian
problem = NonlinearVariationalProblem(dPi, u, bcs, J=ddPi)
solver = NonlinearVariationalSolver(problem)
solver.solve()
Executing the last block is giving me error -
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/tmp/ipykernel_383/4227113402.py in <module>
20 problem = NonlinearVariationalProblem(dPi, u, bcs, J=ddPi)
21 solver = NonlinearVariationalSolver(problem)
---> 22 solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside ./dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------
This issue is arising because energy is constant and differentiating it wrt F return 0 which is similar to the issue seen in this post Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’. Reason: PETSc error code is: 63 (Argument out of range), As suggested by @dokken in the post, you might want to have a look at Hyperelasticity — DOLFIN documentation.
I hope this helps.