DOLFIN

Overwrite vector in assembly

Asked by Marcos Samudio

Hello everyone!

I am trying to overwrite a previously computed vector, say u, and assign some new values to it. My previous question is related to this same subject, but I think this example might provide further understanding!

from dolfin import *

mesh = UnitSquare(5, 5)
mf = FacetFunction("uint", mesh)
mf.set_all(0)
V = VectorFunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
u = Function(V)

#We set u to some value, just so it is initialized
u.vector()[:] = 1.

class DomainBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
DomainBoundary().mark(mf, 1)

#Say I want to set all values on the boundary to (0.2, 0.), I was thinking on doing:
f = Constant((0.2,0.))
L = inner(f, v)*ds(1)
assemble(L, tensor=u.vector(), exterior_facet_domains=mf, reset_sparsity=True)
File('u_end.pvd') << u

But I do not get the expected result... Is there another way I could achieve this?
Note that this is a really simple example, and f would be generic in my case. In fact, it depends on u.

Thanks,

Marcos

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
Last reply:
Revision history for this message
Best Mikael Mortensen (mikael-mortensen) said : 		#1

Hi Marcos,

Den Apr 23, 2012 kl. 9:45 PM skrev Marcos Samudio:

> New question #194519 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/194519
>
> Hello everyone!
>
> I am trying to overwrite a previously computed vector, say u, and assign some new values to it. My previous question is related to this same subject, but I think this example might provide further understanding!
>
> from dolfin import *
>
> mesh = UnitSquare(5, 5)
> mf = FacetFunction("uint", mesh)
> mf.set_all(0)
> V = VectorFunctionSpace(mesh, 'CG', 1)
> v = TestFunction(V)
> u = Function(V)
>
> #We set u to some value, just so it is initialized
> u.vector()[:] = 1.
>
> class DomainBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary
> DomainBoundary().mark(mf, 1)
>
> #Say I want to set all values on the boundary to (0.2, 0.), I was thinking on doing:
> f = Constant((0.2,0.))
> L = inner(f, v)*ds(1)
> assemble(L, tensor=u.vector(), exterior_facet_domains=mf, reset_sparsity=True)
> File('u_end.pvd') << u
>
> But I do not get the expected result... Is there another way I could achieve this?

You could create a DirichletBC and apply it to u.vector(). I think it's something like

bnd = DirichletBC(V, f, mf, 1)
bnd.apply(u.vector())

Mikael

> Note that this is a really simple example, and f would be generic in my case. In fact, it depends on u.
>
> Thanks,
>
> Marcos
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Revision history for this message
Marcos Samudio (marcossamudio) said : 		#2

Hello Mikael!

Thanks for the fast response.

Indeed a Dirichlet BC seems to be the way to go.

Still, I have to relate boundary nodes to inner nodes, so I have something like:

u_j = bound_f(u_i), where j is the boundary node and i is its inner related node.

I tried the following:

        for i in vertices_on_boundary:
            j = boundary_to_inner[i] #this is the mapping
            for k in range(mesh.geometry().dim()):
                u.vector()[i+k*N] = u.vector()[j+k*N] #N is the number of nodes

But I get an error because u is a function that comes out from splitting a MixedFunctionSpace function, as I mentioned in my question https://answers.launchpad.net/dolfin/+question/194259

If this worked, I would simply do:

DirichletBC(V, bound_f, mf, bid)

and forget about indices...

I will keep trying, but if you spot anything, please let me know.

Revision history for this message
Marcos Samudio (marcossamudio) said : 		#3

Thanks Mikael Mortensen, that solved my question.