The unsteady advective nonlinear (mixed) Cahn-Hilliard system after discretization in space and time (theta scheme) reads:
CH_LHS u CH_RHS u0
| M + theta*dt*C theta*dt*m(c)*D | |c ^{n}| | M - (1-theta)*dt*C -(1-theta) dt m(c^{n-1}) D | |c ^{n-1}|
| | * | | = | | * | |
| -eps^2*D - N(c) M | |mu^{n}| | 0 0 | |mu^{n-1}|
C := Convection matrix
N := Nonlinear contribution (gets assembled into a vector)
M := Mass matrix
m(c) := mobility function
theta := Theta time stepping scheme
The residual of the system is computed as: F = CH_LHS * u - CH_RHS * u0
Due to the splitting of the fourth order in space Cahn-Hilliard system and application of integration by parts, the
highest remaining derivative in the weak forms is 1, so one can use linear finite elements as it is in fact done in the demo.
P1 = FiniteElement("Lagrange", cell, 1)
ME = P1*P1
Trying to figure out what exactly fenics assembles, I recompiled the ufl file with the option to use numerical quadrature instead of the tensor representation and didn't use any optimization flags:
$ ffc -r quadrature -l dolfin -f split CahnHilliard2D.ufl
In the generated file CahnHilliard2D.cpp the function in charge with assembling F over a cell is:
void cahnhilliard2d_cell_integral_0_otherwise::tabulate_tensor(double* A,
const double * const * w,
const double* vertex_coordinates,
int cell_orientation) const
Unfortunately the editing system has a 8000 character limitation so I could't list it in this post.
However, the code can be easily generated with ffc as I showed above.
As for the generated code, I would like to know why 6 quadrature points are used and why we have 6 basis function values whilst only linear Lagrange finite elements are used.
Moreover, it would be great to get an explanation for the totally unreadable number crunching block providing values for A[j].
While the code generated for the Poisson demo is understandble, I wonder how the generated code maps to the above depicted discrete Cahn-Hilliard system.
Any help is appreciated.
Kind regards,
Babak S. Hosseini
