coupling 2d and 3d physics

I think that fenics would have no issue with such a situation. One can incorporate both a surface term (an integral over ds) and bulk terms (an integral over dx) in the variational forms. There are examples that include both boundary and bulk terms in the variational forms at, for example,

http://fenicsproject.org/documentation/tutorial/fundamentals.html#multiple-dirichlet-conditions

these are simpler examples than you are considering, but one can imagine incorporating more complicate structure over the surface. There is no need to seperate it into two different problems, you can treat it as a boundary condition, albeit a complicated one.

Thank for your answer N.A. Borggren. Actually, in my problem a boundary condition is not sufficient to model what it is happening on the surface: i need to define a diffusion equation on the surface to model the spreading/absorption of a substance and couple it with the bulk diffusion. It is like having another state variable evolving on the surface, say the concentration of substance on the surface, which is assigned as a Dirichlet condition for the bulk (since the concentration on the surface is equal to the bulk concentration on the boundary). So in my model there are two concentration fields: one for the bulk and one for the balance, governed by different, coupled, diffusion equations defined on the 3d domain and on the boundary of the domain, respectively. Maybe I can add the balance equation for the surface as a weak surface term for the bulk...

Solve your 2d diffusion system for some solution g, then set boundary of 3D solution to g and solve 3D problem.

Construct the variational forms, and just download and run them, its fun. If you have any questions as to getting the forms to work ask, we will be fascinated by any well-defined physics problem that fenics can't solve.

The 2d diffusion problem is coupled with the 3d one... Maybe I should try an iterative solution technique, if I want to solve the 2d and the 3d problem sequentially. I have never used fenics, so it will take me some time to build my model. However I will be happy to share it (and ask you questions to make it work!). Thank you for the suggestions.

On Mon, Apr 16, 2012 at 05:45:49PM -0000, alucantonio wrote:
> Question #193816 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/193816
>
> Status: Answered => Solved
>
> alucantonio confirmed that the question is solved:
> The 2d diffusion problem is coupled with the 3d one... Maybe I should
> try an iterative solution technique, if I want to solve the 2d and the
> 3d problem sequentially. I have never used fenics, so it will take me
> some time to build my model. However I will be happy to share it (and
> ask you questions to make it work!). Thank you for the suggestions.

You will face two challenges if you try this with FEniCS:

1. You cannot (yet) define spaces restricted to a surface. However,
you can use a full function space and set the inactive degrees of
freedom to zero using the ident_zeros function.

2. Solving PDEs on manifolds is not supported so you will need to
map/transform the PDE manually to something that can be solved
directly with FEniCS.

--
Anders

Hi Anders, you perfectly understood my problem. I have no ideas (at the moment) about solving the second point. I have just started doing the tutorials: when I come to the full implementation of my problem, I will surely need help to solve this issue. Thank you.

Hi everybody,
after solving some test problems in fenics, I have eventually come to my actual problem, which couples rubber elasticity with solvent diffusion. The state variables of the problem are: the displacement (u,v,w), the pressure (p), the concentration of the solvent inside the gel (c) and the surface concentration field (cs). I need to introduce cs because in order to assign the concentration c on the boundary I have to solve a nonlinear algebraic equation (see the definition of fcs in the script below) to know the value of the concentration, i.e. I have an implicit boundary condition in c.
So I added two weak terms to enforce the boundary condition: inner(c-cs, d)*ds (d being the testfield associated with c), which imposes c=cs on the boundary in weak sense and inner(fcs,ee)*ds (where ee is the testfield associated with cs) to solve fcs for cs.

I get this error:

Traceback (most recent call last):
  File "/Users/eyelash/Desktop/neohooke_diffusion.py", line 84, in <module>
    chi/pow((1+Omega*Jo*cs),2))+Omega*p
TypeError: a float is required

Maybe I have chosen the wrong strategy to solve this problem: please let me know if there are other ways to assign boundary conditions by solving nonlinear algebraic equations defined on the boundary.

====BEGIN SCRIPT=====

from dolfin import *
from math import *

mesh = UnitCube(5,5,5)
V = FunctionSpace(mesh, 'Lagrange', 2)
R = FunctionSpace(mesh, 'Lagrange', 2)
Q = FunctionSpace(mesh,'Lagrange',1)
W = MixedFunctionSpace([V,V,V,Q,R,Q])

