Hi all, I’m dealing with a plasto-elastic problem. I wrote my dolfin code and at the end I’d like to save the mesh in order to start a new calculation with that result. Due to the complexity, I thought that I should use the format XDMF, I red a lot of documentation and I thougth that my code was ok… But it turned out this error: "AttributeError: module ‘dolfin.mesh’ has no attribute ‘_cpp_object’ ". Can anyone please help me?
Here there is my code:
from __future__ import division
from dolfin import *
from ufl import diag
import h5py
import mshr
import numpy as np
import math
from mesh_cyl import mesh_cyl
from mpi4py import MPI
import os
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["allow_extrapolation"] = True
parameters["form_compiler"]['quadrature_degree'] = 6
params = {'nonlinear_solver': 'snes',
'snes_solver':
{
'linear_solver' : 'mumps',
'absolute_tolerance' : 1e-10,
'relative_tolerance' : 1e-10,
'maximum_iterations' : 20,
# 'error_on_nonconvergence' : False
}
}
class Top(SubDomain):
def inside(self, x, on_boundary):
TOL = 0.001
return on_boundary and near(x[2], 7., TOL)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
TOL = 0.001
return on_boundary and near(x[2], 0., TOL)
class torsion_cyl:
def __init__(self):
self.mesh = self.define_mesh()
self.functionSpace()
self.initialize_functions()
self.gamma = Constant(0)
self.mu = 1
self.bcs = self.boundary_conditions(self.V)
def define_mesh(self):
self.L = 7.
k = (2*DOLFIN_PI)/self.L
m = 1
self.Ri = 0.35
self.Ro = 1
mesh = mesh_cyl(self.Ri, self.Ro, self.L, 20)
return mesh
def functionSpace(self):
Velem = VectorElement("Lagrange", self.mesh.ufl_cell(), 2)
Pelem = FiniteElement("Lagrange", self.mesh.ufl_cell(), 1)
VPRelem = MixedElement([Velem, Pelem])
self.V = FunctionSpace(self.mesh, VPRelem)
boundaries = self.subdomains(self.V)
def initialize_functions(self):
self.up = Function(self.V)
self.u, self.p = split(self.up)
self.up0 = Function(self.V)
def subdomains(self, V):
top = Top()
bottom = Bottom()
boundaries = MeshFunction("size_t", self.mesh, 2)
boundaries.set_all(0)
top.mark(boundaries, 1)
bottom.mark(boundaries, 2)
return boundaries
def boundary_conditions(self, V):
boundaries = self.subdomains(V)
self.ds = Measure("ds", domain=self.mesh, subdomain_data=boundaries)
bcb1 = DirichletBC(V.sub(0).sub(0), 0, boundaries, 1)
bcb2 = DirichletBC(V.sub(0).sub(1), 0, boundaries, 1)
bcb3 = DirichletBC(V.sub(0).sub(2), 0, boundaries, 1)
bct1 = DirichletBC(V.sub(0).sub(0), 0, boundaries, 2)
bct2 = DirichletBC(V.sub(0).sub(1), 0, boundaries, 2)
bct3 = DirichletBC(V.sub(0).sub(2), 0, boundaries, 2)
return [bcb1, bcb2, bcb3, bct1, bct2, bct3]
def FirstAndSecondVariation(self):
boundaries = self.subdomains(self.V)
X = SpatialCoordinate(self.mesh)
R,Theta,Z = X
x = X + self.u
r,theta,z = x
G = as_tensor([[1.,0.,0.],[0.,1.,-self.gamma*R],[0.,0.,1.]])
Ge = det(G)
invG = inv(G)
F = grad(x)
invF = inv(F)
Fe = F * invG
Ce = Fe.T * Fe
self.Je = det(Fe)
I1 = tr(Ce)*self.Je**(-2./3.)
Iden = diag(as_vector([1,1,1]))
psi = (self.mu/2.*(I1-3)+Constant(100)*ln(self.Je)**2-self.p*(self.Je-1))
self.W = psi*dx
FF = derivative(self.W, self.up, TestFunction(self.V))
dF = derivative(FF, self.up, TrialFunction(self.V))
return FF, dF
def monitor(self):
u = Function(self.up,0,name="displacement")
p = Function(self.up,1,name="pressure")
def newton_solver(self):
FF, dF = self.FirstAndSecondVariation()
problem = NonlinearVariationalProblem(FF, self.up, self.bcs, dF)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(params)
solver.solve()
my_problem = torsion_cyl()
T = 1.2 # total simulation time
t = 0.0
dt = 0.004
dt1 = -0.01
gamma0 = 0.0
gamma_max = 1.
test = 0
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
while t < T and dt>1e-6:
t += dt
my_problem.gamma.assign(gamma0 + gamma_max*t)
if rank == 0:
print(f"t: {t}")
print(f"T: {float(my_problem.gamma)}")
ok = 0
while ok == 0 :
try:
if rank == 0:
print("In try")
my_problem.newton_solver()
my_problem.up0.assign(my_problem.up)
my_problem.monitor()
ok = 1
except RuntimeError:
if rank == 0:
print("Error")
dt = dt/2.
t += -dt
my_problem.gamma.assign(gamma0 + gamma_max*t)
print(float(my_problem.gamma))
my_problem.up.assign(my_problem.up0)
encoding = XDMFFile.Encoding.HDF5
file = XDMFFile("mesh.xdmf")
file.write(mesh,encoding)