Extrapolation value for interpolate()

Asked by Peter Maday

I have a 3D scalar field as an output during a solution process defined over an unstructured mesh. I would need to resample it to a regular mesh so that it can be further processed with some existing tools that operate on voxel data. My question is if the interpolate() function is the recommended way to do this, or are there better alternatives?
Using the interpolate() in case the domain of the two grids is not identical certain values will be extrapolated. Is there a way to control this process i.e. to set the extrapolated values to zero instead of obtaining them via some default extrapolation technique? If not is there a simple& fast way to identify cells that are not contained within the other mesh so that manual corrections can be done?

Thanks a lot,
Peter

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Anders Logg
Solved:
Last query:
Last reply:
Revision history for this message
Garth Wells (garth-wells) said :
#1

On Sat, Dec 29, 2012 at 11:41 AM, Peter Maday
<email address hidden> wrote:
> New question #217871 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/217871
>
> I have a 3D scalar field as an output during a solution process defined over an unstructured mesh. I would need to resample it to a regular mesh so that it can be further processed with some existing tools that operate on voxel data. My question is if the interpolate() function is the recommended way to do this, or are there better alternatives?

Yes, if the domain of your second mesh is contained in the the domain
of the first mesh.

> Using the interpolate() in case the domain of the two grids is not identical certain values will be extrapolated.

Do you mean extrapolated beyond the domain boundaries?

> Is there a way to control this process i.e. to set the extrapolated values to zero instead of obtaining them via some default extrapolation technique? If not is there a simple& fast way to identify cells that are not contained within the other mesh so that manual corrections can be done?
>

There are some tools in 'intersection' that provide support for
detecting intersection. Take a look at the interfaces in
dolfin/intersection, and Mesh::intersected_cells

Garth

> Thanks a lot,
> Peter
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.

Revision history for this message
Peter Maday (madapeti) said :
#2

Thanks a lot for your reply,

I believe it would be best if I added some more details about the goal I would like to achieve.
I have an unstructured (tetrahedral) mesh for the interior of some blood vessels as the domain on which the solution is available. I would like to resample the solution to a uniform (voxel) grid, where the domain would be the interior of the bounding box of the blood vessels (to generate projection images of a transported quantity with an existing tool). For this reason the majority of the cells (voxels) would not intersect the original domain.
What I would thus like to achieve is to have a predefined value (0) at these points instead of some value depending on the values of the closest points of the original mesh. (analogous to the "extrapval" value of the interpolation functions in MATLAB).

Considering the intersection tools, I am a little afraid that in case you have mesh sizes of around 100k cells doing pairwise intersection tests would make the conversion extremely slow. I will take a look at the features of the CGAL library to see if there would be some features that could be useful. (especially intersecting a mesh with a plane and an accelerated range based query (KD-tree etc.)).

Thanks a lot for your help,
Peter

Revision history for this message
Best Anders Logg (logg) said :
#3

On Wed, Jan 02, 2013 at 11:16:00AM -0000, Peter Maday wrote:
> Question #217871 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/217871
>
> Status: Answered => Open
>
> Peter Maday is still having a problem:
> Thanks a lot for your reply,
>
> I believe it would be best if I added some more details about the goal I would like to achieve.
> I have an unstructured (tetrahedral) mesh for the interior of some blood vessels as the domain on which the solution is available. I would like to resample the solution to a uniform (voxel) grid, where the domain would be the interior of the bounding box of the blood vessels (to generate projection images of a transported quantity with an existing tool). For this reason the majority of the cells (voxels) would not intersect the original domain.
> What I would thus like to achieve is to have a predefined value (0) at these points instead of some value depending on the values of the closest points of the original mesh. (analogous to the "extrapval" value of the interpolation functions in MATLAB).
>
> Considering the intersection tools, I am a little afraid that in case
> you have mesh sizes of around 100k cells doing pairwise intersection
> tests would make the conversion extremely slow. I will take a look at
> the features of the CGAL library to see if there would be some features
> that could be useful. (especially intersecting a mesh with a plane and
> an accelerated range based query (KD-tree etc.)).
>
> Thanks a lot for your help,
> Peter

The intersection test is quite fast so I wouldn't worry too much, at
least if you do this from C++.

To speed things up, use a bounding box for the small mesh and iterate
only over those points. Then for each point, check whether it
intersects the small domain and if so interpolate.

--
Anders

Revision history for this message
Peter Maday (madapeti) said :
#4

Thanks Anders Logg, that solved my question.