Mercurial > repos > public > sbplib
changeset 1336:0666629aa183 feature/D2_boundary_opt
Add methods for creating grids with different grid point distributions for each coordinate direction, and also supports constructing periodic grids
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Fri, 13 May 2022 13:26:16 +0200 |
parents | 8d9fc7981796 |
children | bf2554f1825d |
files | +grid/generalCurvilinear.m +multiblock/DefCurvilinear.m +util/get_periodic_grid.m |
diffstat | 3 files changed, 71 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/generalCurvilinear.m Fri May 13 13:26:16 2022 +0200 @@ -0,0 +1,54 @@ +% Creates a curvilinear grid of dimension length(m). +% over the logical domain xi_lim, eta_lim, using a provided +% coordinate mapping. The grid point distribution is +% specified is +% Examples: +% g = grid.generalCurvilinear(mapping, [mx, my], xlim, ylim, {'equidist'}, {'boundaryopt', order, 'accurate'}) +function g = generalCurvilinear(mapping, m, varargin) + n = length(m); + + % Check that parameters matches dimensions + matchingParams = false; + if length(varargin) == 2*n + matchingParams = iscell([varargin{1:2*n}]); + end + assert(matchingParams,'grid:generalCurvilinear:NonMatchingParameters','The number of parameters per dimensions do not match.'); + + X = []; + h = []; + inds_periodic = []; + for i = 1:n + lim = varargin{i}; + opts = varargin{i+n}; + gridtype = opts{1}; + switch gridtype + case 'equidist' + gridgenerator = @()util.get_grid(lim{1},lim{2},m(i)); + case 'boundaryopt' + order = opts{2}; + stencil_type = opts{3}; + switch stencil_type + case {'Accurate','accurate','A','acc'} + gridgenerator = @()sbp.grid.accurateBoundaryOptimizedGrid(lim,m(i),order); + case {'Minimal','minimal','M','min'} + gridgenerator = @()sbp.grid.minimalBoundaryOptimizedGrid(lim,m(i),order); + end + case 'periodic' + gridgenerator = @()util.get_periodic_grid(lim{1},lim{2},m(i)); + inds_periodic = [inds_periodic, i]; + otherwise + error("grid type %s not supported. Must be one of 'equidist', 'boundaryopt', 'periodic'",gridtype); + end + try + [X{i},h(i)] = gridgenerator(); + catch exception % Propagate any errors in the grid generation functions. + msgText = getReport(exception); + error('grid:boundaryOptimizedCurvilinear:InvalidParameter',msgText) + end + end + g = grid.Curvilinear(mapping, X{:}); + g.logic.h = h; + for i = inds_periodic + g.logic.lim{i}{2} = g.logic.lim{i}{2}+h(i); + end +end \ No newline at end of file
--- a/+multiblock/DefCurvilinear.m Sat May 07 10:41:22 2022 +0200 +++ b/+multiblock/DefCurvilinear.m Fri May 13 13:26:16 2022 +0200 @@ -52,10 +52,9 @@ % If a scalar is passed, defer to getGridSizes implemented by subclass % TODO: This forces the interface of subclasses. % Should ms be included in varargin? Figure out bow to do it properly - if length(ms) == 1 + if ~iscell(ms) && length(ms) == 1 ms = obj.getGridSizes(ms); end - if isempty(varargin) || strcmp(varargin{1},'equidist') gridgenerator = @(blockMap,m) grid.equidistantCurvilinear(blockMap, m); elseif strcmp(varargin{1},'boundaryopt') @@ -63,13 +62,15 @@ stenciloption = varargin{3}; gridgenerator = @(blockMap,m) grid.boundaryOptimizedCurvilinear(blockMap,m,{0,1},{0,1},... order,stenciloption); - else - error('No grid type supplied!'); - end - grids = cell(1, obj.nBlocks); - for i = 1:obj.nBlocks + elseif strcmp(varargin{1},'general') + gridgenerator = @(blockMap,m) grid.generalCurvilinear(blockMap,m,{0,1},{0,1},varargin{2:end}); + else + error('No grid type supplied!'); + end + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks grids{i} = gridgenerator(obj.blockMaps{i}.S, ms{i}); - end + end g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/get_periodic_grid.m Fri May 13 13:26:16 2022 +0200 @@ -0,0 +1,8 @@ +% Returns the 1D periodic grid where the last grid point at x_r +% is omitted (it is the same as the the first grid point at x_l) +function [x,h] = get_periodic_grid(x_l,x_r,m) + L = x_r-x_l; + h = L/m; + x = linspace(x_l,x_r,m+1)'; + x = x(1:end-1); +end \ No newline at end of file