How can I get, save and load the external forces vector?

variational formulation
1

Hello,

I am trying to obtain the external forces force vector and stiffness matrix of this code. The stiffness matrix works but for obtaining the external forces vector I get stuck in saving the vector and how to load it afterward.

from dolfin import *
import matplotlib.pyplot as plt
import math
import numpy as np
from petsc4py import PETSc
import scipy.sparse as sp_sparse

mesh =UnitCubeMesh(4, 4, 4)
# Defining the function spaces
V_c = FunctionSpace(mesh, 'P', 1)

#STABILIZATION CROSSWIND TERMS

alpha= Constant(0.6)
I= Identity(mesh.geometry().dim())
def q(c): 
    return (sqrt(dot(dot(u,grad(c)),dot(u,grad(c))))/sqrt(dot(grad(c),grad(c))))

# Define variational problem for Picard iteration
c = TrialFunction(V_c)
# Galerkin variational problem
F=(diff*dot(grad(c),grad(v))+dot(u,grad(c))*v)*dx-f*v*dx
# Add SUPG stabilization terms
F += dot(u,grad(v))*tau_supg*dot(u,grad(c))*dx
# Add crosswind stabilization terms
CWD=(I-(outer(u,u)/sqrt(dot(dot(u,u),dot(u,u)))))
F += dot(dot(0.5*alpha*h*q(c_k)*grad(c),CWD),grad(v))*dx

# Create bilinear and linear forms
a = lhs(F)
L = rhs(F)

# Get the stifness matrix and the force vector 


A=assemble(a)
y=assemble(L)
A_mat = as_backend_type(A).mat()
y_array= as_backend_type(y).vec()

from scipy.sparse import csr_matrix
A_sparray = csr_matrix(A_mat.getValuesCSR()[::-1], shape = A_mat.size)


from scipy import sparse
sparse.save_npz("K_sparray.npz", A_sparray)

y.vector().get_local()

I don’t know how to transform the dolfin vector to numpy or neither from petcs to numpy .

2

See for instance: How to save solution in .txt format? - #8 by dokken

3

when I tried to transform it to NumPy I get the following error. That is my main problem here

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-61-2dd9543d2a7b> in <module>
----> 1 y.vector().get_local()

AttributeError: 'dolfin.cpp.la.Vector' object has no attribute 'vector'
4

Simply do y.get_local()