How to define a vector function by explicit values?

Asked by Nico Schlömer

I'd like to manually define a vector field, but I'm not sure how to do it.
I started off with

V1 = VectorFunctionSpace(mesh, 'DG', 0, 3)
J = Function(V1)
for cell_no in xrange(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    cell_midpoint = MeshEntity(mesh, 3, cell_no).midpoint()
    J.vector()[cell_no] = (1,2,3) # ???

Obviously the last thing won't work since a vector has one component only. How can I set the values of J?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Nico Schlömer
Solved:
Last query:
Last reply:

This question was reopened

Revision history for this message
Nico Schlömer (nschloe) said :
#1

Apparently, the vector is ordered such that the first N entries are the X-component, the second for Y, and the last N entries are for Z. Something like

for k in xrange(num_cells)
        J.vector()[k] = 1.0
        J.vector()[num_cells + k] = 0.0
        J.vector()[2*num_cells + k] = 0.0

will define the field (1,0,0).

Revision history for this message
Garth Wells (garth-wells) said :
#2

The above approach is not robust. It makes an assumption on the dof map that is not guaranteed (and which has changed in the dev version), and it will not work in parallel. A better approach is

class MyFunction(Expression):
    def eval(self, values, x):
        values[0] = 1.0
        values[1] = 0.0
        values[2] = 0.0
    def value_shape(self):
        return (3,)

j = MyFunction()
J.interpolate(j)

Compiled Expressions can be used to make it faster. See the DOLFIN chapter in the FEniCS book for more detail on compiled Expressions.

Revision history for this message
Nico Schlömer (nschloe) said :
#3

Is it possible to apply this logic to subdomains that are read from a file as well?
I see that the Expression class is based on geometric information whereas subdomains are really defined by cells, i.e., on the mesh and not the geometry.
There is of course closest_cell() but it would certainly not be the best thing to call every time in eval().

Revision history for this message
Anders Logg (logg) said :
#4

On Sun, Oct 07, 2012 at 11:11:15PM -0000, Nico Schlömer wrote:
> Question #210585 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210585
>
> Status: Solved => Open
>
> Nico Schlömer is still having a problem:
> Is it possible to apply this logic to subdomains that are read from
> a file as well? I see that the Expression class is based on
> geometric information whereas subdomains are really defined by
> cells, i.e., on the mesh and not the geometry. There is of course
> closest_cell() but it would certainly not be the best thing to call
> every time in eval().

You can instead overload the eval_cell function (see docstring for
Expression) and you will get an additional data argument containing
cell index and (possibly, in case you are on the boundary) facet
index.

--
Anders

Revision history for this message
Nico Schlömer (nschloe) said :
#5

Right on. Thanks!