I want to find the largest value of my PDE solution on the boundary of the domain. I know “u.max” gives the largest value of the u on the mesh, but I don’t know how to modify it to get the largest function value on the boundary of the domain.

Thanks for your reply! I understand this strategy but I don’t know what is the syntax for them. For example, how to identify the element that lies on the boundary?

Thank you, it works. My understanding of this code is that, when envaluates onBoundary(0,vh, tgv=-1), it returns a vector that is zero on the interior point and is 1 for all boundary points. Therefore, u=u.*onB makes the function value on the interior point to be zero, which stands out from the boundary points. Is that correct understanding? If so I think an addition modification is required if the function value is negative on the boundary. In addition to that, could I ask why tgv is set to be -1 when constructs onB?

Your understanding is correct. No modification is needed if the function value on the boundary is negative. A negative number multiplied by 1 is still negative.