Указание домена куба без углов в Matlab с использованием цикла for

Я хотел бы указать домен в Matlab, как показано ниже. Это оказалось сложнее, чем я думал. До сих пор я пробовал два типа подходов:

  1. Подход от малого к большому. Здесь я начинаю с определения куба в середине, а затем добавляю 6 торчащих блоков. Таким образом, у меня получается 7 областей, для каждой из которых я определяю основные уравнения и граничные условия. Кажется, это работает, но, очевидно, очень медленно. Я предпочитаю решение, больше похожее на 2, а именно:

  2. Подход от большого к малому. Здесь я начинаю с куба, который больше фактического домена, и я хочу сообщить Matlab, какие части исключить из домена, например. части в углах. Предположим, что углы от i/j/k==1 до i/j/k==2. И вот, я не уверен, как это сделать. Конкретно происходит следующее: если одна из трех координат равна 1 или nx/ny/nz, диапазон двух других координат должен быть равен 2:(nx/ny/nz-1).

В логических терминах этот второй подход будет выглядеть так:

for i=1:nx
for j=1:ny
for k=1:nz

%x
if (j==1||j==ny) for i=2:(nx-1)
elseif (k==1||k==nz) for i=2:(nx-1)
else for i=1:nx
end

%y
if (i==1||i==nx)for j=2:(ny-b1)
elseif (k==1||k==nz)for j=2:(ny-1) 
else  for j=1:ny
end

%z
if (i==1||i==nx) for k=2:(nz-1)
elseif (j==1||j==ny) for k=2:(nz-1)
else for k=1:nz
end

Я знаю, что это недопустимый синтаксис Matlab из-за for после оператора if, но это должно дать четкое представление о том, чего я хочу. Мой вопрос в том, есть ли другой способ написать это так, чтобы он работал в Matlab?!

В качестве альтернативы я думал о чем-то вроде:

for i=1:((nx/3))
for j=1:ny
for k=1:nz

if(i==(2||(nx-2))&&(k==1||k==nz ||j==1||j==ny))  ux(i,j,k)=BC1;
elseif(i==1&&k~=1&&k~=nz&&j~=1&&j~=ny) ux(i,j,k)=BC2;
elseif(i==1) ux(i,j,k)=0; %here also u and uxx need to be set 0. 
else definition ux
end

Проблема в том, что вам все равно нужно указать i==1 в углах, например. вместо того, чтобы исключать углы из домена, вы включаете углы и устанавливаете все значения, связанные с углами, равными 0. Я предпочитаю действительно их исключать.

Ни один из этих методов не работает хорошо, и мне интересно, нет ли лучшего и более простого метода?

Изображение описываемого домена


person Phoneheads    schedule 05.11.2017    source источник


Ответы (1)


вам просто нужно подготовить бинарную маску ребер, если имя вашей матрицы m:

mask=zeros(size(m));
mask(1,:,end)=1;
mask(end,:,1)=1;
mask(1,:,1)=1;
mask(end,:,end)=1;
mask(:,1,1)=1;
mask(:,1,end)=1;
mask(:,end,1)=1;
mask(:,end,end)=1;
mask(1,1,:)=1;
mask(1,end,:)=1;
mask(end,1,:)=1;
mask(end,end,:)=1;

потом:

m.*~mask

доставит тебе то, что ты хотел...

используя blockPlot из fex, с матрицей 7x9x11

эта фигура была создана с использованием blockPlot из fex, с m=ones(7,9,11);

person bla    schedule 05.11.2017