A better way to integrate the gradient along a boundary?

Question

I'm currently integrating the temperature gradient along a boundary by first solving (seems wrong to need a solver) for the gradient of the temperature field and then integrating with Simpson's rule. Is there a cleaner approach with the project I am missing?

I am working in c++, so I have created a ufl file just for calculating the gradient. I've seen some python snippets that do this cleanly with assembling a matrix with 'grad', but I haven't found the c++ equivalent. (I assume its a ufl -> form -> assemble?)

A general way to integrate the gradient over marked facets would be really convenient :)

Example:

UFL file for gradient

vector = VectorElement("Lagrange", triangle, 2)
scalar = FiniteElement("Lagrange", triangle, 2)

q = Coefficient(scalar)  

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

Then I would do a simple solve for grad_q

solve(a_grad == L_grad, grad_q);

With grad_q I would integrate along an edge, in this case the negative x-normal from (1,0) to y=(1,1)

double Qn = 0;
int points = 500;
for(int j = 0; j <= points - 1; j++)
{
    double y0 = j / (double)points;
    double y1 = (j+1) / (double)points;
    double yab = 0.5*(y0+y1);

    Point pa(x,y0);
    Point pb(x,y1);
    Point pab(x, yab);
    Array<double> values(2);

    grad_q(values, pa);
    fa = -values[0]; // x-component

    grad_q(values, pab);
    fab = -values[0]; // x-component

    (grad_q(values, pb);
    fb = -values[0]; // x-component

    Qn += ((y1-y0)/6.0)*(fa + 4*fab+fb);
}

This works, but only because the path is well defined and the normal is also well defined. I'm thinking of an unstructured/curved boundary where I would need to integrate along the facets point to point and calculate the normal components as being painful.

If there was a way to integrate a Function (not necessarily gradient) along a marked boundary that would help, unless I am missing a more direct method?

Comments

Please post simple code that illustrates what you're trying to do.

---

Added, thank you!

Answer 1

Hi, have a look at this demo. There's also a cpp version.

Comments

Excellent, thank you. And just to confirm, is there a better approach to finding the derivative? I would think I could just apply an operator without the need for a solver?

---

I can't think of any better(easier) approach. As far as applying grad goes I think there might be some issues. Suppose you have $ f(x) = \sum_{i} F_i \phi_i(x) $ and then you'd apply grad to get $ grad(f)(x) = \sum_{i} F_i grad(\phi_i)(x) $. Now, this vector can be evaluated in the cell interior just fine, but on the boundary it's discontinuous. You could try weighting contributions from connected cells (see here) to get a unique value, but overall project seems more convenient.

---

Very interesting, and why i've been focusing on the gradient so much (that singularity is making converging on the integral of a boundary's gradient rather painful).

How might I do this in c++?

dp_regular = project(p.dx(0), V)

---

Some digging, and actually finding your previous answer on projection, it seems to not be possible in c++.

---

Hi, I am a bit lost. Are you trying to do in cpp what project does in python?

---

Apparently I am a bit lost as well ;)

Your results hinted that I might solve an issue with the gradient in a corner of my unit cell by using project, so I was trying to find out how I might do that without doing it through ffc/solver. I don't think there is an equivalent in c++ of projecting an expression with UFL terms onto a Function, just interpolating of Expressions with c++ or other Functions. (Is there a way to compile UFL expressions to c++ that can be interpolated? That seems to be whats missing)

I started digging into your weighted gradient matrix to see how I might do that in c++, but I'm not sure how I would assemble this since its just a UFL expression:

dP = assemble(TrialFunction(S).dx(i)*TestFunction(DG)*dx)

Looking into how the weighting works, I think solving for the gradient on a DG space will help me. (So far testing seems okay)