This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

How to define a vector function depends on the solution?

0 votes

Basically, referring the hyper-elasticity example, I wanna add a body force B which is a 2D vector in the 2d domain. Its components depend on the displacement. For example, let's say

body force

   B = [bx, by]
   bx = 1/(x+ux)
   by = 1/(x+uy)

In the python code, I write

Define functions

  V = VectorFunctionSpace(mesh, "Lagrange", 1)
  du = TrialFunction(V)            # Incremental displacement
  v  = TestFunction(V)             # Test function
  u  = Function(V)                 # Displacement from previous iteration
  ...
  def B(u):
      return Expression(("1/(x[0]+u[0])","1/(x[1]+u[1])"))
  ...

  # Total potential energy
  Pi = psi*dx - dot(B(u), u)*dx - dot(T, u)*ds

and I get the error, which says the way of defining the body force B is not correct:

In instant.recompile: The module did not compile with command 'make VERBOSE=1', see '/home/abc/.instant/error/dolfin_compile_code_87e04d513466c00e045d483572e6d4f47c45a4a0/compile.log'
Traceback (most recent call last):
File "hyper_elast.py", line 95, in
Pi = psidx - dot(B(u), u)dx - dot(T, u)ds
File "hyper_elast.py", line 83, in B
return Expression(("1/(x[0]+u[0])","1/(x1+u1)"))
File "/usr/lib/python2.7/dist-packages/dolfin/functions/expression.py", line 602, in new
mpi_comm=mpi_comm)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/expressions.py", line 217, in compile_expressions
mpi_comm=mpi_comm)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/expressions.py", line 145, in compile_expression_code
mpi_comm=mpi_comm)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 64, in mpi_jit
return local_jit(
args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/compilemodule.py", line 458, in compile_extension_module
**instant_kwargs)
File "/usr/lib/python2.7/dist-packages/instant/build.py", line 563, in build_module
recompile(modulename, module_path, new_compilation_checksum, build_system)
File "/usr/lib/python2.7/dist-packages/instant/build.py", line 165, in recompile
instant_error(msg % (cmd, compile_log_filename_dest))
File "/usr/lib/python2.7/dist-packages/instant/output.py", line 85, in instant_error
raise RuntimeError(text)
RuntimeError: In instant.recompile: The module did not compile with command 'make VERBOSE=1', see '/home/abc/.instant/error/dolfin_compile_code_87e04d513466c00e045d483572e6d4f47c45a4a0/compile.log'

Any one can help me with this? Thank you a lot.

asked May 30, 2016 by Wilhelm FEniCS Novice (410 points)

Did you see /home/abc/.instant/error/dolfin_compile_code_87e04d513466c00e045d483572e6d4f47c45a4a0/compile.log as the error message told you?

1 Answer

+1 vote
 
Best answer

I think

Expression(("1/(x[0]+u0)","1/(x[1]+u1)"),u0=u.sub(0),u1=u.sub(1))

might work.

answered May 30, 2016 by KristianE FEniCS Expert (12,900 points)
selected Jun 1, 2016 by Wilhelm

Thanks a lot. That works for me. A further issue is that, based on the current code, I have a conditional statement for the body force dependent on the current coordinate of the particle , for example:

body force

B = [b_x, b_y]
if y+u_y > 0.5:
    b_x = 1/(x+u_x)
    b_y = 1/(y+u_y)
else:
    b_x = 2/(x+u_x)
    b_y = 2/(y+u_y)

In the python code, I have

python

V = VectorFunctionSpace(mesh, "Lagrange", 1)
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
...
def B(u,x):
    u0=u.sub(0)
    u1=u.sub(1)
    if x[1]+u1 > 0.5:
        return Expression(("1/(x[0]+u[0])","1/(x[1]+u[1])"))
    else:
        return Expression(("2/(x[0]+u[0])","2/(x[1]+u[1])"))
...

# Total potential energy
Pi = psi*dx - dot(B(u,x), u)*dx - dot(T, u)*ds

Run this I get the following error:

UFL conditions cannot be evaluated as bool in a Python context.
Traceback (most recent call last):
File "linear_elast.py", line 87, in
Pi = psidx - dot(B(u,x), u)dx - dot(T, u)ds
File "linear_elast.py", line 59, in B
if x[1] > 0.5:
File "/usr/lib/python2.7/dist-packages/ufl/conditional.py", line 42, in bool
error("UFL conditions cannot be evaluated as bool in a Python context.")
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(
message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

which is related with the conditional statement. Can you help me with this? Thank you.

You can do

Expression(("x[1] + u1 ? 1/(x[0]+u0) : 2./(x[0]+u0) ","x[1] + u1 ? 2./(x[1]+u1)"),u0=u.sub(0),u1=u.sub(1))

Thank you KristantE. But where is the condition x[1]+u1 > 0.5? Should it be like this:

Expression(("x[1] + u1<0.5 ? 1/(x[0]+u0) : 2./(x[0]+u0) ","x[1] + u1<0.5 ? 2./(x[1]+u1)"),u0=u.sub(0),u1=u.sub(1))

?

No, like

Expression(("x[1] + u1 > 0.5 ? 1/(x[0]+u0) : 2./(x[0]+u0) ","x[1] + u1 > 0.5 ? 2./(x[1]+u1)"),u0=u.sub(0),u1=u.sub(1))

You are right. Thank you so much!!

I am trying to do a similar thing, but my u is a scalar function of 2 variables. I have something like

Expression(("1/(x[0]+z)","1/(x[1]+z)"),z=u)

but I get the error

The module did not compile with command 'make VERBOSE=1'

What should I do differently?

You should post a minimum working example.

I figured it out, it was because I had z**2 instead of pow(z,2).

...