Solving linear system of uBLAS data types

Or, if I have the nodal values of a Function object as a numpy array, can I incorporate those values to the Function object? I wish to plot the solution using the plot() function!

On Sat, Mar 17, 2012 at 04:55:44PM -0000, Anders Pettersson wrote:
> Question #190860 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/190860
>
> Anders Pettersson gave more information on the question:
> Or, if I have the nodal values of a Function object as a numpy array,
> can I incorporate those values to the Function object? I wish to plot
> the solution using the plot() function!

Try

u.vector()[:] = x

where x is the array.

--
Anders

The problem is that the uBLASKrylovSolver does not support dense
matrices.

--
Anders

On Fri, Mar 16, 2012 at 03:30:49PM -0000, Anders Pettersson wrote:
> New question #190860 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/190860
>
> Hi!
> The program below crashes when I run it. My question is how I could solve for the unknown vector (x in Ax=b) when I have A and b as uBLASDenseMatrix and uBLASVector? What kind of object should x be? I am now trying to solve it with an uBLASKrylovSolver. The parts of interest for the question are inside the ##### ##### below:
>
> from dolfin import *
> from numpy import array, linspace
> at=Expression("1") #a(t) in the problem
> T=1 #end time
> nx=4 #number of interior cells
> kn=float(T)/nx #step size
> f=Constant("0")
> x_0=Constant("1")
>
> #defines time mesh and and normal component of the nodes
> mesh=Interval(nx, 0, T)
> normal=FacetNormal(mesh) #outer normal to nodes in this case
>
> #creates a mesh function to represent subdomains of boundary
> boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)
>
> #defines function spaces
> W_h=FunctionSpace(mesh, 'DG', 0)
> V_h=FunctionSpace(mesh, 'CG', 1)
>
> #function to determine wether x=0 or not, where the Dirichlet bc should be applied
> """def Dirichlet_boundary(x, on_boundary):
> tol=1E-9 #tolerance for round off errors
> return on_boundary and abs(x[0]<tol)"""
>
> #classes for marking different parts of the boundary
> class InitialBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[0]) < DOLFIN_EPS
>
> #creates instance of the boundary classes
> lower_boundary=InitialBoundary()
> lower_boundary.mark(boundary_parts, 0)
>
> class EndBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[0]-T) < DOLFIN_EPS
>
> upper_boundary=EndBoundary()
> upper_boundary.mark(boundary_parts, 1)
>
> #sets the Dirichlet boundary conditions on the initial boundary
> bc=DirichletBC(W_h, x_0, lower_boundary, "pointwise")
>
> #defines test and trial functions
> x=TrialFunction(W_h)
> v=TestFunction(V_h)
>
> #defines the abstract formulation
> a=(Dx(x,0) - at*x)*v*dx - jump(x)*avg(v)*dS - x*v*ds(1) + x*v*ds(0)
>
> L=f*v*dx
>
> #assembles system and connects the ds(i) symbols with A and L
> A=assemble(a, exterior_facet_domains=boundary_parts)
> b=assemble(L, exterior_facet_domains=boundary_parts)
>
> numRows = A.size(0)
> numCols = A.size(1)
>
> ###################################################
> newA=uBLASDenseMatrix(numRows+1, numRows+1)
> newA.zero()
> pos=array([1,2,3,4],dtype='I')
> newA.set(A.array(), pos, pos)
> pos2=array([0],dtype='I')
> newA.set(array([1.0]),pos2, pos2)
> pos3=array([1],dtype='I')
> pos4=array([numRows], dtype='I')
> newA.set(array([1.0]),pos4, pos4)
> newA.set(array([-1.0]),pos3, pos2)
>
> bNew=uBLASVector(numRows+1)
> bNew.zero()
> bPos=array([0], dtype='I')
> bNew.add(array([1.0]),array([0], dtype='I'))
>
> X=uBLASVector(numRows+1) #The problem might be here
>
> solver=uBLASKrylovSolver()
> solver.solve(newA, X, bNew)
>
> #####################################################
>
>
>