# Parameters
Omega = 1.0
D = 1E-1
Gdry = 4E5
Rg = 8.314
T = 293.0
lamo = 1.001
Jo = lamo**3
cod = (lamo**3-1.0)/Omega
co = cod/Jo
chi = 0.2
po = Gdry
muo = Rg*T*(log(Omega*cod/(1+Omega*cod),e)+1/(1+Omega*cod)+ chi/pow((1+Omega*cod),2)) \
+Omega*po
mus = muo
mutop = -10000

# Create test and trial functions
vqde = TestFunction(W) # Test function in the mixed space
v1,v2,v3,q,d,ee = split(vqde)
dupccs = TrialFunction(W) # Trial function in the mixed space

upccs = Function(W) # Function in the mixed space to store the solution
u,v,w,p,c,cs = split(upccs) # Function in each subspace to write the functional

# Definitions
F11 = 1 + u.dx(0)
F12 = u.dx(1)
F13 = u.dx(2)
F21 = v.dx(0)
F22 = 1 + v.dx(1)
F23 = v.dx(2)
F31 = w.dx(0)
F33 = 1 + w.dx(2)
F32 = w.dx(1)

J = F11*F22*F33+F12*F23*F31+F13*F21*F32-F13*F22*F31-F12*F21*F33-F11*F23*F32

FStar11 = F22*F33-F32*F23
FStar12 = F31*F23-F21*F33
FStar13 = F21*F32-F22*F31
FStar21 = F13*F32-F12*F33
FStar22 = F11*F33-F31*F13
FStar23 = F31*F12-F11*F32
FStar31 = F12*F23-F13*F22
FStar32 = F21*F13-F11*F23
FStar33 = F11*F22-F21*F12

# S = Gdry*F-p*cofac(F) # automatic differentiation not working

S11 = Gdry*F11-p*FStar11
S12 = Gdry*F12-p*FStar12
S13 = Gdry*F13-p*FStar13
S21 = Gdry*F21-p*FStar21
S22 = Gdry*F22-p*FStar22
S23 = Gdry*F23-p*FStar23
S31 = Gdry*F31-p*FStar31
S32 = Gdry*F32-p*FStar32
S33 = Gdry*F33-p*FStar33

fcs = -mus + Rg*T*(log(Omega*Jo*cs/(1+Omega*Jo*cs),e) + 1/(1+Omega*Jo*cs)+ \
chi/pow((1+Omega*Jo*cs),2))+Omega*p

a = (inner(S11,v1.dx(0)) + inner(S12,v1.dx(1)) + inner(S13,v1.dx(2)) \
+ inner(S21,v2.dx(0)) + inner(S22,v2.dx(1)) + inner(S23,v2.dx(2)) \
+ inner(S31,v3.dx(0)) + inner(S32,v3.dx(1)) + inner(S33,v3.dx(2)))*dx \
+q*(J-(1/Jo+Omega*c))*dx + inner((c-cs),d)*ds + inner(fcs,ee)*ds

def left_boundary(x, on_boundary):
 return on_boundary and abs(x[0]) < DOLFIN_EPS

def back_boundary(x, on_boundary):
 return on_boundary and abs(x[1]) < DOLFIN_EPS

def bottom_boundary(x, on_boundary):
 return on_boundary and abs(x[2]) < DOLFIN_EPS

# Boundary conditions on displacement (remove rigid body motions)
bcl = DirichletBC(W.sub(0), 0, left_boundary)
bcb = DirichletBC(W.sub(1), 0, back_boundary)
bcbo = DirichletBC(W.sub(2), 0, bottom_boundary)
bc = [bcl,bcb,bcbo]

Jac = derivative(a,upccs,dupccs)

c = interpolate(co, R)
p = interpolate(po, Q)
cs = interpolate(co, Q)

# chemical potential increment
dmu = 1000

while mus <= mutop:
 print "mus = ", mus
 problem = NonlinearVariationalProblem(a, upccs, bc, Jac)
 solver = NonlinearVariationalSolver(problem)

# set_log_level(PROGRESS)
 solver.parameters["linear_solver"]="petsc" # use PETSc builtin LU solver

 solver.solve()

 mus += dmu

=== END SCRIPT ===

I suspect the problem is importing math after dolfin. That way you
will get pow and log functions in your namespace which are not the
FEniCS (UFL) functions but standard math functions that only work with
floats (not test and trial functions).

--
Anders

