Dear all,
Firstly, I must say FEniCS is really an awesome software. Thanks the developer!!!
I'm using the new version FEniCS, but there are few tutorials about the newest version. Somethings, especially for the coupled initial value problems with multiple dofs for each node.
I see the tutorial on Github for the reaction system(which will deal with multiple dofs, but that doesn't help me)
For example, for some simple diffusion problem coupled with poisson equation.
Let's consider that there are three species,$c_{1}$,$c_{2}$ and $c_{3}$. They are coupled with a potential $\phi$
The governing equations for these four dofs(c1,c2,c3,phi) are list as below:
$\frac{\partial c_{1}}{\partial t}=\nabla\cdot(D_{1}\nabla c_{1}+D_{1}c_{1}\nabla\Phi)$
$\frac{\partial c_{2}}{\partial t}=\nabla\cdot(D_{2}\nabla c_{2}+\frac{D_{2}c_{2}}{1-c_{3}}\nabla c_{3}+\frac{D_{2}c_{2}}{1-c_{3}}\nabla\Phi)$
$\frac{\partial c_{3}}{\partial t}=\nabla\cdot(D_{3}\nabla c_{3}+\frac{D_{3}c_{3}}{1-c_{2}}\nabla c_{2}+\frac{D_{3}c_{3}}{1-c_{2}}\nabla\Phi)$
$\sigma\nabla^{2}\Phi+\sum_{i=1}^{3}c_{i}=0$
In the old version FEniCS, I defined the function space as follow:
# Mesh info
w=10;h=0.1*w
mesh=RectangleMesh(Point(0,0),Point(w,h),100,10)
P1=FiniteElement('P',triangle,1)
element=MixedElement([P1,P1,P1,P1]))
# Is it right that how many dofs per node your are using,how many terms
# will be used in the MixedElement?
# Function space
V=FunctionSpace(mesh,element)
du=TrialFunction(V)
v1,v2,v3,v4=TestFunctions(V)
u=Function(V)
u0=Function(V)
# In the old version, I defined the initial class as follow to compute the initial
# value, but it dosen't work in the newest version. Why?
# Initial Condition
class InitialConditions(Expression):
def eval(self, values, x):
values[0] = value1
values[1] = value2
values[2] = value3
values[3] = value4
if x[0]<w/2.0:
values[0] = value5
values[1] = value6
values[2] = value7
values[3] = value8
def value_shape(self):
return (4,)
# Then use the interpolation, I can obtain the initial u0 as follow(but it can't work in the newest version, I don't know why?):
u0.interpolate(InitialConditions())
So, my question is:
1. How to define the function space for a multiple dofs problem(each node has multiple dofs), not a vector field but a multiple scalar components(For example: c1,c2,c3...,cn).
2. How to use some high order element or other different element for these c1,c2,c3,phi? For my case, I don't need to consider VectorFunctionspace or other things. Just some simple 2D-4Node rectangle element or 2D-8Nodes rectangle element for all the c1,c2,c3,phi.
I tried use higher order, but FEniCS given me some warning like "too much gauss points...", some thing like that. How to solve it?
3. How to use Functionspace to interpolate the initial value for each dof. In the old version, I can use the self-defined class (class InitialConditions(...):) and use
u0.interpolate(InitialConditions())
but now it gives me the errors.
4. To my understanding, for a multiple dofs problem,For example 4 dofs each node. Then I can define the test function as follow:
v1,v2,v3,v4=Functions(V)
du=TrialFunction(V)
u0=Function(V)
u0.interpolate(InitialCondition())
c10,c20,c30,phi0=split(u0)
u=Function(V)
c1,c2,c3,phi=split(u)
then after I defined the weak form $L$ for this problem use (v1,v2,v3,v4),(c10,c20,c30,phi0) and
(c1,c2,c3,phi). I can compute the Jacobian of this system in each time step as follow:
a=derivative(L,u,du)
problem=NonlinearVariationalProblem(L,u,bcs=bc,J=a)
solver=NonlinearVariationalSolver(problem)
solver.solve()
u0.assign(u)
Is this understanding right? Because I did some simple test and it seems every is OK. But for a complex case, it's not so easy to do the simulation.
Any helpful suggestion is welcome.
Thank you!
