ListPlot3D данных о «правильной» решетке, содержащейся в «неправильной» области - PullRequest
2 голосов
/ 11 октября 2010

У меня возникает следующая проблема при попытке построить некоторые трехмерные данные в Mathematica.Данные первоначально вычисляются на регулярной решетке, то есть я вычисляю

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

Проблема в том, что f принимает действительные значения только в невыпуклом подмножестве U на плоскости.U на самом деле довольно неприятно: «треугольная» область, где все углы являются острыми (представьте область между тремя одинаковыми кругами, где любые два из них касаются друг друга).Когда я пытаюсь построить data с помощью ListPlot3D, последний наносит на поверхность поверхность над выпуклой оболочкой U.

Мне было интересно, есть ли способ сказать Mathematica ограничить себя только внутри U.Я думал, что, поскольку мои точки уже лежат на «правильной» решетке, это не должно быть слишком сложно, но я пока не нашел решения.

Мне известен вариант RegionFunction, но он очень дорогойвычислить в моем случае, поэтому я пытаюсь найти способ, который использует только уже вычисленные точки в data.

Ответы [ 2 ]

3 голосов
/ 13 октября 2010

Чтобы создать действительно красивую картинку, без использования миллиардов точек, вы можете использовать неравномерное распределение точек, которое включает в себя границы области, которую вы хотите.Вот пример, похожий на то, что вы описываете.Мы начинаем с трех касательных друг к другу окружностей.

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

Mathematica graphics

Мы пишем булеву функцию, которая определяет, будет ли точка в прямоугольнике {-1 / 2,1 / 2}{0, Sqrt [3] / 2} лежит за пределами всех окружностей и используйте это для генерации некоторых точек в интересующей области.

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];

Теперь мы сгенерируем границу.Параметр n определяет, сколько точек мы размещаем на границе.

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];

Давайте посмотрим.

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

Mathematica graphics

Теперь мы определимиспользуйте и ListPlot3D, чтобы попытаться построить его.

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]

Mathematica graphics

Каким-то образом мы должны удалить этот материал, который находится за пределами региона.В этом конкретном примере мы можем использовать тот факт, что функция равна нулю на границе.

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]

Mathematica graphics

Довольно хорошо, но есть несколько проблем рядом с остриямии это определенно не очень общее, так как использовало определенное свойство функции.Если вы исследуете структуру pic, вы обнаружите, что она содержит GraphicsComplex, а первые n точек в первом аргументе этого GraphicsComplex являются именно границей.Вот доказательство:

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

Теперь граница состоит из трех компонентов, и мы хотим удалить любой треугольник, образованный точками, выбранными только из одного из этих компонентов.Следующий код делает это.Обратите внимание, что borderComponents содержит индексы тех точек, которые образуют границу, и someSubsetQ [A, Bs] возвращает true, если A является подмножеством любого из B.Мы хотим удалить те индексы треугольника в многополигоне, которые являются подмножествами одного из borderComponents.Это достигается в следующем коде командой DeleteCases.

О, и давайте добавим немного украшений.

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]]}]]}]]
0 голосов
/ 12 октября 2010

Возможно, это не оптимальное решение, поэтому я буду держать вопрос открытым, если у кого-то есть идея получше.

Вот как я решил проблему, описанную в моем вопросе.Сначала я заменил точки {x,y,f[x,y]} в data, для которых f[x,y] был сложным, на {x,y,None}.Затем следующая функция создает трехмерную поверхность из моих точек данных.Обратите внимание, что data - это результат

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

, то есть без выравнивания (что лучше работает для следующей функции).Функция:

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 &]]]]

Приведенный выше код, вероятно, можно оптимизировать, но он работал достаточно хорошо для меня.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...