On Mon, Apr 23, 2012 at 06:41:05PM -0000, alucantonio wrote:
> Question #193816 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/193816
>
> Status: Solved => Open
>
> alucantonio is still having a problem:
> Hi everybody,
> after solving some test problems in fenics, I have eventually come to my actual problem, which couples rubber elasticity with solvent diffusion. The state variables of the problem are: the displacement (u,v,w), the pressure (p), the concentration of the solvent inside the gel (c) and the surface concentration field (cs). I need to introduce cs because in order to assign the concentration c on the boundary I have to solve a nonlinear algebraic equation (see the definition of fcs in the script below) to know the value of the concentration, i.e. I have an implicit boundary condition in c.
> So I added two weak terms to enforce the boundary condition: inner(c-cs, d)*ds (d being the testfield associated with c), which imposes c=cs on the boundary in weak sense and inner(fcs,ee)*ds (where ee is the testfield associated with cs) to solve fcs for cs.
>
> I get this error:
>
> Traceback (most recent call last):
> File "/Users/eyelash/Desktop/neohooke_diffusion.py", line 84, in <module>
> chi/pow((1+Omega*Jo*cs),2))+Omega*p
> TypeError: a float is required
>
> Maybe I have chosen the wrong strategy to solve this problem: please let
> me know if there are other ways to assign boundary conditions by solving
> nonlinear algebraic equations defined on the boundary.
>
> ====BEGIN SCRIPT=====
>
> from dolfin import *
> from math import *
>
> mesh = UnitCube(5,5,5)
> V = FunctionSpace(mesh, 'Lagrange', 2)
> R = FunctionSpace(mesh, 'Lagrange', 2)
> Q = FunctionSpace(mesh,'Lagrange',1)
> W = MixedFunctionSpace([V,V,V,Q,R,Q])
>
> # Parameters
> Omega = 1.0
> D = 1E-1
> Gdry = 4E5
> Rg = 8.314
> T = 293.0
> lamo = 1.001
> Jo = lamo**3
> cod = (lamo**3-1.0)/Omega
> co = cod/Jo
> chi = 0.2
> po = Gdry
> muo = Rg*T*(log(Omega*cod/(1+Omega*cod),e)+1/(1+Omega*cod)+ chi/pow((1+Omega*cod),2)) \
> +Omega*po
> mus = muo
> mutop = -10000
>
> # Create test and trial functions
> vqde = TestFunction(W) # Test function in the mixed space
> v1,v2,v3,q,d,ee = split(vqde)
> dupccs = TrialFunction(W) # Trial function in the mixed space
>
> upccs = Function(W) # Function in the mixed space to store the solution
> u,v,w,p,c,cs = split(upccs) # Function in each subspace to write the functional
>
> # Definitions
> F11 = 1 + u.dx(0)
> F12 = u.dx(1)
> F13 = u.dx(2)
> F21 = v.dx(0)
> F22 = 1 + v.dx(1)
> F23 = v.dx(2)
> F31 = w.dx(0)
> F33 = 1 + w.dx(2)
> F32 = w.dx(1)
>
> J =
> F11*F22*F33+F12*F23*F31+F13*F21*F32-F13*F22*F31-F12*F21*F33-F11*F23*F32
>
> FStar11 = F22*F33-F32*F23
> FStar12 = F31*F23-F21*F33
> FStar13 = F21*F32-F22*F31
> FStar21 = F13*F32-F12*F33
> FStar22 = F11*F33-F31*F13
> FStar23 = F31*F12-F11*F32
> FStar31 = F12*F23-F13*F22
> FStar32 = F21*F13-F11*F23
> FStar33 = F11*F22-F21*F12
>
> # S = Gdry*F-p*cofac(F) # automatic differentiation not working
>
> S11 = Gdry*F11-p*FStar11
> S12 = Gdry*F12-p*FStar12
> S13 = Gdry*F13-p*FStar13
> S21 = Gdry*F21-p*FStar21
> S22 = Gdry*F22-p*FStar22
> S23 = Gdry*F23-p*FStar23
> S31 = Gdry*F31-p*FStar31
> S32 = Gdry*F32-p*FStar32
> S33 = Gdry*F33-p*FStar33
>
> fcs = -mus + Rg*T*(log(Omega*Jo*cs/(1+Omega*Jo*cs),e) + 1/(1+Omega*Jo*cs)+ \
> chi/pow((1+Omega*Jo*cs),2))+Omega*p
>
> a = (inner(S11,v1.dx(0)) + inner(S12,v1.dx(1)) + inner(S13,v1.dx(2)) \
> + inner(S21,v2.dx(0)) + inner(S22,v2.dx(1)) + inner(S23,v2.dx(2)) \
> + inner(S31,v3.dx(0)) + inner(S32,v3.dx(1)) + inner(S33,v3.dx(2)))*dx \
> +q*(J-(1/Jo+Omega*c))*dx + inner((c-cs),d)*ds + inner(fcs,ee)*ds
>
> def left_boundary(x, on_boundary):
> return on_boundary and abs(x[0]) < DOLFIN_EPS
>
> def back_boundary(x, on_boundary):
> return on_boundary and abs(x[1]) < DOLFIN_EPS
>
> def bottom_boundary(x, on_boundary):
> return on_boundary and abs(x[2]) < DOLFIN_EPS
>
> # Boundary conditions on displacement (remove rigid body motions)
> bcl = DirichletBC(W.sub(0), 0, left_boundary)
> bcb = DirichletBC(W.sub(1), 0, back_boundary)
> bcbo = DirichletBC(W.sub(2), 0, bottom_boundary)
> bc = [bcl,bcb,bcbo]
>
> Jac = derivative(a,upccs,dupccs)
>
> c = interpolate(co, R)
> p = interpolate(po, Q)
> cs = interpolate(co, Q)
>
> # chemical potential increment
> dmu = 1000
>
> while mus <= mutop:
> print "mus = ", mus
> problem = NonlinearVariationalProblem(a, upccs, bc, Jac)
> solver = NonlinearVariationalSolver(problem)
>
> # set_log_level(PROGRESS)
> solver.parameters["linear_solver"]="petsc" # use PETSc builtin LU solver
>
> solver.solve()
>
> mus += dmu
>
> === END SCRIPT ===
>

