## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.


# we generate a structured mesh of size [0, 4]x[0, 1] 
# with a square hole at [1, 1.25]x[0.375, 0.625].

# WARNING: the following code is an attempt to create a structured mesh with 
# an hole using the pkg msh.
# Unfortunately at the moment I've not yet been able to do it an a good way, 
# and the following code will raise an error with the version 1.0.8 of msh.
# To avoid this error, you should comment the line number 132 in the file 
# msh2m_join_structured_mesh.m.

pkg load msh

r = 0.25;
nel = 2;
n1 = nel*8 + 1;
n2 = nel*2 + 1;
n3 = nel*22 + 1;
n4 = nel*3 + 1;
n5 = nel*1 + 1;

x = linspace (0, 1, n1);
y = linspace (0, 0.5 - r/2, n4);
msh1 = msh2m_structured_mesh (x, y, 1, 1:4);

x = linspace (1, 1+r, n2);
msh2 = msh2m_structured_mesh (x, y, 1, 5:8);

x = linspace (1+r, 4, n3);
msh3 = msh2m_structured_mesh (x, y, 1, 9:12);

msh12 = msh2m_join_structured_mesh (msh1, msh2, 2, 8);
msh13 = msh2m_join_structured_mesh (msh12, msh3, 6, 12);


####
x = linspace (0, 1, n1);
y = linspace (0.5 + r/2, 1, n4);
msh7 = msh2m_structured_mesh (x, y, 1, 25:28);

x = linspace (1, 1+r, n2);
msh8 = msh2m_structured_mesh (x, y, 1, 29:32);
msh78 = msh2m_join_structured_mesh (msh7, msh8, 26, 32);

x = linspace (1+r, 4, n3);
msh9 = msh2m_structured_mesh (x, y, 1, 33:36);
msh79 = msh2m_join_structured_mesh (msh78, msh9, 30, 36);


###
x = linspace (0, 1, n1);
y = linspace (0.5 - r/2, 0.5, n5);
msh4 = msh2m_structured_mesh (x, y, 1, 13:16);
msh14 = msh2m_join_structured_mesh (msh13, msh4, 3, 13);

x = linspace (1, 1+r, n2);
msh5 = msh2m_structured_mesh (x, y, 1, 17:20);

x = linspace (1+r, 4, n3);
msh6 = msh2m_structured_mesh (x, y, 1, 21:24);
msh16 = msh2m_join_structured_mesh (msh14, msh6, 10, 21);

for i = 1:size(msh16.e (5, :),2)
  if (msh16.e(5, i) == 15)
    msh16.e(5, i) = 12;
  endif
endfor

####
x = linspace (0, 1, n1);
y = linspace (0.5, 0.5 + r/2, n5);
msh4 = msh2m_structured_mesh (x, y, 1, 37:40);
msh74 = msh2m_join_structured_mesh (msh79, msh4, 25, 39);

x = linspace (1, 1+r, n2);
msh5 = msh2m_structured_mesh (x, y, 1, 41:44);

x = linspace (1+r, 4, n3);
msh6 = msh2m_structured_mesh (x, y, 1, 45:48);
msh76 = msh2m_join_structured_mesh (msh74, msh6, 32, 47);

for i = 1:size(msh76.e (5, :),2)
  if (msh76.e(5, i) == 38)
    msh76.e(5, i) = 35;
  endif
endfor

###########
msh19 = msh2m_join_structured_mesh (msh16, msh76, 12, 35);