Hi all, I’m struggling with fenics because of an attribute error. I have already tried with the solution proposed in a similar problem (AttributeError (object has no attribute '_cpp_object')) but it didn’t work.
My problem occurs when I try to save a function (u) in an xdmf file. This procedure works fine with the mesh and the boundary savings, so I can’t understard why it makes problems…
Here my code, thank for help!
`from future import division
from dolfin import *
import math
import mshr
import numpy as np
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
}
}
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1.)
class Left(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],0.)
class Bottom(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],0.)
class deformation:
def __init__(self,**kwargs):
super().__init__(**kwargs)
self.mesh_xdmf = XDMFFile("xdmf_files/mesh.xdmf")
self.boundary_xdmf = XDMFFile("xdmf_files/boundary.xdmf")
self.save_file_xdmf = XDMFFile("disp.xdmf")
self.save_file = File("output/disp.pvd")
self.mesh = self.define_mesh()
self.functionSpace()
self.initialize_functions()
self.gamma = Constant(0)
self.mu = Constant(1)
self.k = Constant(20)
self.bcs = self.boundary_conditions(self.V)
def define_mesh(self):
mesh = UnitSquareMesh(10,10)
encoding = XDMFFile.Encoding.HDF5
self.mesh_xdmf.write(mesh, encoding)
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)
print('Dofs: {}'.format(self.V.dim()))
boundaries = self.subdomains(self.V)
encoding = XDMFFile.Encoding.HDF5
self.boundary_xdmf.write(boundaries,encoding)
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):
left = Left()
bottom = Bottom()
right = Right()
boundaries = MeshFunction("size_t",self.mesh,1)
boundaries.set_all(0)
left.mark(boundaries,1)
bottom.mark(boundaries,2)
right.mark(boundaries,3)
return boundaries
def boundary_conditions(self,V):
boundaries = self.subdomains(V)
bcb1 = DirichletBC(V.sub(0).sub(0), 0, boundaries,2)
bcb2 = DirichletBC(V.sub(0).sub(1), 0, boundaries,2)
return [bcb1,bcb2]
def FirstAndSecondVariation(self):
boundaries = self.subdomains(self.V)
self.ds = Measure("ds",domain=self.mesh,subdomain_data=boundaries)
X = SpatialCoordinate(self.mesh)
x = X + self.u
I = Identity(2)
F = grad(x)
C = F.T*F
I1 = tr(C)
self.J = det(F)
self.con_neu = as_vector([self.gamma, 0.0])
psi = self.mu/2. * (I1 - 2) - self.mu*ln(self.J) + self.k*ln(self.J)**2
self.W = psi*dx + inner(self.u, self.con_neu)*ds(3) + inner(self.u, self.con_neu)*ds(1)
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 = deformation()
T = 1.2
t = 0.0
dt = 0.004
dt1 = -0.01
gamma0 = 0.0
gamma_max = 1.
test = 0
while t < T and dt>1e-6:
t += dt
my_problem.gamma.assign(gamma0+gamma_max*t)
print("t: ",t)
print("T: ",float(my_problem.gamma))
ok = 0
while ok ==0:
try:
print("In try...")
my_problem.newton_solver()
my_problem.up0.assign(my_problem.up)
my_problem.monitor()
ok = 1
except RuntimeError:
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)
if dt < 1e-6:
print("Error: dt < 1e-6")
print("Finished")
encoding = XDMFFile.Encoding.HDF5
my_problem.save_file_xdmf.write(my_problem.u)
`