Dividing by a function in C++

Question

Dear All,

I would like to define a Function in C++.
The question is: How can we divide by a function?

Here is some code to explain the issue.

First some declarations:

// Create mesh and function space
int nx = 10;
IntervalMesh mesh(nx, 0, 1);
Band::FunctionSpace V(mesh);

// Set parameter values
double g0 = 5;
double eta = .05;

// Create functions
Function u(V); 
Function f(V); 
Function g(V);  

Now initialisation

//Create intial conditions and interpolate
SourceIni u_init; //u_init is zero plus a small gaussian noise
u.interpolate(u_init); 

//initial condition for function f

//compute <u>
double umean=0;
for (MeshEntityIterator e(mesh, 0); !e.end(); ++e){
    umean+=(double)(*u.vector())[e->index()];
}
umean=umean/nx;

Now I want to create g(x) = g0 + (<u> – u(x))/eta.
I do it in this way, which works, even though it doesn't look pretty (any comments on that?):

*g.vector() = *u.vector();
*g.vector() -= umean;
*g.vector() *= -1./eta; 
*g.vector() += g0; 

Here is where I need help:

I want to define f(x) = g(x)/(1 + g(x)). Note that g represents a physical quantity that needs to be positive.

//THIS DOES NOT WORK:
*temp.vector() = *g.vector();
*temp.vector() += 1.;  
f = project(g/temp, V); 
L.f = f; 

Later on, there is a time loop, in which the function u is calculated at each time step, and g and f need to be recomputed at each time step too, as I need to update L.

Thanks for any hints!

---

PS the .ufl file is:

element = FiniteElement("Lagrange", interval, 1)
k  = Constant("interval")
D = 0.001
# Define variational problem
u = TrialFunction(element)
v = TestFunction(element)
u1 = Coefficient(element)
f = Coefficient(element)
a = ((1 + k)*inner(u,v) + D*k*inner(grad(u),grad(v)))*dx
L = inner((u1 + k*f),v)*dx 

Answer 1

g can be defined and updated, in an elegant way, using a function to calculate the operations (like in the elastodynamics demo).

For the other hand, i think that the best way to define f is do it in the .ufl file, for example:

# Define variational problem
g  = Coefficient(element)
f  = g/(1 + g)
a = ((1 + k)*inner(u,v) + D*k*inner(grad(u),grad(v)))*dx
L = inner((u1 + k*f),v)*dx 

of this manner you only need to recompute and/or update g.

Comments

1. Thank you for you advice. Indeed, the elastodynamics demo is a good suggestion as to how I can define the function g in a more elegant way.

2. As you suggested, I included into my .ufl file:

   g  = Coefficient(element)
   f  = g/(1 + g)
   

   This does not seem to work as I get an error at compilation ("error: no member named ' f ' " ) at the line:

   L.f = f
   

   Do you have any idea how to fix this?

---

Also, when you say:

you only need to recompute and/or update g.

Do you imply that the function f would automatically be updated anytime g changes?

Thank your time and any further hints!

---

The solution to the error that you get is to write

L.g = g;

instead of

L.f = f;

(remember that f is equal to g/(1 + g), i.e., f isn't an Expression that must be initialized in your .cpp code).

Answering your last question: yes, because f has been written explicitly in the .ufl code (and f depend on g).

---

Thank you for the precisions you gave.
Using L.g = g is indeed a good solution
It is now compiling and running with no errors.

Thank you for your help in closing this Question.