Using eval function to define a vector

Question

Hello
I want to define a 2*1 vector inside a class and use the eval function to define the components of this vector, like:

class BodyForce(Expression):

    def __init__(self, **kwargs):

        self.value = kwargs['value']

    def eval(self, value,x):

        value[0] = self.value[0]
        value[1] = self.value[1]

    def value_shape(self):
        return (2,)

GG = BodyForce(value=(Constant(A)*(c - 1.), 0.0), degree=1)

The thing is the first component of this vector is with respect to "c" which is an unknown that is needed to be calculated. In addition this vector should be applied in all over the domain. I have to use this vector in my variational problem. The contributing term from this vector in the variational problem is:

inner(GG,v_3) * dx

I get this error:

raceback (most recent call last):
  File "/home/hp/Desktop/ert/23.py", line 171, in <module>
    solve(F == 0, z, bcs)
  File "/home/hp/Desktop/ert/23.py", line 151, in eval
    value[0] = self.value[0]
ValueError: setting an array element with a sequence.

Any idea how to fix it?
Thanks!

Answer 1

I'm not sure what you want to achieve, but it seems that you want to define a constant coefficient (A (c - 1), 0), where c varies, I guess in an iterative method.

If you want to define an expression which depends on some coefficient and you want to update this coefficient, you can just access it as a property of the object:

mesh = UnitSquareMesh(1,1)
V = FunctionSpace(mesh, 'Lagrange', 1)
v = TestFunction(V)
GG = Expression(("A*(c-1)","0"), A=7.0, c=42.0, degree=1)
a = inner(GG, grad(v))*dx
M1 = assemble(a)
GG.A = 14.0
GG.c = -7.0
M2 = assemble(a)
(M1-M2).array()  # we did something

Comments

Hello
Thanks for your response but the thing that I am looking for is different.
As I explained I want to define a Vector. The first component of this vector is : A(c-1) where "A" is constant but "c" is unknown and is present in the variational form. The second component of this vector is zero.
I have already solved the problem by defining this vector in another way :

# Define vectorial body force
def F():
    return as_vector([Constant(A)*(c - 1.)] +[0.0])

I just wanna know if it is possible to define the "F" in this fashion:

class BodyForce(Expression):

    def __init__(self, **kwargs):

        self.value = kwargs['value']

    def eval(self, value,x):

        value[0] = self.value[0]
        value[1] = self.value[1]

    def value_shape(self):
        return (2,)

F = BodyForce(value=(Constant(A)*(c - 1.), 0.0), degree=1)

However I could not get it to work!
I want it to play the same role as a vector like F() as I mentioned above. I need to mention that I got this idea from your answer in :

https://fenicsproject.org/qa/13149/clamped-beam-under-the-end-load-problem

In other words I want to know if it is possible to define a vector by combining:

1. class
2. def init(self, **kwargs),
3. def eval(self, value,x).

Thanks in advance for your time.

---

I see... well the approach with eval in the class cannot work, or at least not in that form, since value is a numpy array, and you cannot assign Constant(A)*(c-1) to its first component. What about this:

class BodyForce(Expression):

    def __init__(self, **kwargs):
        self.A = kwargs['A']
        self.c = kwargs['c']
        self.x1 = kwargs['x1']

        assert iscallable(self.c), 'I need to call c(x)'
        assert isintance(self.A, float) and isinstance(self.x1, float), 'I need scalar constants'

    def eval(self, value,x):
        value[0] = self.A * (self.c(x)-1)
        value[1] = self.x1

    def value_shape(self):
        return (2,)

And use it with

GG = BodyForce(A=A, c=c, x1=0, degree=1)

---

I changed what you suggested a little bit. I think there was a typo in there. So this is what I replaced in my code:

class BodyForce(Expression):

    def __init__(self, **kwargs):
        self.A = kwargs['A']
        self.c = kwargs['c']
        self.x1 = kwargs['x1']

        assert callable(self.c)
        assert isinstance(self.A, float) and isinstance(self.x1, float)

    def eval(self, value,x):
        value[0] = self.A * (self.c(x)-1)
        value[1] = self.x1

    def value_shape(self):
        return (2,)

But then I got a new error:

assert isinstance(self.A, float) and isinstance(self.x1, float)
AssertionError

Then I removed the lines :

assert callable(self.c)
assert isinstance(self.A, float) and isinstance(self.x1, float)

This time it worked but it increased the time of running dramatically. I guess it is not a right approach to define a vector this way. Anyway, thanks again for providing the support.

---

My pleasure. The assertions were there just to make it clear that we wanted to be multiplying numbers, not dolfin Constants or stuff like that. And yes, I can see it taking a lot longer, being pure python instead of whatever it is that dolfin does under the hood. Happy hacking!