This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

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

+2 votes

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=",")
asked May 18, 2017 by caterinabig FEniCS User (1,460 points)

1 Answer

+2 votes
 
Best answer

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()
answered May 18, 2017 by nate FEniCS Expert (17,050 points)
selected May 19, 2017 by caterinabig

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 ;)

...