Thread Subject: data gridding on a 2-d mesh with voids

Hi all,

I'm wrestling with something but simply unfamiliar enough with either the background math or matlab itself to figure it out.

I have a dataset of 2D triangular mesh data from an external solver (actually a slice of an asymmetric 3D FEM simulation), which I'm attempting to plot. I don't have a problem griddding and plotting it with say the method below:

tx = -100:1:100;
ty = -100:1:100;
[XI,YI] = meshgrid(tx,ty);
ZI = griddata(x,y,Efield,XI,YI,'cubic');

What concerns me is that the original mesh itself contains a number of voids, empty regions with no points which define certain objects with fixed boundary conditions. When I grid, of course, it also goes into those regions, so the interpolated values at the boundary are of necessity slightly off (the interior is irrelevant, but the most important information for me lies right near the boundary of the objects, and if that is artificially reduced by averaging with the inside (where its arbitrarily zero), then its a problem.

Suggestions, ideas where to start? From what I understand, griddata can't handle convex regions, so if I set up a 2d grid excluding them then it will surely choke. I imagine the problem must be common but I haven't had luck finding the solution as yet.

Cheers!

Shawn wrote:
...
> I have a dataset of 2D triangular mesh data from an external solver
> (actually a slice of an asymmetric 3D FEM simulation), which I'm
> attempting to plot. ...
> What concerns me is that the original mesh itself contains a number of
> voids, empty regions with no points which define certain objects with
> fixed boundary conditions. When I grid, of course, it also goes into
> those regions, so the interpolated values at the boundary are of
> necessity slightly off (the interior is irrelevant, but the most
> important information for me lies right near the boundary of the
> objects, and if that is artificially reduced by averaging with the
> inside (where its arbitrarily zero), then its a problem.
> Suggestions, ideas where to start? From what I understand, griddata
> can't handle convex regions, so if I set up a 2d grid excluding them
> then it will surely choke. I imagine the problem must be common but I
> haven't had luck finding the solution as yet.

AFAIK ML doesn't deal w/ the problem of arbitrary meshes well at all.
My solution for such problems has always been Tecplot from Tecplot, Inc.
(formerly Amtec Engineering) -- <www.tecplot.com>

Not the c.s-s.m way, but is a tool specifically by FEM'ers for FEM
modelers...

--

dpb wrote:
> Shawn wrote:
> ...
>> I have a dataset of 2D triangular mesh data from an external solver
>> (actually a slice of an asymmetric 3D FEM simulation), which I'm
>> attempting to plot. ...
>> What concerns me is that the original mesh itself contains a number of
>> voids, empty regions with no points which define certain objects with
>> fixed boundary conditions. ...

Something like this, maybe???

<http://www.tecplot.com/showcase/gallery/mesh/18.asp>

--

In article <hgvv10$pos$1@fred.mathworks.com>, sfostner@physics.mcgill.ca
says...
> Hi all,
>
> I'm wrestling with something but simply unfamiliar enough with either the background math or matlab itself to figure it out.
>
> I have a dataset of 2D triangular mesh data from an external solver (actually a slice of an asymmetric 3D FEM simulation), which I'm attempting to plot. I don't have a problem griddding and plotting it with say the method below:
>
> tx = -100:1:100;
> ty = -100:1:100;
> [XI,YI] = meshgrid(tx,ty);
> ZI = griddata(x,y,Efield,XI,YI,'cubic');
>
> What concerns me is that the original mesh itself contains a number of voids, empty regions with no points which define certain objects with fixed boundary conditions. When I grid, of course, it also goes into those regions, so the interpolated values at the boundary are of necessity slightly off (the interior is irrelevant, but the most important information for me lies right near the boundary of the objects, and if that is artificially reduced by averaging with the
inside (where its arbitrarily zero), then its a problem.
>
> Suggestions, ideas where to start? From what I understand, griddata can't handle convex regions, so if I set up a 2d grid excluding them then it will surely choke. I imagine the problem must be common but I haven't had luck finding the solution as yet.
>
> Cheers!

Consider some of these MATLAB tools: trimesh, trisurf
Or perhaps some of the delaunay triangulation tools can help set up the
data in the right form first.

--
Loren
http://blogs.mathworks.com/loren

"Shawn " <sfostner@physics.mcgill.ca> wrote in message <hgvv10$pos$1@fred.mathworks.com>...
> Hi all,
>
> I'm wrestling with something but simply unfamiliar enough with either the background math or matlab itself to figure it out.
>
> I have a dataset of 2D triangular mesh data from an external solver (actually a slice of an asymmetric 3D FEM simulation), which I'm attempting to plot. I don't have a problem griddding and plotting it with say the method below:
>
> tx = -100:1:100;
> ty = -100:1:100;
> [XI,YI] = meshgrid(tx,ty);
> ZI = griddata(x,y,Efield,XI,YI,'cubic');
>
> What concerns me is that the original mesh itself contains a number of voids, empty regions with no points which define certain objects with fixed boundary conditions. When I grid, of course, it also goes into those regions, so the interpolated values at the boundary are of necessity slightly off (the interior is irrelevant, but the most important information for me lies right near the boundary of the objects, and if that is artificially reduced by averaging with the inside (where its arbitrarily zero), then its a problem.
>
> Suggestions, ideas where to start? From what I understand, griddata can't handle convex regions, so if I set up a 2d grid excluding them then it will surely choke. I imagine the problem must be common but I haven't had luck finding the solution as yet.
>
> Cheers!

Griddata will span the interior holes with triangles,
and the triangles must of course be large ones.

If you already have the triangulation at hand, then
the simple way to plot it is to use trisurf. trisurf will
work even if the region is not convex. The holes
will be no problem, since you have the triangles.

If you don't have a triangulation, then the simplest
thing is to use gridfit to fit the entire region. It
too will span the interior holes, but it will interpolate
the edges of the holes far more smoothly than
griddata.

I do have a fully general (and adaptive) code to
model non-convex regions, but I have not yet
gotten it ready for general dispersal. Gridfit should
be sufficient here.

http://www.mathworks.com/matlabcentral/fileexchange/8998

HTH,
John
John

"Shawn " <sfostner@physics.mcgill.ca> wrote in message <hgvv10$pos$1@fred.mathworks.com>...
> Hi all,
>
> I'm wrestling with something but simply unfamiliar enough with either the background math or matlab itself to figure it out.
>
> I have a dataset of 2D triangular mesh data from an external solver (actually a slice of an asymmetric 3D FEM simulation), which I'm attempting to plot. I don't have a problem griddding and plotting it with say the method below:
>
> tx = -100:1:100;
> ty = -100:1:100;
> [XI,YI] = meshgrid(tx,ty);
> ZI = griddata(x,y,Efield,XI,YI,'cubic');
>
> What concerns me is that the original mesh itself contains a number of voids, empty regions with no points which define certain objects with fixed boundary conditions. When I grid, of course, it also goes into those regions, so the interpolated values at the boundary are of necessity slightly off (the interior is irrelevant, but the most important information for me lies right near the boundary of the objects, and if that is artificially reduced by averaging with the inside (where its arbitrarily zero), then its a problem.
>
> Suggestions, ideas where to start? From what I understand, griddata can't handle convex regions, so if I set up a 2d grid excluding them then it will surely choke. I imagine the problem must be common but I haven't had luck finding the solution as yet.
>
> Cheers!

The TriScatteredInterp tool (R2009a and later) was specifically designed to handle this usecase. You first need to create a constrained Delaunay triangulation (DelaunayTri), where the constraints represent the boundary that you wish to respect. You then use the constrained DelaunayTri to create a TriScatteredInterp (interpolant).
This example desmonstrates how to create a constrained Delaunay triangulation.
http://www.mathworks.com/products/demos/shipping/matlab/demoDelaunayTri.html#18

If you haven't used TriScatteredInterp before, here's an overview on Loren's Art of MATLAB blog.
http://blogs.mathworks.com/loren/2009/07/22/computational-geometry-in-matlab-r2009a-part-ii/

When you create the constrained Delaunay you could, if you wish, replicate your FE mesh by constraining all edges and not just the edges on the free boundary. Either way, you should create a TriRep (Triangulation Representation) to represent your FE mesh as it can provide you with the Free Boundary edges or all edges which you will need to constrain your triangulation.

Damian

Thanks for the great replies all.

I think the constrained delaunay mesh will work best, I actually didn't realize that I didn't have that, as I've been using 2008b up until now. As I have the option, I upgraded to 2009b and am playing with it now.

Tecplot loots like a good tool, but not one I have access to unfortunately.

Shawn wrote:
> Thanks for the great replies all.
> I think the constrained delaunay mesh will work best, I actually didn't
> realize that I didn't have that, as I've been using 2008b up until now.
> As I have the option, I upgraded to 2009b and am playing with it now.
> Tecplot loots like a good tool, but not one I have access to unfortunately.

I'm stuck w/ a much earlier release of ML -- glad to see they apparently
have finally gotten more flexibility for mesh plotting than this release
and earlier. I'll not be giving up Tecplot, though, methinks... :)

--

Ok, so I have tried a few things, but not sure yet what the best solution is and some input would be nice:

1) I can take the original point cloud and similar to outlined in the original example (7), I can do the triangulation based on that. (in this case its a 2d slice of what was a 3d matrix, but it still works for the most part). I'm still having trouble with the core of the procedure though:
% The Delaunay edges that connect pairs of sample points represent the
% boundary.
delEdges = dt.edges();
validx = delEdges(:,1) <= numpts;
validy = delEdges(:,2) <= numpts;
boundaryEdges = delEdges((validx & validy), :)';
xb = x(boundaryEdges);
yb = y(boundaryEdges);

Not entirely sure what the point of this is, if I read it correctly its really only finding those points with indices < the number of original points on the curve (not the triangular mesh which has much more edges). Does it actually care about the BOUNDARY, or just the original polygon?

I've gotten that far with my data, but from this triangulated cloud I don't know where to go. Now I have a 2d delaunay mesh, but an unconstrained one at that, still. I think.

2) I can also take the ORIGINAL data file, ie a vtk mesh with both 3d tetras and 2d face triangles, and import it into matlab. This gives me 3 matrices, points, cells (however some have 3 vertices specified only, where they are part of the surface mesh), and the actual POINT DATA, eg the relevant quantity specified on each vertex.

However, I'm not so sure about getting that into a delaunaytri structure, I tried creating a dummy structure to do it, then inserting them, interesting enough it removed all duplicate points for me automatically. However, that means that my point data matrix is no longer quite correct, as some points/cells were removed, and I don't know how to deal with that. Well that and the fact that it crashed matlab badly when I tried to plot it, rather than just the boundary elements. :) My bad.

Thoughts, suggestions?

Cheers
Shawn

Hello Damian
I have a question on constrained Delaunay:
There are say 600 vetices and the boundary is only consisted of 300 points, and it comes up this error message saying "??? Error using ==> DelaunayTri
The column-vectors defining the edges have inconsistent dimensions."

Could you give some advice?

Thank you!
-Hualei
"Damian " <damian.sheehy@mathworks.com> wrote in message <hhacn4$ddq$1@fred.mathworks.com>...
> "Shawn " <sfostner@physics.mcgill.ca> wrote in message <hgvv10$pos$1@fred.mathworks.com>...
> > Hi all,
> >
> > I'm wrestling with something but simply unfamiliar enough with either the background math or matlab itself to figure it out.
> >
> > I have a dataset of 2D triangular mesh data from an external solver (actually a slice of an asymmetric 3D FEM simulation), which I'm attempting to plot. I don't have a problem griddding and plotting it with say the method below:
> >
> > tx = -100:1:100;
> > ty = -100:1:100;
> > [XI,YI] = meshgrid(tx,ty);
> > ZI = griddata(x,y,Efield,XI,YI,'cubic');
> >
> > What concerns me is that the original mesh itself contains a number of voids, empty regions with no points which define certain objects with fixed boundary conditions. When I grid, of course, it also goes into those regions, so the interpolated values at the boundary are of necessity slightly off (the interior is irrelevant, but the most important information for me lies right near the boundary of the objects, and if that is artificially reduced by averaging with the inside (where its arbitrarily zero), then its a problem.
> >
> > Suggestions, ideas where to start? From what I understand, griddata can't handle convex regions, so if I set up a 2d grid excluding them then it will surely choke. I imagine the problem must be common but I haven't had luck finding the solution as yet.
> >
> > Cheers!
>
> The TriScatteredInterp tool (R2009a and later) was specifically designed to handle this usecase. You first need to create a constrained Delaunay triangulation (DelaunayTri), where the constraints represent the boundary that you wish to respect. You then use the constrained DelaunayTri to create a TriScatteredInterp (interpolant).
> This example desmonstrates how to create a constrained Delaunay triangulation.
> http://www.mathworks.com/products/demos/shipping/matlab/demoDelaunayTri.html#18
>
> If you haven't used TriScatteredInterp before, here's an overview on Loren's Art of MATLAB blog.
> http://blogs.mathworks.com/loren/2009/07/22/computational-geometry-in-matlab-r2009a-part-ii/
>
> When you create the constrained Delaunay you could, if you wish, replicate your FE mesh by constraining all edges and not just the edges on the free boundary. Either way, you should create a TriRep (Triangulation Representation) to represent your FE mesh as it can provide you with the Free Boundary edges or all edges which you will need to constrain your triangulation.
>
> Damian