Computing derivatives in C++

Question

Hello All,

I am wondering how to calculate the derivative of solution in C++.

auto wp = std::make_shared<Function>(V);


// Solve for the weighting potential
auto wp_problem = std::make_shared<LinearVariationalProblem>(a, L, w, wbc);
LinearVariationalSolver solver(wp_problem);
solver.parameters["linear_solver"]="mumps";
solver.solve();

My w is fine, and it the right solution. What I need is the gradient of w. I tried this unsuccessfully:

// Now I need to calculate the gradient of w
auto wgrad = std::make_shared<Poisson::LinearForm>(V);
solve(a==w,wgrad)

my Poisson.ufl:

element = FiniteElement("Lagrange", tetrahedron, 1)
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
a = inner(nabla_grad(u),nabla_grad(v))*dx
L = f*v*dx

Can anyone help?

Thanks
Victor.

Answer 1

you would need to do it using variational forms:

element = VectorElement("Lagrange", tetrahedron, 1) # or perhaps DG
u = TrialFunction(element)
v = TestFunction(element)

a = inner(u, v)dx
L = inner(grad(solution),v)dx

Comments

Thanks Mardal.

I tried that.

Here are my Gradient.ufl files:

vector_element = VectorElement("Lagrange",tetrahedron , 1)
element = FiniteElement("Lagrange",tetrahedron , 1)

u = Coefficient(element)
w = TrialFunction(vector_element)
v = TestFunction(vector_element)


a = inner(w, v)*dx
L = inner(grad(u), v)*dx

and my code is as follow:

auto L1    = std::make_shared<Gradient::LinearForm>(V);
auto a1    = std::make_shared<Gradient::BilinearForm>(V,V);
auto w1   = std::make_shared<Function>(V);
L1->u=w ;  // w is my solution from earlier
solve( a1 == L1 , w1 );

However this generates quite a bit of errors.

error: no matching function for call to ‘solve(bool, std::shared_ptr&)’
solve( a1 == L1 , w1 );

Any suggestions?
Thanks

Answer 2

I wanted to do the same thing. Here's how it is working for me:
The jpcGrad.ufl file:

element     = FiniteElement("Lagrange", triangle, 2)
velement    = VectorElement("Lagrange", triangle, 2)

solution    = Coefficient(element)   
v           = TestFunction(velement)  # Test Function
u           = TrialFunction(velement) # Trial Function

a = inner(u, v)*dx 
L = inner(grad(solution),v)*dx

Now, in the .cpp file:

#include "jpcGrad.h"

...

class GradClass : public Expression
{
public:
    GradClass(Function &cref)
    {
        c = &cref;
    }
    void eval(Array<double>& values, const Array<double>& x) const
    {
        calledGradEval();
        Array<double> myval(1);
        c->eval(myval, x);
        values[0] = myval[0];
    }
private:
    Function *c;
};

...

jpcGrad::FunctionSpace G(mesh, PBC);
jpcGrad::LinearForm L(G);
jpcGrad::BilinearForm a(G,G);
auto gradf = GradClass(c);
L.solution = gradf;
Function u(G);

...//now in  the time loop:

//this is the solver that gives the solution c
Solver.solve();
//now solve for grad c:    
solve(a == L, u);

Thus, I had to subclass Expression so that the solution (a reference to which is saved in the private data of the subclass) can be evaluated in the gradient variational form assembly.