Calculate values from function

Linear Algebra
1

Hi,

I’m solving an advection-diffusion-equation in a mixed system from which I obtain the concentrations for three substances, c1, c2 and c3. From these concentrations I want to calculate saturation functions with the following formula:

Capture

Where Ct is a threshold, so a constant, and phi is the saturation function, and C is the concentration c1,c2 or c3. I tried this (for c3 in this example):

element = FiniteElement('P', triangle, 1)  
Vphi = FunctionSpace(mesh, element)
phi = Function(Vphi)
ptest = TestFunction(Vphi)

thr = Constant(0.002)
form = ((c3**2)/((c3**2)+(thr**2)))*ptest*dx

phi = assemble(form)

This works, but in this case phi is of the type ‘dolfin.cpp.la.Vector’ which I want to transfer back to a Function. The other option I tried is:

thr = 0.002
c3_vec = c3.vector().get_local()
phi_vec = (c3**2)/((c3**2)+(thr**2))
phi.vector().set_local(phi_vec)

This does work normally, but since c3 comes from a mixed functionspace I can’t convert it to c3_vec. The error I get is: ‘Unable to access vector of degrees of freedom. cannot access a non-const vector from a subfunction.’

So my question is: 1. Which of the two methods is correct, and 2. How do I solve the described problems?

Thanks in advance

2

Does phi = project( ((c3**2)/((c3**2)+(thr**2))) , Vphi) provide a good approximation in your application?

(If there are sharp layers in the quantity to be projected, the L^2 projection onto CG1 may have some oscillations, but you could instead try a lumped-mass projection, as demonstrated here, which also avoids doing a linear solve.)

3

Hi,

of course, you could also work with Expressions, and then you can interpolate them into an arbitrary FunctionSpace (i.e. one that fits your needs). This could look like this:

expr = Expression('pow(c, 2) / (pow(c, 2) + pow(thr, 2))', degree=k, c=c3, thr=thr)

In regards to the first idea, the projection is exactly what you want to do since this is equivalent to solve the following variational problem (with “form” defined as in your example, and I write down the discretized version since this is easier to do)

u = TrialFunction(Vphi)
a = u*ptest*dx
phi = Function(Vphi)
solve(a == form, phi)

Note that you may or may not want to impose additional Dirichlet boundary conditions here.

For your second idea: Have you tried to extract the solution from the mixed space by using

c3.sub(i, True)

where i is the index of the component of the mixed space?

Personally, I would prefer using an Expression and interpolating it, since it should give you the nodally exact values for Lagrangian finite elements.

1 Like