ListPlot3D de données sur un réseau "régulier" contenu dans une région "irrégulière"

Question

J'ai le problème suivant lorsque j'essaie de tracer des données 3D dans Mathematica. Les données sont d'abord calculés sur la base d'un réseau régulier, qui est je calculer

data = Flatten[Table[{x,y,f[x,y]},{x,x0,x1,dx},{y,y0,y1,dy}],1] 

Le problème est que f prend des valeurs réelles que dans un sous-ensemble non convexe U sur le plan. U est en fait assez méchant: une région "triangulaire" où tous les coins sont cuspides (imaginez la région entre trois cercles identiques avec deux d'entre eux tangents l'un à l'autre). Lorsque j'essaie de tracer data avec ListPlot3D, ce dernier trace une surface sur la coque convexe de U.

Je me demandais s'il y avait un moyen de dire à Mathematica de se limiter uniquement à U. Je pensais que puisque mes points reposent déjà sur un treillis "régulier" cela ne devrait pas être trop difficile mais je n'ai pas encore trouvé de solution.

Je suis conscient de l'option RegionFunction mais il est très coûteux de calculer dans mon cas, donc j'essaie de trouver une façon qui utilise uniquement les points déjà calculés en data.

Answer 1

Pour générer une très belle image, sans utiliser des zillions de points, vous pouvez utiliser une distribution non uniforme des points qui inclut la limite de la région que vous voulez. Voici un exemple un peu comme ce que vous décrivez. Nous commençons avec trois cercles mutuellement tangents.

circPic = Graphics[{Circle[{0, Sqrt[3]}, 1], 
    Circle[{-1, 0}, 1], Circle[{1, 0}, 1]}] 

nous écrivons une fonction booléenne qui détermine si un point dans le rectangle {-1/2,1/2} par {0, Sqrt [3]/2} se situe en dehors de toute les cercles et l'utiliser pour générer des points dans la région d'intérêt.

inRegionQ[p:{x_, y_}] := Norm[p - {1, 0}] > 1 && 
    Norm[p + {1, 0}] > 1 && Norm[p - {0, Sqrt[3]}] > 1; 
rectPoints = N[Flatten[Table[{x, y}, 
    {x, -1/2, 1/2, 0.02}, {y, 0.05, Sqrt[3]/2, 0.02}], 1]]; 
regionPoints = Select[rectPoints, inRegionQ]; 

Maintenant, nous générons la frontière. Le paramètre n détermine combien de points nous plaçons sur la frontière.

n = 120; 
boundary = N[Join[ 
    Table[{1 - Cos[t], Sin[t]}, {t, Pi/n, Pi/3, Pi/n}], 
    Table[{Cos[t], Sqrt[3] - Sin[t]}, {t, Pi/3 + Pi/n, 2 Pi/3, Pi/n}], 
    Table[{Cos[t] - 1, Sin[t]}, {t, Pi/3 - Pi/n, 0, -Pi/n}]]]; 
points = Join[boundary, regionPoints]; 

Jetons un coup d'œil.

Show[circPic, Graphics[Point[points]], 
    PlotRange -> {{-3/4, 3/4}, {-0.3, 1.3}}] 

Maintenant, nous définissons une fonction et d'utiliser ListPlot3D pour tenter de le tracer.

f[x_, y_] := -(1 - Norm[{x - 1, y}]) (1 - Norm[{x + 1, y}])* 
    (1 - Norm[{x, y - Sqrt[3]}]); 
points3D = {#[[1]], #[[2]], f[#[[1]], #[[2]]]} & /@ points; 
pic = ListPlot3D[points3D, Mesh -> All] 

D'une certaine façon, nous devons supprimer ce genre de choses qui se trouve en dehors de la région. Dans cet exemple particulier, nous pouvons utiliser le fait que la fonction est nulle sur la frontière.

DeleteCases[Normal[pic], Polygon[{ 
    {x1_, y1_, z1_?(Abs[#] < 1/10.0^6 &)}, 
    {x2_, y2_, z2_?(Abs[#] < 1/10.0^6 &)}, 
    {x3_, y3_, z3_?(Abs[#] < 1/10.0^6 &)}}, ___], 
    Infinity] 

Très bon, mais il y a quelques problèmes à proximité des points de rebroussement et il est certainement pas très général car il a utilisé une propriété spécifique de la fonction. Si vous examinez la structure de pic, vous trouverez qu'il contient un GraphicsComplex et les premiers n points dans le premier argument de ce GraphicsComplex est exactement la limite. En voici la preuve:

Most /@ pic[[1, 1, 1 ;; n]] == boundary 

Maintenant, la frontière se décline en trois composantes et nous voulons supprimer un triangle formé par les points choisis d'un seul de ces composants. Le code suivant le fait. Notez que boundaryComponents contient les indices des points qui forment la frontière et someSubsetQ [A, Bs] renvoie vrai si A est un sous-ensemble de l'un des Bs. Nous voulons supprimer les indices triangulaires dans le multi-polygone qui sont des sous-ensembles de l'un des boundaryComponents. Cela est accompli dans le code suivant par la commande DeleteCases.

Oh, et ajoutons de la décoration aussi.

subsetQ[A_, B_] := Complement[A, B] == {}; 
someSubsetQ[A_, Bs_] := Or @@ Map[subsetQ[A, #] &, Bs]; 
boundaryComponents = Partition[Prepend[Range[n], n], 1 + n/3, n/3]; 
Show[pic /. Polygon[triangles_] :> {EdgeForm[Opacity[0.3]], 
    Polygon[DeleteCases[triangles, _?(someSubsetQ[#, boundaryComponents] &)]]}, 
Graphics3D[{Thick, Line[Table[Append[pt, 0], 
    {pt, Prepend[boundary, Last[boundary]]}]]}]] 

Answer 2

Ce n'est probablement pas la solution optimale, donc je laisserai la question ouverte au cas où quelqu'un aurait une meilleure idée.

Voici comment j'ai résolu le problème que j'ai décrit dans ma question. Tout d'abord, j'ai remplacé les points {x,y,f[x,y]} dans data pour lesquels f[x,y] était complexe par {x,y,None}. Ensuite, la fonction suivante crée une surface 3D hors de mes points de données. Notez ici que data est le résultat de

data = Table[{x,y,f[x,y]},{x,x0,x1,dx},{y,y0,y1,dy}] 

qui est, sans aplatissement (qui fonctionne mieux pour la fonction suivante). La fonction est:

makeSurface[data_] := Module[{len1, len2, polys, a, b, c, d, p}, 
    len1 = Length[data]; 
    len2 = Length[data[[1]]]; 
    polys = Table[ 
    a = data[[i, j]]; 
    b = data[[i + 1, j]]; 
    c = data[[i + 1, j + 1]]; 
    d = data[[i, j + 1]]; 
    p = Select[{a, b, c, d}, #[[3]] =!= None &]; 
    If[Length[p] >= 2, Polygon[p], None], 
    {i, 1, len1 - 1}, {j, 1, len2 - 1}]; 
    Graphics3D[Join[{EdgeForm[None],FaceForm[Directive[GrayLevel[0.5]]]}, 
    Select[Flatten[polys, 1], # =!= None &]]]] 

Le code ci-dessus peut être probablement optimisé mais cela a fonctionné assez bien pour moi.