Here is a trapezoid within which I am solving a simple Laplace equation with the Dirichlet boundary conditions set at its top and bottom:

coord = {{0, 0}, {1/3, 1}, {2/3, 1}, {1, 0}};
incidents = {{1, 2}, {2, 3}, {3, 4}, {4, 1}};
bm = ToBoundaryMesh["Coordinates" -> coord2, 
   "BoundaryElements" -> {LineElement[incidents]}];
mesh = ToElementMesh[bm];
sigma = 1;
nds = NDSolveValue[{Inactive[Div][
     sigma*Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0,
   DirichletCondition[u[x, y] == 0, y == 0], 
   DirichletCondition[u[x, y] == 1, y == 1]}, 
  u[x, y], {x, y} ∈ mesh]

enter image description here

Now I would like to calculate the integral of the square of the result gradient: Evaluate[Grad[nds, {x, y}].Grad[nds, {x, y}]] over x lying within the mesh, and study the dependence of the result on y.

How can I do that?

The result should probably be a nested list of with the sublists structure {yi, resi}, where yiis the point coordinate along y, and resi is the integral value in this point.

What I do understand myself, is that I could divide the initial trapezoid into several parts, and then integrate within each such a part and, thus, make a list. I am curious, if there is a more regular approach?


1 Answer 1


With two minor changes to the code

bm = ToBoundaryMesh["Coordinates" -> coord, 
     "BoundaryElements" -> {LineElement[incidents]}];
nds = NDSolveValue[{Inactive[Div][sigma*Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0, 
        DirichletCondition[u[x, y] == 0, y == 0], 
        DirichletCondition[u[x, y] == 1, y == 1]}, u, {x, y} ∈ mesh];

the first to correct a typo, and the second for convenience, we obtain

Plot3D[nds[x, y], {x, y} ∈ mesh]

enter image description here


Plot3D[Evaluate[Grad[nds[x, y], {x, y}].Grad[nds[x, y], {x, y}]], {x, y} ∈ 
    mesh, PlotRange -> All]

enter image description here

To obtain the integral described in the question, define

f[y0_] := Quiet@NIntegrate[Evaluate[Grad[nds[x, y], {x, y}].Grad[nds[x, y], {x, y}]] 
    /. y -> y0, {x, y0/3, 1 - y0/3}]

which gives, for instance,

(* 0.672611 *)
  • $\begingroup$ I see the typo you corrected (coord instead of coord2), but cannot find, what was the second correction. Please kindly let me know. $\endgroup$ Aug 5, 2015 at 8:18
  • $\begingroup$ I changed u[x. y], {x, y} ∈ mesh to u, {x, y} ∈ mesh at the end of NDSolveValue. $\endgroup$
    – bbgodfrey
    Aug 5, 2015 at 13:28
  • $\begingroup$ Thank you, I appreciate your above answer very much. $\endgroup$ Aug 5, 2015 at 13:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.