Hi everyone, lately I’ve been using Julia to solve a stochastic partial differential equation for a 2D problem. I used a refined triangular mesh for my mesh generation, but I encountered some issues while implementing my log-likelihood for the Bayesian estimation of the range correlation parameter and the standard deviation. Specifically, I’m having trouble performing the integral over the triangular mesh. If anyone can help me with that, I would greatly appreciate it. Thank you!
You will probably get faster help if you include the mathematical details / the code you have thus far
Thank you for your reply,
So basically i used the function Trianglemesh.jl to generate a refined triangulation over a space, now i just need to estimate the numerical integral of my values at each points (red circle) in this mesh like presented here
If you have values only at the red dots, then you’ll need to decide how you want to interpolate them in order to integrate.
One option would be to create a triangulation using only the red dots (no refinement), and then do bilinear interpolation on each triangle. This corresponds to multiplying the area of each triangle by the average of the values at its three vertices.
(I don’t see how your “refined” triangulation will be relevant, however, unless you define some other interpolation scheme based on the refined mesh.)
Another option would be to interpolate using some other method, e.g. radial basis functions, and then to integrate the resulting smooth interpolant with a generic quadrature scheme.
It’s also not clear to me over what domain you want to integrate. If you want to integrate over the whole rectangular domain that you showed, you’ll have to decide on some extrapolation scheme for how to define the value of your integrand outside of the convex hull of your red dots.
3 Likes
I’ll note that another good smooth interpolant would be most of the interpolants from NaturalNeighbours.jl (e.g. Farin(), just not Nearest() or probably Triangle()) that are defined directly using triangulations (technically it’s actually based on Voronoi tessellations)
Indeed though you might need to be a bit clearer on what you’re interpolating over / your domain if @stevengj’s approach is not clear.
1 Like
Would the package MeshIntegrals.jl be of use here?
1 Like
Thanks for the suggestion 
It’s not clear to me exactly what @sofiane is trying to accomplish here based on the description and image, i.e. is this:
- calculate the value of a surface integral \int_0^5\int_0^5 f(x,y) \text{d}x\,\text{d}y where only the red points have known values?
- at each red point, formulate some integral problem and solve it?
In this first case, it sounds like most of the work will be in figuring out how to calculate/estimate values for other points in the domain, whether by some kind of interpolation scheme, direct calculation, etc. Some kind of post-hoc Monte Carlo integration scheme might also be possible, but results might not be particularly accurate, especially given how non-uniform the sampling distribution appears to be.
In any event, if you have an integral problem where the domain can be defined via a Meshes.jl geometry, and you have a field-like integrand that can be defined as a function of position within that geometry, then MeshIntegrals.jl might be helpful.
2 Likes
Are you implementing the SPDE method to approximate a Matern process (Lindgren et al. 2011)? If so, one strategy I’ve used is to find the Voronoi polygons of the mesh and use these as weights for the integration. This is also the strategy recommended in Simpson et al. (2012), where the motivation for the integral is to estimate a log-Gaussian Cox process.
On a side note, have you seen SPDEPrecisionMatrices.jl?
For others’ context (in case there is a better way to do this integration!), this model uses the mesh with linear elements to sparsely approximate the precision matrix of a Matern Gaussian random field. The value at each vertex of the mesh is a parameter, so the values for the entire domain will be available for the integral during estimation.
Thank you for your Replay and yes it’s exactly what i’am trying to do. Thank you for the paper but Can you please maybe send a piece of code that Can Do that. If you Can of course
Unfortunately I don’t have an example. It looks TriangleMesh.jl can create a Voronoi Diagram from your mesh. The trick then is to get the area of each of the Voronoi polygons and use these as the integration weights.
GeoStats.jl already provides interpolation/simulation methods on top of Meshes.jl including the LindgrenProcess based on SPDEs.
Unlike the alternatives it works with all kinds of meshes, not just triangle meshes.