I added the import math line because without that I got this error:

Traceback (most recent call last):
  File "/Users/eyelash/Desktop/neohooke_diffusion.py", line 22, in <module>
    muo = Rg*T*(log(Omega*cod/(1+Omega*cod))+1/(1+Omega*cod)+ chi/pow((1+Omega*cod),2)) \
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/cpp/common.py", line 2355, in log
    return _common.log(*args)
TypeError: log expected 2 arguments, got 1

Note that the log in the muo expression contains only constants, while the log in the fcs definition contains cs which is a function. I think that the problem is the log function because if I change the expressions containing the log with other expressions, the compilation is ok and the solution process starts.
So with math import the muo line works but not fcs, while without math import the muo line does not work. I would like to understand how to use the log function correctly in fenics, both applied to floats and to fenics functions.

Use ln instead of log.

Martin
Den 24. apr. 2012 00:55 skrev "alucantonio" <
<email address hidden>> følgende:

> Question #193816 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/193816
>
> Status: Answered => Open
>
> alucantonio is still having a problem:
> I added the import math line because without that I got this error:
>
> Traceback (most recent call last):
> File "/Users/eyelash/Desktop/neohooke_diffusion.py", line 22, in <module>
> muo = Rg*T*(log(Omega*cod/(1+Omega*cod))+1/(1+Omega*cod)+
> chi/pow((1+Omega*cod),2)) \
> File
> "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/cpp/common.py",
> line 2355, in log
> return _common.log(*args)
> TypeError: log expected 2 arguments, got 1
>
> Note that the log in the muo expression contains only constants, while the
> log in the fcs definition contains cs which is a function. I think that the
> problem is the log function because if I change the expressions containing
> the log with other expressions, the compilation is ok and the solution
> process starts.
> So with math import the muo line works but not fcs, while without math
> import the muo line does not work. I would like to understand how to use
> the log function correctly in fenics, both applied to floats and to fenics
> functions.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

And in general in such a situation (namespace collisions), you can do like this to differentiate:

import math
import ufl
import dolfin
x = math.log(y)
muo = ... ufl.ln(uflexpression) ...
dolfin.log(log_message)

Thanks Martin Sandve Alnæs, that solved my question.

Using math.log for muo (float constant) and ufl.ln for fcs is ok, but now the newton solver gives NaN for both the relative and the absolute error. Maybe there is a problem with the automatic differentiation of the expression containing the log term (fcs), because when I remove the log term I get no NaNs. I would like to avoid to calculate the Jacobian manually...

=== BEGIN SCRIPT ===

from dolfin import *
import math
import ufl

mesh = UnitCube(5,5,5)
V = FunctionSpace(mesh, 'Lagrange', 2)
R = FunctionSpace(mesh, 'Lagrange', 2)
Q = FunctionSpace(mesh,'Lagrange',1)
W = MixedFunctionSpace([V,V,V,Q,R,Q])

# Parameters
Omega = 6.023e-5
D = 1E-1
Gdry = 4E5
Rg = 8.314
T = 293.0
lamo = 1.001
Jo = lamo**3
cod = (lamo**3-1.0)/Omega
co = cod/Jo
chi = 0.2
po = Gdry
muo = Rg*T*(math.log(Omega*cod/(1+Omega*cod))+1/(1+Omega*cod)+ chi/pow((1+Omega*cod),2)) \
+Omega*po
mus = muo
mutop = -10000

