AttributeError: 'Sum' object has no attribute 'vector'

Question

Hello,

I would like to solve a second order hyperbolic equation using Newmark method for the time integration on a 2D domain.

I would like to save both the solution and its second derivative (both .csv and .pvd files). How can I do this?
- For the .csv I get an error saying that u1 has no attribute 'vector' as it is a sum of different objects.
- For the .pvd, do I need to project the solution over the functionspace? I still get troubles when I open the time series in paraview.

Here is the code I'm using:

from dolfin import *
from mshr import *
import numpy  as np


# Create mesh
size = 10
ylim = pi 
pdim = 2

domain = Rectangle(Point(0., 0.), Point(ylim, ylim))
mesh = generate_mesh(domain,size)


# Newmark Constants
beta_ = 0.35
gamma_ = 0.5

# Time Constants
dt =  1e-5
T = 5e-5
t = 0.


# Boundaries
class Horizontal(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - ylim) < DOLFIN_EPS))

class Vertical(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[0] - 0.) < DOLFIN_EPS) or (abs(x[0] - ylim) < DOLFIN_EPS))


boundaries = FacetFunction("size_t", mesh, 0) 

Horizontal().mark(boundaries, 1) 
Vertical().mark(boundaries, 2) 


# Define Function space, trial and test functions
V = FunctionSpace(mesh, 'CG', pdim)
a, phi = TrialFunction(V), TestFunction(V)

# Define Dirichlet BC
bcV = DirichletBC(V, Constant(0.), boundaries, 2)


# Fields from previous time step (displacement, velocity, acceleration)
u0IC = Expression('sin(x[0])', degree = pdim)
v0IC = Expression('-sin(x[0])', degree = pdim)

u0 = interpolate(u0IC, V)
v0 = interpolate(v0IC, V)


# Define time-dependent source
c = Constant(0.2) # velocity
source = Expression('( sin(x[0])/(1.+t) + (2.*sin(x[0]))/( pow(c,2)*(pow(1.+t,3)) ) )', c = c, t = t, degree = pdim)


# Blocks definition
def k_block(value):
    return inner(grad(value),grad(phi)) * dx 


m = 1/pow(c,2) * a * phi * dx 
f = source * phi * dx

# Assemble left hand side 
A = assemble(m + beta_*pow(dt,2)*k_block(a))
# Assemble right hand side
b0 = assemble(f-k_block(u0))
# Apply BC
bcV.apply(A,b0) 

# Compute a0
a0 = Function(V)
solve (A,a0.vector(),b0)


# Save initial displacement (actually pressure here) and acceleration  
vtk_file_a = File("results/p"+str(pdim)+"_"+str(size)+"_acc.pvd")
#vtk_file_a << a0

# vtk_file_u = File("results/p"+str(pdim)+"_"+str(size)+"_displacement.pvd")
# vtk_file_u << u0

# Define right hand side in general (Newmark method)
L = f - (k_block(u0) + dt*k_block(v0) + pow(dt,2)*(0.5-beta_)*k_block(a0))

# Initialize a1 
a1 = Function(V)

# Time loop
while t <= T:

    # Update time 
    t += dt
    print "t="+str(t)

    # Update the source that is time depedent and save it
    source.t = t

    # Assemble the RHS & apply BCs
    b = assemble(L)
    bcV.apply(A,b) 
    # Solve system for displacement
    solve(A, a1.vector(), b)


    # Newmark method to compute displacement and velocity
    u1 = u0 + dt*v0 + pow(dt,2)*((0.5-beta_)*a0 + beta_*a1)
    v1 = v0 + dt*((1-gamma_)*a0 + gamma_*a1)

    # Save solution to VTK format
    vtk_file_a << a1
    # vtk_file_u << project(u1,V)

    # Update solution, velocity and acceleration:
    u0.assign(u1)
    v0.assign(v1)
    a0.assign(a1)


X = V.tabulate_dof_coordinates()
X.resize((V.dim(), 2))
x = X[:,0]
y = X[:,1]
a_values  = a1.vector().array()
u_values  = u1.vector().array()
np.savetxt("ex_res/a_p"+str(pdim)+"_"+str(size)+".csv", np.c_[x, y, a_values], delimiter=",")
np.savetxt("ex_res/u_p"+str(pdim)+"_"+str(size)+".csv", np.c_[x, y, u_values], delimiter=",")

Answer 1

u1 is currently a ufl expression. You need to project it onto an FE space. Consider the following:

u1_proj = project(u1, V)
u_values  = u1_proj.vector().array()

Comments

Thank you nate!

This solves the problem for saving the csv, but it is still not possible to open the time series of the solution in paraview: for each time step I get different "types" f_69, f_123, ... so it is impossible to have a "smooth video" of the solution.

Also, how precise is project? Can I expect the solution u1 to have the same accuracy of a1?

---

To solve the problem related to time-series try something like:

u1_proj = Function(V)
while t <= T:
  .
  .
  .
  u1_proj.vector()[:] = project(u1, V).vector()

Also, if you don't wan't to project you can do this:

# Extract vector components
u0_vec = u0.vector()
v0_vec = v0.vector()
a0_vec = a0.vector()
a1_vec = a0.vector()

# Calculate displacement and velocity
u1.vector()[:] = u0_vec + dt*v0_vec + pow(dt,2)*((0.5-beta_)*a0_vec + beta_*a1_vec)
v1.vector()[:] = v0_vec + dt*((1-gamma_)*a0_vec + gamma_*a1_vec)

---

Thank you hernan_mella! It works just fine ;)