Mercurial > repos > public > sbplib
view +multiblock/+domain/Rectangle.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Thu, 17 Jan 2019 10:25:06 +0100 |
parents | 4bb298faa8dc |
children | 375f73edbbd4 b0208b130880 |
line wrap: on
line source
classdef Rectangle < multiblock.Definition properties blockTi % Transfinite interpolation objects used for plotting xlims ylims blockNames % Cell array of block labels nBlocks connections % Cell array specifying connections between blocks boundaryGroups % Structure of boundaryGroups end methods % Creates a divided rectangle % x and y are vectors of boundary and interface positions. % blockNames: cell array of labels. The id is default. function obj = Rectangle(x,y,blockNames) default_arg('blockNames',[]); n = length(y)-1; % number of blocks in the y direction. m = length(x)-1; % number of blocks in the x direction. N = n*m; % number of blocks if ~issorted(x) error('The elements of x seem to be in the wrong order'); end if ~issorted(flip(y)) error('The elements of y seem to be in the wrong order'); end % Dimensions of blocks and number of points blockTi = cell(N,1); xlims = cell(N,1); ylims = cell(N,1); for i = 1:n for j = 1:m p1 = [x(j), y(i+1)]; p2 = [x(j+1), y(i)]; I = flat_index(m,j,i); blockTi{I} = parametrization.Ti.rectangle(p1,p2); xlims{I} = {x(j), x(j+1)}; ylims{I} = {y(i+1), y(i)}; end end % Interface couplings conn = cell(N,N); for i = 1:n for j = 1:m I = flat_index(m,j,i); if i < n J = flat_index(m,j,i+1); conn{I,J} = {'s','n'}; end if j < m J = flat_index(m,j+1,i); conn{I,J} = {'e','w'}; end end end % Block names (id number as default) if isempty(blockNames) obj.blockNames = cell(1, N); for i = 1:N obj.blockNames{i} = sprintf('%d', i); end else assert(length(blockNames) == N); obj.blockNames = blockNames; end nBlocks = N; % Boundary groups boundaryGroups = struct(); nx = m; ny = n; E = cell(1,ny); W = cell(1,ny); S = cell(1,nx); N = cell(1,nx); for i = 1:ny E_id = flat_index(m,nx,i); W_id = flat_index(m,1,i); E{i} = {E_id,'e'}; W{i} = {W_id,'w'}; end for j = 1:nx S_id = flat_index(m,j,ny); N_id = flat_index(m,j,1); S{j} = {S_id,'s'}; N{j} = {N_id,'n'}; end boundaryGroups.E = multiblock.BoundaryGroup(E); boundaryGroups.W = multiblock.BoundaryGroup(W); boundaryGroups.S = multiblock.BoundaryGroup(S); boundaryGroups.N = multiblock.BoundaryGroup(N); boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]); boundaryGroups.WS = multiblock.BoundaryGroup([W,S]); boundaryGroups.WN = multiblock.BoundaryGroup([W,N]); boundaryGroups.ES = multiblock.BoundaryGroup([E,S]); boundaryGroups.EN = multiblock.BoundaryGroup([E,N]); obj.connections = conn; obj.nBlocks = nBlocks; obj.boundaryGroups = boundaryGroups; obj.blockTi = blockTi; obj.xlims = xlims; obj.ylims = ylims; end % Returns a multiblock.Grid given some parameters % ms: cell array of [mx, my] vectors % For same [mx, my] in every block, just input one vector. function g = getGrid(obj, ms, varargin) default_arg('ms',[21,21]) % Extend ms if input is a single vector if (numel(ms) == 2) && ~iscell(ms) m = ms; ms = cell(1,obj.nBlocks); for i = 1:obj.nBlocks ms{i} = m; end end grids = cell(1, obj.nBlocks); for i = 1:obj.nBlocks grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); end g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); end % label is the type of label used for plotting, % default is block name, 'id' show the index for each block. function show(obj, label, gridLines, varargin) default_arg('label', 'name') default_arg('gridLines', false); if isempty('label') && ~gridLines for i = 1:obj.nBlocks obj.blockTi{i}.show(2,2); end axis equal return end if gridLines m = 10; for i = 1:obj.nBlocks obj.blockTi{i}.show(m,m); end end switch label case 'name' labels = obj.blockNames; case 'id' labels = {}; for i = 1:obj.nBlocks labels{i} = num2str(i); end otherwise axis equal return end for i = 1:obj.nBlocks parametrization.Ti.label(obj.blockTi{i}, labels{i}); end axis equal end % Returns the grid size of each block in a cell array % The input parameters are determined by the subclass function ms = getGridSizes(obj, varargin) end end end