7
$\begingroup$

I want to interpolate a type of data which is on a triangular lattice in order to make DensityPlot faster(ListDensityPlot is so slow). However, the interoplation failed with error, no matter setting InterpolationOrder->1 or InterpolationOrder->All

Here is a less dense data

data = Import["https://pastebin.com/raw/6XmFDzmf"]; plotData = Partition[ StringCases[data, x : NumberString :> Internal`StringToDouble@x], 3]; Interpolation[plotData]

it will raise several errors

Interpolation::udeg: Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will be reduced to 1.

Interpolation::femimq: The element mesh has insufficient quality of 0.`. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements.

Interpolation::fememtlq: The quality 0.of the underlying mesh is too low. The quality needs to be larger than 0..

What does it mean? How to interpolate such data or general non-rectangular data?

PS: the density plot is

enter image description here

Improve this question
$\endgroup$

2 Answers 2

Reset to default
7
$\begingroup$

The problem is that for regular (but not rectangular) meshes the Delaunay mesh is unstable. It's a bug, also mentioned in this question, which was about rectangular grids. The same workaround works here -- just jiggle the points a tiny bit around their perfect lattice positions:

epsilon = 10^-7; jiggledPlotData = {#1 + RandomReal[epsilon {-1, 1}], #2 + RandomReal[epsilon {-1, 1}], #3} & @@@ plotData; reg = ConvexHullMesh[jiggledPlotData[[All, 1 ;; 2]]]; f = Interpolation[jiggledPlotData, InterpolationOrder -> 1]; DensityPlot[f[x, y], {x, y} \[Element] reg]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you so much! You saved my day! $\endgroup$ Commented Dec 16, 2017 at 6:00
3
$\begingroup$

A better method than jiggle original grid which I learned recently is dealing with mesh explicitly.

First, we need a package

Needs["NDSolve`FEM`"];

We can generate Delaunay mesh by

originalMesh = ToElementMesh[plotData[[;; , 1 ;; 2]]]

but this will throw a warning message

ToElementMesh::femimq: The element mesh has insufficient quality of 0.`. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements.

this is the same femimq problem, though the mesh looks fine at first sight as(enlarged)

enter image description here

The ElementMesh provides powerful tools to examing what is going on.

originalMesh["Quality"][[1]] // Histogram

enter image description here

The histogram shows most of the triangle has good quality close to 1. But there are plenty of triangle has quality near 0 which are bad enough.

using highlightBadMeshElement(code at the end), we can show the location of these bad triangles.

highlightBadMeshElement[originalMesh, 0.2]

enter image description here

They are all located at the edge shown by red marker.

We can delete these redundent bad triangles and form a refined mesh using refineMesh(code at the end)

mesh = refineMesh[originalMesh, 0.2]; f = ElementMeshInterpolation[{mesh}, plotData[[;; , -1]]]; DensityPlot[f[x, y], {x, y} \[Element] mesh]

gives

enter image description here

functions:

highlightBadMeshElement[mesh_, qualityThreshold_] := Module[{}, pos = Position[mesh["Quality"], _?(# <= qualityThreshold &)]; Show[ mesh["Wireframe"[ "ElementMeshDirective" -> Directive[EdgeForm[GrayLevel[.6]], FaceForm[]]]], mesh["Wireframe"[pos, "MeshElement" -> "MeshElements", "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]] , Boxed -> False]] refineMesh[mesh_, qualityThreshold_] := Module[{}, qualityList = First@mesh["Quality"]; pos = Flatten@Position[qualityList, _?(# > qualityThreshold &)]; If[Length@mesh["MeshElements"] != 1, Message[refineMesh::MoreThanOneTypeOfElement, Head /@ mesh["MeshElements"]];Abort[], elementHead = Head@First@mesh["MeshElements"]; ToElementMesh["Coordinates" -> mesh["Coordinates"], "MeshElements" -> {elementHead[ mesh["MeshElements"][[1, 1, pos]]]}]]]; refineMesh::MoreThanOneTypeOfElement = "More than one type of element. `1`";
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.