Variable vector field

Question

I am wondering if it is possible to define a variable vector field. I give the most important code below:

mesh = UnitCubeMesh(10,10,10) 

V1 = VectorFunctionSpace(mesh,'P',2) 

class Conductor(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[2], (0.4, 0.6)) and between(x[1], (0.4, 0.6)) and between(x[0], (0.4, 0.6)))

conductor = Conductor()

domain_marker =  CellFunction("size_t", mesh) 
domain_marker.set_all(0)
conductor.mark(domain_marker, 1)
dx = Measure('dx', domain=mesh, subdomain_data=domain_marker)

f1 = Expression(('sin(x[0])','sin(x[1])','sin(x[2])')), degree=4)
f2 = Expression(('sin(x[2])','sin(x[0])','sin(x[1])')), degree=4)

class Source(Expression):
    def __init__(self,domain_marker, f1, f2, **kwargs):
        self.domain_marker = domain_marker
        self.f1 = f1
        self.f2 = f2

    def eval_cell(self, values, x, cell):
        if self.domain_marker[cell.index] == 0:
            values[0] =  self.f1(x)
        else:
            values[0] =  self.f2(x)

f = Source(domain_marker,f1,f2)  

plot(interpolate(f,V1),interactive = True)

This is working for scalar functions, but not for vector fields. I obtain the following error: Rank of function (0) does not match rank of function space (1). I appreciate any help.

Answer 1

Hi, the way you have defined Source it is a scalar expression, see here. The code needs to be modified as follows

f1 = Expression(('sin(x[0])','sin(x[1])','sin(x[2])'), degree=4)
f2 = Expression(('sin(x[2])','sin(x[0])','sin(x[1])'), degree=4)

class Source(Expression):
    def __init__(self,domain_marker, f1, f2, **kwargs):
        self.domain_marker = domain_marker
        self.f1 = f1
        self.f2 = f2

    def eval_cell(self, values, x, cell):
        if self.domain_marker[cell.index] == 0:
            values[:] =  self.f1(x)
        else:
            values[:] =  self.f2(x)

    def value_shape(self): return (3, ) 

f = Source(domain_marker,f1,f2, degree=4) 

Comments

Thank you both for the solution.

Answer 2

If your custom Expression is not a scalar function, you need to define it's shape by overriding Expression.value_shape(). E.g.:

class Source(Expression):
    def __init__(self,domain_marker, f1, f2, **kwargs):
        self.domain_marker = domain_marker
        self.f1 = f1
        self.f2 = f2

    def eval_cell(self, values, x, cell):
        if self.domain_marker[cell.index] == 0:
            values[0] =  self.f1(x)[0]
        else:
            values[0] =  self.f2(x)[0]

    def value_shape(self):
        return (3,)

Note also that f1 and f2 are vector valued functions, and therefore you need to populate the values array correctly. For brevity I simply chose the first element.