# Create test and trial functions
vqde = TestFunction(W) # Test function in the mixed space
v1,v2,v3,q,d,ee = split(vqde)
dupccs = TrialFunction(W) # Trial function in the mixed space

upccs = Function(W) # Function in the mixed space to store the solution
u,v,w,p,c,cs = split(upccs) # Function in each subspace to write the functional

# Definitions
F11 = 1 + u.dx(0)
F12 = u.dx(1)
F13 = u.dx(2)
F21 = v.dx(0)
F22 = 1 + v.dx(1)
F23 = v.dx(2)
F31 = w.dx(0)
F33 = 1 + w.dx(2)
F32 = w.dx(1)

J = F11*F22*F33+F12*F23*F31+F13*F21*F32-F13*F22*F31-F12*F21*F33-F11*F23*F32

FStar11 = F22*F33-F32*F23
FStar12 = F31*F23-F21*F33
FStar13 = F21*F32-F22*F31
FStar21 = F13*F32-F12*F33
FStar22 = F11*F33-F31*F13
FStar23 = F31*F12-F11*F32
FStar31 = F12*F23-F13*F22
FStar32 = F21*F13-F11*F23
FStar33 = F11*F22-F21*F12

S11 = Gdry*F11-p*FStar11
S12 = Gdry*F12-p*FStar12
S13 = Gdry*F13-p*FStar13
S21 = Gdry*F21-p*FStar21
S22 = Gdry*F22-p*FStar22
S23 = Gdry*F23-p*FStar23
S31 = Gdry*F31-p*FStar31
S32 = Gdry*F32-p*FStar32
S33 = Gdry*F33-p*FStar33

fcs = -mus + Rg*T*(ufl.ln(Omega*Jo*cs/(1+Omega*Jo*cs)) + 1/(1+Omega*Jo*cs)+ \
chi/pow((1+Omega*Jo*cs),2))+Omega*p

a = (inner(S11,v1.dx(0)) + inner(S12,v1.dx(1)) + inner(S13,v1.dx(2)) \
+ inner(S21,v2.dx(0)) + inner(S22,v2.dx(1)) + inner(S23,v2.dx(2)) \
+ inner(S31,v3.dx(0)) + inner(S32,v3.dx(1)) + inner(S33,v3.dx(2)))*dx \
+q*(J-(1/Jo+Omega*c))*dx + inner(fcs,ee)*ds + inner((c-cs),d)*ds

def left_boundary(x, on_boundary):
 return on_boundary and abs(x[0]) < DOLFIN_EPS

def back_boundary(x, on_boundary):
 return on_boundary and abs(x[1]) < DOLFIN_EPS

def bottom_boundary(x, on_boundary):
 return on_boundary and abs(x[2]) < DOLFIN_EPS

bcl = DirichletBC(W.sub(0), 0, left_boundary)
bcb = DirichletBC(W.sub(1), 0, back_boundary)
bcbo = DirichletBC(W.sub(2), 0, bottom_boundary)
bc = [bcl,bcb,bcbo]

Jac = derivative(a,upccs,dupccs)

w = Function(V)
cc = Function(R)

dmu = 1000

c = interpolate(Constant(co), R)
p = interpolate(Constant(po), Q)
cs = interpolate(Constant(co), Q)

while mus <= mutop:
 print "mus = ", mus
 problem = NonlinearVariationalProblem(a, upccs, bc, Jac, form_compiler_parameters={})
 solver = NonlinearVariationalSolver(problem)

 solver.parameters["linear_solver"]="petsc" # use PETSc builtin LU solver

 solver.solve()

 mus += dmu

=== END SCRIPT ===

A part from the NaNs, when the ln term is present in the expression of mu, I get this warning:

WARNING: The number of integration points for each cell will be: 125
           Consider using the option 'quadrature_degree' to reduce the number of points

Please anybody can explain to me what this means?

On 04/29/2012 03:10 PM, al6 wrote:
> Question #193816 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/193816
>
> al6 gave more information on the question:
> A part from the NaNs, when the ln term is present in the expression of
> mu, I get this warning:
>
> WARNING: The number of integration points for each cell will be: 125
> Consider using the option 'quadrature_degree' to reduce the number of points
>
> Please anybody can explain to me what this means?

FFC is conservative with the number of quadrature points automatically
used for the generated code. This can be controlled by setting:

parameters.form_compiler.quadrature_degree = 2 # or some higher number

Johan