Mercurial > repos > public > sbplib
changeset 423:a2cb0d4f4a02 feature/grids
Merge in default.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 07 Feb 2017 15:47:51 +0100 |
| parents | da058ce66876 (diff) b6a5dc423990 (current diff) |
| children | e56dbd9e4196 |
| files | +parametrization/place_label.m +sbp/BlockNorm.m +sbp/Higher.m +sbp/HigherCompatible.m +sbp/HigherCompatibleVariable.m +sbp/HigherPeriodic.m +sbp/Ordinary.m +sbp/Upwind.m +sbp/Variable.m +sbp/blocknorm10.m +sbp/blocknorm4.m +sbp/blocknorm6.m +sbp/blocknorm8.m +sbp/higher2_compatible_halfvariable.m +sbp/higher4.m +sbp/higher4_compatible_halfvariable.m +sbp/higher6.m +sbp/higher6_compatible_halfvariable.m +sbp/higher_compatible2.m +sbp/higher_compatible4.m +sbp/higher_compatible6.m +sbp/ordinary10.m +sbp/ordinary2.m +sbp/ordinary4.m +sbp/ordinary6.m +sbp/ordinary8.m +sbp/upwind2.m +sbp/upwind3.m +sbp/upwind4.m +sbp/upwind5.m +sbp/upwind6.m +sbp/upwind7.m +sbp/upwind8.m +sbp/upwind9.m +sbp/variable4.m |
| diffstat | 81 files changed, 3499 insertions(+), 1319 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/fromMatrix.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,24 @@ +function bm = fromMatrix(A, div) + d1 = div{1}; + d2 = div{2}; + [n, m] = size(A); + if n ~= sum(d1) || m ~= sum(d2) + error('blockmatrix:fromMatrix:NonMatchingDim','The dimensions in div does not sum to the dimensions in A.'); + end + + bm = cell(length(d1), length(d2)); + I = 1; + for i = 1:length(d1) + J = 1; + for j = 1:length(d2) + Asub = A(I:(I + d1(i)-1), J:(J + d2(j)-1)); + if nnz(Asub) == 0 + bm{i,j} = []; + else + bm{i,j} = Asub; + end + J = J + d2(j); + end + I = I + d1(i); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/fromMatrixTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,64 @@ +function tests = fromMatrixTest() + tests = functiontests(localfunctions); +end + +function testErrorNonMatchingDim(testCase) + in = { + {magic(5), {[1 2 3], [4]}}, + {magic(5), {[1 1 1 1 1 1], [5]}}, + {magic(5), {[5], [1 1 1 1 1 1]}}, + {ones(4,2),{[2 3],[2]}}, + {ones(4,2),{[2 2],[3]}}, + }; + + for i = 1:length(in) + testCase.verifyError(@()blockmatrix.fromMatrix(in{i}{:}),'blockmatrix:fromMatrix:NonMatchingDim'); + end +end + +function testFromMatrix(testCase) + cases = { + { + {[],{[],[]}}, + {} + }, + { + { + magic(3), + {[3],[3]} + }, + {magic(3)} + }, + { + { + magic(3), + {[1 1 1],[1 1 1]} + }, + mat2cell(magic(3),[1 1 1],[1 1 1]) + }, + { + { + [17 24 1 8 15; 23 5 7 14 16; 4 6 13 20 22; 10 12 19 21 3; 11 18 25 2 9], + {[1 4],[2 3]} + }, + { + [17 24], [1 8 15]; + [23 5; 4 6; 10 12; 11 18], [7 14 16; 13 20 22; 19 21 3; 25 2 9]; + }; + }, + }; + for i = 1:length(cases) + out = convertToFull(blockmatrix.fromMatrix(cases{i}{1}{:})); + expected = cases{i}{2}; + testCase.verifyEqual(out,expected); + end +end + +function C = convertToFull(C) + [N,M] = size(C); + for i = 1:N + for j = 1:M + C{i,j} = full(C{i,j}); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/getDivision.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,37 @@ +function div = getDivision(bm) + if ~blockmatrix.isBlockmatrix(bm) + error('blockmatrix:getDivision:NotABlockmatrix', 'Input is not a blockmatrix'); + end + + if isempty(bm) + div = {[],[]}; + return + end + + div = {row_height(bm),col_width(bm)}; +end + + +function m = col_width(C) + m = zeros(1,size(C,2)); + for j = 1:size(C,2) + for i = 1:size(C,1) + if isempty(C{i,j}) + continue + end + m(j) = size(C{i,j},2); + end + end +end + +function n = row_height(C) + n = zeros(1,size(C,1)); + for i = 1:size(C,1) + for j = 1:size(C,2) + if isempty(C{i,j}) + continue + end + n(i) = size(C{i,j},1); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/getDivisionTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,69 @@ +function tests = getDivisionTest() + tests = functiontests(localfunctions); +end + +function testError(testCase) + cases = { + magic(3), + {[2 2 2];{1,2}}, + {[2 2 2];[1 2]}, + {[2; 2; 2], [1; 2]}, + }; + + for i =1:length(cases) + testCase.verifyError(@()blockmatrix.getDivision(cases{i}), 'blockmatrix:getDivision:NotABlockmatrix') + end +end + +function testGetDivision(testCase) + cases = { + { + {}, + {[],[]}; + }, + { + { + [2 2; 2 1], [1; 2]; + [2 2], [1] + }, + {[2 1], [2 1]} + }, + { + { + [2 2; 2 1], []; + [2 2], [1] + }, + {[2 1], [2 1]} + }, + { + { + [2 2; 2 1], []; + [2 2], [] + }, + {[2 1], [2 0]} + }, + { + { + [2 2; 2 1], [1; 2]; + [], [] + }, + {[2 0], [2 1]} + }, + { + { + [2 2; 2 1]; + [2 2] + }, + {[2 1], 2} + }, + }; + + for i = 1:length(cases) + in = cases{i}{1}; + out = blockmatrix.getDivision(in); + expected = cases{i}{2}; + testCase.verifyEqual(out, expected); + end +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/isBlockmatrix.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,59 @@ +function b = isBlockmatrix(bm) + if ~iscell(bm) + b = false; + return + end + + % Make sure all blocks are numerica matrices + for i = 1:length(bm) + if ~isnumeric(bm{i}) + b = false; + return + end + end + + [N,M] = size(bm); + % Make sure column dimensions agree + for i = 1:N + d = []; + for j = 1:M + d_ij = size(bm{i,j},1); + if d_ij == 0 + continue + end + + if isempty(d) + d = d_ij; + continue + end + + if d ~= d_ij + b = false; + return + end + end + end + + % Make sure row dimensions agree + for j = 1:M + d = []; + for i = 1:N + d_ij = size(bm{i,j},2); + if d_ij == 0 + continue + end + + if isempty(d) + d = d_ij; + continue + end + + if d ~= d_ij + b = false; + return + end + end + end + + b = true; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/isBlockmatrixTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,63 @@ +function tests = isBlockmatrixTest() + tests = functiontests(localfunctions); +end + +function testIsBlockmatrix(testCase) + cases = { + { + magic(3), + false % Must be a cell array + } + { + {[2 2 2];{1,2}}, + false % All elements of the cell matrix must be regular matrices + }, + { + {[2 2 2];[1 2]}, + false % Row dimensions must match + }, + { + {[2; 2; 2], [1; 2]}, + false % Column dimensions must match + }, + { + {}, + true % An empty matrix is a matrix too + }, + { + { + [2 2; 2 1], [1; 2]; + [2 2], [1] + }, + true % A simple valid one + }, + { + { + [2 2; 2 1], []; + [2 2], [1] + }, + true % Empty blocks assumed to be zero and match dimensions + }, + { + { + [2 2; 2 1], []; + [2 2], [] + }, + true % Empty blocks allowed. + }, + { + { + [2 2; 2 1], [1; 2]; + [], [] + }, + true % Empty blocks allowed. + }, + }; + + for i = 1:length(cases) + in = cases{i}{1}; + out = blockmatrix.isBlockmatrix(in); + expected = cases{i}{2}; + testCase.verifyEqual(out, expected, sprintf('Should return %d for %s', expected, toString(in))); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/isDivision.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,34 @@ +function b = isDivision(div) + % Make sure it is a cellarray + if ~iscell(div) + b = false; + return + end + + % Make sure it has the right shape + if numel(div) ~= 2 + b = false; + return + end + + if ~isDivisionVector(div{1}) || ~isDivisionVector(div{2}) + b = false; + return + end + + b = true; +end + +function b = isDivisionVector(v) + if isempty(v) + b = false; + return + end + + if any(v <= 0) + b = false; + return + end + + b = true; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/isDivisionTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,23 @@ +function tests = isDivisionTest() + tests = functiontests(localfunctions); +end + +function testIsDivision(testCase) + cases = { + {{[2 2 2],[1 2]} ,true}, + {{[1 2],[1 0]} ,false}, + {{[0 2],[1 1]} ,false}, + {{[1 2],[]} ,false}, + {{[1 2],[1]} ,true}, + {{[1 2],[1], [1 2 3]} ,false}, + {{[1 2 3]} ,false}, + {[1 2] ,false}, + }; + + for i = 1:length(cases) + in = cases{i}{1}; + out = blockmatrix.isDivision(in); + expected = cases{i}{2}; + testCase.verifyEqual(out, expected, sprintf('Should return %d for %s', expected, toString(in))); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/toMatrix.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,26 @@ +function A = toMatrix(bm) + if ~blockmatrix.isBlockmatrix(bm) + error('blockmatrix:toMatrix:NotABlockmatrix', 'Input is not a blockmatrix'); + end + + div = blockmatrix.getDivision(bm); + n = div{1}; + m = div{2}; + + N = sum(n); + M = sum(m); + + A = sparse(N,M); + + n_ind = [0 cumsum(n)]; + m_ind = [0 cumsum(m)]; + + for i = 1:size(bm,1) + for j = 1:size(bm,2) + if isempty(bm{i,j}) + continue + end + A(n_ind(i)+1:n_ind(i+1),m_ind(j)+1:m_ind(j+1)) = bm{i,j}; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/toMatrixTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,62 @@ +function tests = toMatrixTest() + tests = functiontests(localfunctions); +end + +function testError(testCase) + testCase.verifyError(@()blockmatrix.toMatrix([]), 'blockmatrix:toMatrix:NotABlockmatrix') +end + +function testToMatrix(testCase) + cases = { + { + {}, + [], + }, + { + {1, 2; 3, 4}, + [1,2; 3,4], + } + { + { + [2 2; 2 1], [1; 2]; + [2 2], [1] + }, + [2 2 1; + 2 1 2; + 2 2 1], + }, + { + { + [2 2; 2 1], []; + [2 2], [1] + }, + [2 2 0; + 2 1 0; + 2 2 1], + }, + { + { + [2 2; 2 1], []; + [2 2], [] + }, + [2 2; + 2 1; + 2 2], + }, + { + { + [2 2; 2 1], [1; 2]; + [], [] + }, + [2 2 1; + 2 1 2], + }, + }; + + for i = 1:length(cases) + in = cases{i}{1}; + out = full(blockmatrix.toMatrix(in)); + expected = cases{i}{2}; + testCase.verifyEqual(out, expected); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/zero.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,16 @@ +% Creates a block matrix according to the division with zeros everywhere. +function bm = zero(div) + n = div{1}; + m = div{2}; + + N = length(n); + M = length(m); + + bm = cell(N,M); + + for i = 1:N + for j = 1:M + bm{i,j} = sparse(n(i),m(j)); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+blockmatrix/zeroTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,52 @@ +function tests = zeroTest() + tests = functiontests(localfunctions); +end + +function testZero(testCase) + cases = { + { + {[],[]}, + {}, + }, + { + {0,0}, + {[]}; + }, + { + {1,1}, + {0}; + }, + { + {2,1}, + {[0; 0]}; + }, + { + {1,2}, + {[0 0]}; + }, + { + {[1 2],2}, + {[0 0];[0 0; 0 0]}; + }, + { + {[1 2],[2 1]}, + {[0 0],[0];[0 0; 0 0],[0; 0]}; + }, + }; + + for i = 1:length(cases) + out = convertToFull(blockmatrix.zero(cases{i}{1})); + expected = cases{i}{2}; + testCase.verifyEqual(out,expected); + end +end + + +function C = convertToFull(C) + [N,M] = size(C); + for i = 1:N + for j = 1:M + C{i,j} = full(C{i,j}); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Cartesian.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,189 @@ +classdef Cartesian < grid.Structured + properties + n % Number of points in the grid + d % Number of dimensions + m % Number of points in each direction + x % Cell array of vectors with node placement for each dimension. + h % Spacing/Scaling + end + + % General d dimensional grid with n points + methods + % Creates a cartesian grid given vectors conatining the coordinates + % in each direction + function obj = Cartesian(varargin) + obj.d = length(varargin); + + for i = 1:obj.d + obj.x{i} = varargin{i}; + obj.m(i) = length(varargin{i}); + end + + obj.n = prod(obj.m); + if obj.n == 0 + error('grid:Cartesian:EmptyGrid','Input parameter gives an empty grid.') + end + + obj.h = []; + end + % n returns the number of points in the grid + function o = N(obj) + o = obj.n; + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.d; + end + + function m = size(obj) + m = obj.m; + end + + % points returns a n x d matrix containing the coordianets for all points. + % points are ordered according to the kronecker product with X*Y*Z + function X = points(obj) + X = zeros(obj.n, obj.d); + + for i = 1:obj.d + if iscolumn(obj.x{i}) + c = obj.x{i}; + else + c = obj.x{i}'; + end + + m_before = prod(obj.m(1:i-1)); + m_after = prod(obj.m(i+1:end)); + + X(:,i) = kr(ones(m_before,1),c,ones(m_after,1)); + end + end + + % matrices returns a cell array with coordinates in matrix form. + % For 2d case these will have to be transposed to work with plotting routines. + function X = matrices(obj) + + if obj.d == 1 % There is no 1d matrix data type in matlab, handle special case + X{1} = reshape(obj.x{1}, [obj.m 1]); + return + end + + X = cell(1,obj.d); + for i = 1:obj.d + s = ones(1,obj.d); + s(i) = obj.m(i); + + t = reshape(obj.x{i},s); + + s = obj.m; + s(i) = 1; + X{i} = repmat(t,s); + end + end + + function h = scaling(obj) + if isempty(obj.h) + error('grid:Cartesian:NoScalingSet', 'No scaling set') + end + + h = obj.h; + end + + % Restricts the grid function gf on obj to the subgrid g. + % Only works for even multiples + function gf = restrictFunc(obj, gf, g) + m1 = obj.m; + m2 = g.m; + + % Check the input + if prod(m1) ~= numel(gf) + error('grid:Cartesian:restrictFunc:NonMatchingFunctionSize', 'The grid function has to few or too many points.'); + end + + if ~all(mod(m1-1,m2-1) == 0) + error('grid:Cartesian:restrictFunc:NonMatchingGrids', 'Only integer downsamplings are allowed'); + end + + % Calculate stride for each dimension + stride = (m1-1)./(m2-1); + + % Create downsampling indecies + I = {}; + for i = 1:length(m1) + I{i} = 1:stride(i):m1(i); + end + + gf = reshapeRowMaj(gf, m1); + gf = gf(I{:}); + gf = reshapeRowMaj(gf, prod(m2)); + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + error('grid:Cartesian:NotImplemented') + end + + % Return the names of all boundaries in this grid. + function bs = getBoundaryNames(obj) + switch obj.D() + case 1 + bs = {'l', 'r'}; + case 2 + bs = {'w', 'e', 's', 'n'}; + case 3 + bs = {'w', 'e', 's', 'n', 'd', 'u'}; + otherwise + error('not implemented'); + end + end + + % Return coordinates for the given boundary + function X = getBoundary(obj, name) + % In what dimension is the boundary? + switch name + case {'l', 'r', 'w', 'e'} + D = 1; + case {'s', 'n'} + D = 2; + case {'d', 'u'} + D = 3; + otherwise + error('not implemented'); + end + + % At what index is the boundary? + switch name + case {'l', 'w', 's', 'd'} + index = 1; + case {'r', 'e', 'n', 'u'} + index = obj.m(D); + otherwise + error('not implemented'); + end + + + + I = cell(1, obj.d); + for i = 1:obj.d + if i == D + I{i} = index; + else + I{i} = ':'; + end + end + + % Calculate size of result: + m = obj.m; + m(D) = []; + N = prod(m); + + X = zeros(N, obj.d); + + coordMat = obj.matrices(); + for i = 1:length(coordMat) + Xtemp = coordMat{i}(I{:}); + X(:,i) = reshapeRowMaj(Xtemp, [N,1]); + end + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/CartesianTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,318 @@ +function tests = CartesianTest() + tests = functiontests(localfunctions); +end + + +function testWarningEmptyGrid(testCase) + in = { + {[]}, + {[],[1]}, + {[1],[2], []}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.Cartesian(in{i}{:}),'grid:Cartesian:EmptyGrid'); + end +end + +function testN(testCase) + in = { + {[1 2 3]}, + {[1 2 3],[1 2]}, + {[1 2 3],[1 2 3]}, + {[1 2 3],[1 2 3], [1]}, + {[1 2 3],[1 2 3], [1 3 4]}, + }; + + out = [3,6,9,9,27]; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.N(),out(i)); + end +end + + +function testD(testCase) + in = { + {[1 2 3]}, + {[1 2 3],[1 2]}, + {[1 2 3],[1 2 3]}, + {[1 2 3],[1 2 3], [1]}, + {[1 2 3],[1 2 3], [1 3 4]}, + }; + + out = [1,2,2,3,3]; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.D(),out(i)); + end +end + +function testSize(testCase) + in = { + {[1 2 3]}, + {[1 2 3],[1 2]}, + {[1 2 3],[1 2 3]}, + {[1 2 3],[1 2 3], [1]}, + {[1 2 3],[1 2 3], [1 3 4]}, + }; + + out = { + [3], + [3 2], + [3 3], + [3 3 1], + [3 3 3], + }; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.size(),out{i}); + end +end + +function testPoints(testCase) + in = { + {[1 2]}, + {[1 2],[3 4]}, + {[1 2],[3 4], [5 6]}, + }; + + out = { + [[1; 2]], + [[1; 1; 2; 2],[3; 4; 3; 4]], + [[1; 1; 1; 1; 2; 2; 2; 2],[3; 3; 4; 4; 3; 3; 4; 4],[ 5; 6; 5; 6; 5; 6; 5; 6]], + }; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}); + end +end + +function testMatrices(testCase) + in = { + {[1 2]}, + {[1 2],[3 4]}, + {[1 2],[3 4], [5 6]}, + }; + + out{1}{1} = [1; 2]; + + out{2}{1} = [1, 1; 2, 2]; + out{2}{2} = [3, 4; 3, 4]; + + out{3}{1}(:,:,1) = [1, 1; 2, 2]; + out{3}{1}(:,:,2) = [1, 1; 2, 2]; + + out{3}{2}(:,:,1) = [3, 4; 3, 4]; + out{3}{2}(:,:,2) = [3, 4; 3, 4]; + + out{3}{3}(:,:,1) = [5, 5; 5, 5]; + out{3}{3}(:,:,2) = [6, 6; 6, 6]; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.matrices(),out{i}); + end +end + + +function testRestrictFuncInvalidInput(testCase) + inG1 = { + {[1 2 3 4 5]}, + {[1 2 3],[4 5 6 7 8]}, + {[1 2 3],[4 5 6 7 8]}, + {[1 2 3],[4 5 6 7 8]}, + }; + + inG2 = { + {[1 3 4 5]}, + {[1 3],[4 5 6 8]}, + {[1 3],[4 6 8]}, + {[1 3],[4 6 8]}, + }; + + inGf = { + [1; 2; 3; 4; 5], + [14; 15; 16; 17; 18; 24; 25; 26; 27; 28; 34; 35; 36; 37; 38]; + [14; 15; 16; 17; 18; 24; 25; 26; 27; 28; 34; 35; 36]; + [14; 15; 16; 17; 18; 24; 25; 26; 27; 28; 34; 35; 36; 37; 38; 39; 40]; + }; + + out = { + 'grid:Cartesian:restrictFunc:NonMatchingGrids', + 'grid:Cartesian:restrictFunc:NonMatchingGrids', + 'grid:Cartesian:restrictFunc:NonMatchingFunctionSize', + 'grid:Cartesian:restrictFunc:NonMatchingFunctionSize', + }; + + for i = 1:length(inG1) + g1 = grid.Cartesian(inG1{i}{:}); + g2 = grid.Cartesian(inG2{i}{:}); + testCase.verifyError(@()g1.restrictFunc(inGf{i},g2),out{i}); + end +end + +function testRestrictFunc(testCase) + inG1 = { + {[1 2 3 4 5]}, + {[1 2 3],[4 5 6 7 8]}, + }; + + inG2 = { + {[1 3 5]}, + {[1 3],[4 6 8]}, + }; + + inGf = { + [1; 2; 3; 4; 5], + [14; 15; 16; 17; 18; 24; 25; 26; 27; 28; 34; 35; 36; 37; 38]; + }; + + outGf = { + [1; 3; 5], + [14; 16; 18; 34; 36; 38]; + }; + + for i = 1:length(inG1) + g1 = grid.Cartesian(inG1{i}{:}); + g2 = grid.Cartesian(inG2{i}{:}); + testCase.verifyEqual(g1.restrictFunc(inGf{i}, g2), outGf{i}); + end +end + +function testScaling(testCase) + in = {[1 2 3], [1 2]}; + g = grid.Cartesian(in{:}); + + testCase.verifyError(@()g.scaling(),'grid:Cartesian:NoScalingSet'); + + g.h = [2 1]; + testCase.verifyEqual(g.scaling(),[2 1]); + +end + + +function testGetBoundaryNames(testCase) + in = { + {[1 2 3]}, + {[1 2 3], [4 5]}, + {[1 2 3], [4 5], [6 7 8]}, + }; + + out = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.getBoundaryNames(), out{i}); + end +end + +function testGetBoundary(testCase) + grids = { + {[1 2 3]}, + {[1 2 3], [4 5]}, + {[1 2 3], [4 5], [6 7 8]}, + }; + + boundaries = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + + % 1d + out{1,1} = 1; + out{1,2} = 3; + + % 2d + out{2,1} = [ + 1,4; + 1,5; + ]; + out{2,2} = [ + 3,4; + 3,5; + ]; + out{2,3} = [ + 1,4; + 2,4; + 3,4; + ]; + out{2,4} = [ + 1,5; + 2,5; + 3,5; + ]; + + % 3d + out{3,1} = [ + 1,4,6; + 1,4,7; + 1,4,8; + 1,5,6; + 1,5,7; + 1,5,8; + ]; + out{3,2} = [ + 3,4,6; + 3,4,7; + 3,4,8; + 3,5,6; + 3,5,7; + 3,5,8; + ]; + out{3,3} = [ + 1,4,6; + 1,4,7; + 1,4,8; + 2,4,6; + 2,4,7; + 2,4,8; + 3,4,6; + 3,4,7; + 3,4,8; + ]; + out{3,4} = [ + 1,5,6; + 1,5,7; + 1,5,8; + 2,5,6; + 2,5,7; + 2,5,8; + 3,5,6; + 3,5,7; + 3,5,8; + ]; + out{3,5} = [ + 1,4,6; + 1,5,6; + 2,4,6; + 2,5,6; + 3,4,6; + 3,5,6; + ]; + out{3,6} = [ + 1,4,8; + 1,5,8; + 2,4,8; + 2,5,8; + 3,4,8; + 3,5,8; + ]; + + for ig = 1:length(grids) + g = grid.Cartesian(grids{ig}{:}); + for ib = 1:length(boundaries{ig}) + testCase.verifyEqual(g.getBoundary(boundaries{ig}{ib}), out{ig,ib}); + end + end +end
--- a/+grid/Curve.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,388 +0,0 @@ -classdef Curve - properties - g - gp - transformation - end - - methods - %TODO: - % Concatenation of curves - % Subsections of curves - % Stretching of curve parameter - done for arc length. - % Curve to cell array of linesegments - - % Returns a curve object. - % g -- curve parametrization for parameter between 0 and 1 - % gp -- parametrization of curve derivative - function obj = Curve(g,gp,transformation) - default_arg('gp',[]); - default_arg('transformation',[]); - p_test = g(0); - assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.'); - - if ~isempty(transformation) - transformation.base_g = g; - transformation.base_gp = gp; - [g,gp] = grid.Curve.transform_g(g,gp,transformation); - end - - obj.g = g; - obj.gp = gp; - obj.transformation = transformation; - - end - - function n = normal(obj,t) - assert(~isempty(obj.gp),'Curve has no derivative!'); - deriv = obj.gp(t); - normalization = sqrt(sum(deriv.^2,1)); - n = [-deriv(2,:)./normalization; deriv(1,:)./normalization]; - end - - - % Plots a curve g(t) for 0<t<1, using n points. Returns a handle h to the plotted curve. - % h = plot_curve(g,n) - function h = plot(obj,n) - default_arg('n',100); - - t = linspace(0,1,n); - p = obj.g(t); - h = line(p(1,:),p(2,:)); - end - - function h= plot_normals(obj,l,n,m) - default_arg('l',0.1); - default_arg('n',10); - default_arg('m',100); - t_n = linspace(0,1,n); - - normals = obj.normal(t_n)*l; - - n0 = obj.g(t_n); - n1 = n0 + normals; - - h = line([n0(1,:); n1(1,:)],[n0(2,:); n1(2,:)]); - set(h,'Color',Color.red); - obj.plot(m); - end - - function h= show(obj,name) - p = obj.g(1/2); - n = obj.normal(1/2); - p = p + n*0.1; - - % Add arrow - - h = text(p(1),p(2),name); - h.HorizontalAlignment = 'center'; - h.VerticalAlignment = 'middle'; - - obj.plot(); - end - % Shows curve with name and arrow for direction. - - % Length of arc from parameter t0 to t1 (which may be vectors). - % Computed using derivative. - function L = arcLength(obj,t0,t1) - assert(~isempty(obj.gp),'Curve has no derivative!'); - speed = @(t) sqrt(sum(obj.gp(t).^2,1)); - L = integral_vec(speed,t0,t1); - end - - % Creates the arc length parameterization of a curve. - % N -- number of points used to approximate the arclength function - function curve = arcLengthParametrization(obj,N) - default_arg('N',100); - assert(~isempty(obj.gp),'Curve has no derivative!'); - - % Construct arcLength function using splines - tvec = linspace(0,1,N); - arcVec = obj.arcLength(0,tvec); - tFunc = spline(arcVec,tvec); % t as a function of arcLength - L = obj.arcLength(0,1); - arcPar = @(s) tFunc(s*L); - - % New function and derivative - g_new = @(t)obj.g(arcPar(t)); - gp_new = @(t) normalize(obj.gp(arcPar(t))); - curve = grid.Curve(g_new,gp_new); - - end - - % how to make it work for methods without returns - function p = subsref(obj,S) - %Should i add error checking here? - %Maybe if you want performance you fetch obj.g and then use that - switch S(1).type - case '()' - p = obj.g(S.subs{1}); - % case '.' - - % p = obj.(S.subs); - otherwise - p = builtin('subsref',obj,S); - % error() - end - end - - - %% TRANSFORMATION OF A CURVE - function D = reverse(C) - % g = C.g; - % gp = C.gp; - % D = grid.Curve(@(t)g(1-t),@(t)-gp(1-t)); - D = C.transform([],[],-1); - end - - function D = transform(C,A,b,flip) - default_arg('A',[1 0; 0 1]); - default_arg('b',[0; 0]); - default_arg('flip',1); - if isempty(C.transformation) - g = C.g; - gp = C.gp; - transformation.A = A; - transformation.b = b; - transformation.flip = flip; - else - g = C.transformation.base_g; - gp = C.transformation.base_gp; - A_old = C.transformation.A; - b_old = C.transformation.b; - flip_old = C.transformation.flip; - - transformation.A = A*A_old; - transformation.b = A*b_old + b; - transformation.flip = flip*flip_old; - end - - D = grid.Curve(g,gp,transformation); - - end - - function D = translate(C,a) - g = C.g; - gp = C.gp; - - % function v = g_fun(t) - % x = g(t); - % v(1,:) = x(1,:)+a(1); - % v(2,:) = x(2,:)+a(2); - % end - - % D = grid.Curve(@g_fun,gp); - - D = C.transform([],a); - end - - function D = mirror(C, a, b) - assert_size(a,[2,1]); - assert_size(b,[2,1]); - - g = C.g; - gp = C.gp; - - l = b-a; - lx = l(1); - ly = l(2); - - - % fprintf('Singular?\n') - - A = [lx^2-ly^2 2*lx*ly; 2*lx*ly ly^2-lx^2]/(l'*l); - - % function v = g_fun(t) - % % v = a + A*(g(t)-a) - % x = g(t); - - % ax1 = x(1,:)-a(1); - % ax2 = x(2,:)-a(2); - % v(1,:) = a(1)+A(1,:)*[ax1;ax2]; - % v(2,:) = a(2)+A(2,:)*[ax1;ax2]; - % end - - % function v = gp_fun(t) - % v = A*gp(t); - % end - - % D = grid.Curve(@g_fun,@gp_fun); - - % g = A(g-a)+a = Ag - Aa + a; - b = - A*a + a; - D = C.transform(A,b); - - end - - function D = rotate(C,a,rad) - assert_size(a, [2,1]); - assert_size(rad, [1,1]); - g = C.g; - gp = C.gp; - - - A = [cos(rad) -sin(rad); sin(rad) cos(rad)]; - - % function v = g_fun(t) - % % v = a + A*(g(t)-a) - % x = g(t); - - % ax1 = x(1,:)-a(1); - % ax2 = x(2,:)-a(2); - % v(1,:) = a(1)+A(1,:)*[ax1;ax2]; - % v(2,:) = a(2)+A(2,:)*[ax1;ax2]; - % end - - % function v = gp_fun(t) - % v = A*gp(t); - % end - - % D = grid.Curve(@g_fun,@gp_fun); - - - % g = A(g-a)+a = Ag - Aa + a; - b = - A*a + a; - D = C.transform(A,b); - end - end - - methods (Static) - - % Computes the derivative of g: R -> R^2 using an operator D1 - function gp_out = numericalDerivative(g,D1) - m = length(D1); - t = linspace(0,1,m); - gVec = g(t)'; - gpVec = (D1*gVec)'; - - gp1_fun = spline(t,gpVec(1,:)); - gp2_fun = spline(t,gpVec(2,:)); - gp_out = @(t) [gp1_fun(t);gp2_fun(t)]; - end - - function obj = line(p1, p2) - - function v = g_fun(t) - v(1,:) = p1(1) + t.*(p2(1)-p1(1)); - v(2,:) = p1(2) + t.*(p2(2)-p1(2)); - end - g = @g_fun; - - obj = grid.Curve(g); - end - - function obj = circle(c,r,phi) - default_arg('phi',[0; 2*pi]) - default_arg('c',[0; 0]) - default_arg('r',1) - - function v = g_fun(t) - w = phi(1)+t*(phi(2)-phi(1)); - v(1,:) = c(1) + r*cos(w); - v(2,:) = c(2) + r*sin(w); - end - - function v = g_fun_deriv(t) - w = phi(1)+t*(phi(2)-phi(1)); - v(1,:) = -(phi(2)-phi(1))*r*sin(w); - v(2,:) = (phi(2)-phi(1))*r*cos(w); - end - - obj = grid.Curve(@g_fun,@g_fun_deriv); - end - - function obj = bezier(p0, p1, p2, p3) - function v = g_fun(t) - v(1,:) = (1-t).^3*p0(1) + 3*(1-t).^2.*t*p1(1) + 3*(1-t).*t.^2*p2(1) + t.^3*p3(1); - v(2,:) = (1-t).^3*p0(2) + 3*(1-t).^2.*t*p1(2) + 3*(1-t).*t.^2*p2(2) + t.^3*p3(2); - end - - function v = g_fun_deriv(t) - v(1,:) = 3*(1-t).^2*(p1(1)-p0(1)) + 6*(1-t).*t*(p2(1)-p1(1)) + 3*t.^2*(p3(1)-p2(1)); - v(2,:) = 3*(1-t).^2*(p1(2)-p0(2)) + 6*(1-t).*t*(p2(2)-p1(2)) + 3*t.^2*(p3(2)-p2(2)); - end - - obj = grid.Curve(@g_fun,@g_fun_deriv); - end - - - function [g_out,gp_out] = transform_g(g,gp,tr) - A = tr.A; - b = tr.b; - flip = tr.flip; - - function v = g_fun_noflip(t) - % v = A*g + b - x = g(t); - - v(1,:) = A(1,:)*x+b(1); - v(2,:) = A(2,:)*x+b(2); - end - - function v = g_fun_flip(t) - % v = A*g + b - x = g(1-t); - - v(1,:) = A(1,:)*x+b(1); - v(2,:) = A(2,:)*x+b(2); - end - - - switch flip - case 1 - g_out = @g_fun_noflip; - gp_out = @(t)A*gp(t); - case -1 - g_out = @g_fun_flip; - gp_out = @(t)-A*gp(1-t); - end - end - - end -end - - - -function g_norm = normalize(g0) - g1 = g0(1,:); - g2 = g0(2,:); - normalization = sqrt(sum(g0.^2,1)); - g_norm = [g1./normalization; g2./normalization]; -end - -function I = integral_vec(f,a,b) -% Wrapper around the built-in function integral that -% handles multiple limits. - - Na = length(a); - Nb = length(b); - assert(Na == 1 || Nb == 1 || Na==Nb,... - 'a and b must have same length, unless one is a scalar.'); - - if(Na>Nb); - I = zeros(size(a)); - for i = 1:Na - I(i) = integral(f,a(i),b); - end - elseif(Nb>Na) - I = zeros(size(b)); - for i = 1:Nb - I(i) = integral(f,a,b(i)); - end - else - I = zeros(size(b)); - for i = 1:Nb - I(i) = integral(f,a(i),b(i)); - end - end -end - -% Returns a function handle to the spline. -function f = spline(tval,fval,spline_order) - default_arg('spline_order',4); - [m,~] = size(tval); - assert(m==1,'Need row vectors.'); - - f_spline = spapi( optknt(tval,spline_order), tval, fval ); - f = @(t) fnval(f_spline,t); -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Curvilinear.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,169 @@ +classdef Curvilinear < grid.Structured & grid.Mapped + properties + logic % Grid of Logical domain + coords % N x D matrix with coordinates of each point in the physical domain + end + + methods + % Creates a curvilinear grid. + % Ex: grid.Curvilinear(mapping, xi, eta, ...) + % mapping -- either a matrix or a cell array with physical coordinates. + % A matrix should be a grid function (N*D x 1 vector) or a N x D + % A cell array should be a 1 x D cell array with either N x 1 vectors + % or matrices of the same dimesions as the logical grid. + % xi, eta, ... -- are the coordinate positions of the cartesian logical grid. + function obj = Curvilinear(mapping, varargin) + xi = varargin; + obj.logic = grid.Cartesian(xi{:}); + + % If mapping is a function evaluate it + if isa(mapping, 'function_handle') + mapping = grid.evalOn(obj.logic, mapping); + end + + D = obj.logic.D(); + N = obj.logic.N(); + + obj.coords = zeros(N,D); + + if iscell(mapping) + obj.coords = cellMappingToCoords(mapping, N, D, obj.logic.m); + elseif isnumeric(mapping) + obj.coords = matrixMappingToCoords(mapping, N, D); + else + error('grid:Curvilinear:Curvilinear','mapping must be a matrix or a cell array.'); + end + end + + function m = size(obj) + m = obj.logic.size(); + end + + % logicalGrid returns the domain grid of the mapping. + function g = logicalGrid(obj) + g = obj.logic; + end + + % mapping returns the mapped coordinates as a grid.Function + function m = mapping(obj); + m = obj.coords; + end + + % n returns the number of points in the grid + function o = N(obj) + o = obj.logic.N(); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.logic.D(); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = obj.coords; + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + gf = obj.logic.restrictFunc(gf, g.baseGrid()); + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + gf = obj.logic.projectFunc(gf,g.baseGrid()); + end + + function h = scaling(obj) + if isempty(obj.logic.h) + error('grid:Curvilinear:NoScalingSet','No scaling set'); + end + h = obj.logic.h; + end + + % Return the names of all boundaries in this grid. + function bs = getBoundaryNames(obj) + bs = obj.logic.getBoundaryNames(); + end + + % Return coordinates for the given boundary + function X = getBoundary(obj, name) + % In what dimension is the boundary? + switch name + case {'l', 'r', 'w', 'e'} + D = 1; + case {'s', 'n'} + D = 2; + case {'d', 'u'} + D = 3; + otherwise + error('not implemented'); + end + + % At what index is the boundary? + switch name + case {'l', 'w', 's', 'd'} + index = 1; + case {'r', 'e', 'n', 'u'} + index = obj.logic.m(D); + otherwise + error('not implemented'); + end + + + + I = cell(1, obj.D); + for i = 1:obj.D + if i == D + I{i} = index; + else + I{i} = ':'; + end + end + + % Calculate size of result: + m = obj.logic.m; + m(D) = []; + N = prod(m); + + X = zeros(N, obj.D); + + p = obj.points; + for i = 1:obj.D() + coordMat{i} = reshapeRowMaj(p(:,i), obj.logic.m); + end + + for i = 1:length(coordMat) + Xtemp = coordMat{i}(I{:}); + X(:,i) = reshapeRowMaj(Xtemp, [N,1]); + end + end + end +end + + +function coords = cellMappingToCoords(mapping, N, D, m) + if ~isequal(size(mapping),[1 D]) + error('grid:Curvilinear:Curvilinear','The cell array must be a 1xD array.'); + end + + if isequal(size(mapping{1}),[N 1]) + coords = cell2mat(mapping); + elseif isequal(size(mapping{1}), m) + for i = 1:length(mapping) + coords(:,i) = reshapeRowMaj(mapping{i}, [N 1]); + end + else + error('grid:Curvilinear:Curvilinear','The matrix must have size [N 1] or the same dimension as the grid. Actual: %s', toString(m)); + end +end + +function coords = matrixMappingToCoords(mapping, N, D) + if isequal(size(mapping), [N, D]) + coords = mapping; + elseif isequal(size(mapping), [N*D, 1]) + coords = reshapeRowMaj(mapping,[N D]); + else + error('grid:Curvilinear:Curvilinear','A matrix mapping must be of size [N D] or [N*D 1].'); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/CurvilinearTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,131 @@ +function tests = CurvilinearTest() + tests = functiontests(localfunctions); +end + +function testMappingInputGridFunction(testCase) + in = { + {{1:10}, @(x) exp(x)}, + {{1:10,1:6}, @(x,y) [exp(x+y); exp(x-y)]}, + {{1:10,1:5,1:7}, @(x,y,z)[exp(x+y+z); exp(x-y-z); 2+x+y-z]}, + }; + + out = { + [10, 1]; + [10*6, 2]; + [10*5*7, 3]; + }; + + + % How to test this? Just make sure it runs without errors. + + for i = 1:length(in) + g = grid.Curvilinear(in{i}{2},in{i}{1}{:}); + testCase.verifyEqual(size(g.coords),out{i}); + end +end + +function testMappingInputComponentMatrix(testCase) + in = { + {{1:3}, [1 2 3]'}, + {{1:2, 1:3}, [1 2 3 4 5 6; 7 8 9 10 11 12]'}, + }; + + for i = 1:length(in) + g = grid.Curvilinear(in{i}{2},in{i}{1}{:}); + testCase.verifyEqual(g.coords,in{i}{2}); + end +end + +function testMappingInputCellOfMatrix(testCase) + + in = { + {{1:3}, {[1 2 3]'}}, + {{1:2, 1:3}, {[1 2 3; 4 5 6], [7 8 9; 10 11 12]}}, + }; + + out = { + [1 2 3]', + [1 2 3 4 5 6; 7 8 9 10 11 12]', + }; + + for i = 1:length(in) + g = grid.Curvilinear(in{i}{2},in{i}{1}{:}); + testCase.verifyEqual(g.coords,out{i}); + end +end + +function testMappingInputCellOfVectors(testCase) + in = { + {{1:3}, {[1 2 3]'}}, + {{1:2, 1:3}, {[1 2 3 4 5 6]', [7 8 9 10 11 12]'}}, + }; + + out = { + [1 2 3]', + [1 2 3 4 5 6; 7 8 9 10 11 12]', + }; +end + +function testMappingInputError(testCase) + testCase.verifyFail(); +end + +function testScaling(testCase) + in = {{1:2, 1:3}, {[1 2 3 4 5 6]', [7 8 9 10 11 12]'}}; + g = grid.Curvilinear(in{2},in{1}{:}); + + testCase.verifyError(@()g.scaling(),'grid:Curvilinear:NoScalingSet'); + + g.logicalGrid.h = [2 1]; + testCase.verifyEqual(g.scaling(),[2 1]); +end + +function testGetBoundaryNames(testCase) + in = { + {{1:10}, @(x) exp(x)}, + {{1:10,1:6}, @(x,y) [exp(x+y); exp(x-y)]}, + {{1:10,1:5,1:7}, @(x,y,z)[exp(x+y+z); exp(x-y-z); 2+x+y-z]}, + }; + + out = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + for i = 1:length(in) + g = grid.Curvilinear(in{i}{2},in{i}{1}{:}); + testCase.verifyEqual(g.getBoundaryNames(), out{i}); + end +end + +function testGetBoundary(testCase) + grids = { + {{1:10}, @(x) exp(x)}, + {{1:10,1:6}, @(x,y) [exp(x+y); exp(x-y)]}, + {{1:10,1:5,1:7}, @(x,y,z)[exp(x+y+z); exp(x-y-z); 2+x+y-z]}, + }; + + boundaries = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + + for ig = 1:length(grids) + g = grid.Curvilinear(grids{ig}{2},grids{ig}{1}{:}); + + logicalGrid = grid.Cartesian(grids{ig}{1}{:}); + + for ib = 1:length(boundaries{ig}) + + logicalBoundary = logicalGrid.getBoundary(boundaries{ig}{ib}); + + x = num2cell(logicalBoundary',2); + expectedBoundary = grids{ig}{2}(x{:})'; + testCase.verifyEqual(g.getBoundary(boundaries{ig}{ib}), expectedBoundary); + end + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Grid.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,25 @@ +classdef Grid < handle + % General d dimensional grid with n points + methods (Abstract) + % n returns the number of points in the grid + o = N(obj) + + % d returns the spatial dimension of the grid + o = D(obj) + + % points returns a n x d matrix containing the coordinates for all points. + X = points(obj) + + % Restricts the grid function gf on obj to the subgrid g. + gf = restrictFunc(obj, gf, g) + + % Projects the grid function gf on obj to the grid g. + gf = projectFunc(obj, gf, g) + + % Return the names of all boundaries in this grid. + bs = getBoundaryNames(obj) + + % Return coordinates for the given boundary + b = getBoundary(obj, name) + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Mapped.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,10 @@ +classdef Mapped < grid.Grid + % General grid mapping + methods (Abstract) + % logicalGrid returns the domain grid of the mapping. + g = logicalGrid(obj); + + % mapping returns the mapped coordinates as a N x D component matrix + m = mapping(obj); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Structured.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,7 @@ +classdef Structured < grid.Grid + methods (Abstract) + % Returns the size of the grid in each dimension m = [mx my mz ...] + m = size(obj); % Is this a good idea? Isn't immersed a structured grid? + h = scaling(obj); + end +end
--- a/+grid/Ti.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,201 +0,0 @@ -classdef Ti - properties - gs % {4}Curve - S % FunctionHandle(u,v) - end - - methods - % TODO function to label boundary names. - % function to find largest and smallest delta h in the grid. Maybe shouldnt live here - function obj = Ti(C1,C2,C3,C4) - obj.gs = {C1,C2,C3,C4}; - - g1 = C1.g; - g2 = C2.g; - g3 = C3.g; - g4 = C4.g; - - A = g1(0); - B = g2(0); - C = g3(0); - D = g4(0); - - function o = S_fun(u,v) - x1 = g1(u); - x2 = g2(v); - x3 = g3(1-u); - x4 = g4(1-v); - o1 = (1-v).*x1(1,:) + u.*x2(1,:) + v.*x3(1,:) + (1-u).*x4(1,:) ... - -((1-u)*(1-v).*A(1,:) + u*(1-v).*B(1,:) + u*v.*C(1,:) + (1-u)*v.*D(1,:)); - o2 = (1-v).*x1(2,:) + u.*x2(2,:) + v.*x3(2,:) + (1-u).*x4(2,:) ... - -((1-u)*(1-v).*A(2,:) + u*(1-v).*B(2,:) + u*v.*C(2,:) + (1-u)*v.*D(2,:)); - - o = [o1;o2]; - end - - obj.S = @S_fun; - end - - function [X,Y] = map(obj,u,v) - default_arg('v',u); - - if isscalar(u) - u = linspace(0,1,u); - end - - if isscalar(v) - v = linspace(0,1,v); - end - - S = obj.S; - - nu = length(u); - nv = length(v); - - X = zeros(nv,nu); - Y = zeros(nv,nu); - - u = rowVector(u); - v = rowVector(v); - - for i = 1:nv - p = S(u,v(i)); - X(i,:) = p(1,:); - Y(i,:) = p(2,:); - end - end - - function h = plot(obj,nu,nv) - S = obj.S; - - default_arg('nv',nu) - - u = linspace(0,1,nu); - v = linspace(0,1,nv); - - m = 100; - - X = zeros(nu+nv,m); - Y = zeros(nu+nv,m); - - - t = linspace(0,1,m); - for i = 1:nu - p = S(u(i),t); - X(i,:) = p(1,:); - Y(i,:) = p(2,:); - end - - for i = 1:nv - p = S(t,v(i)); - X(i+nu,:) = p(1,:); - Y(i+nu,:) = p(2,:); - end - - h = line(X',Y'); - end - - - function h = show(obj,nu,nv) - default_arg('nv',nu) - S = obj.S; - - if(nu>2 || nv>2) - h_grid = obj.plot(nu,nv); - set(h_grid,'Color',[0 0.4470 0.7410]); - end - - h_bord = obj.plot(2,2); - set(h_bord,'Color',[0.8500 0.3250 0.0980]); - set(h_bord,'LineWidth',2); - end - - - % TRANSFORMATIONS - function ti = translate(obj,a) - gs = obj.gs; - - for i = 1:length(gs) - new_gs{i} = gs{i}.translate(a); - end - - ti = grid.Ti(new_gs{:}); - end - - % Mirrors the Ti so that the resulting Ti is still left handed. - % (Corrected by reversing curves and switching e and w) - function ti = mirror(obj, a, b) - gs = obj.gs; - - new_gs = cell(1,4); - - new_gs{1} = gs{1}.mirror(a,b).reverse(); - new_gs{3} = gs{3}.mirror(a,b).reverse(); - new_gs{2} = gs{4}.mirror(a,b).reverse(); - new_gs{4} = gs{2}.mirror(a,b).reverse(); - - ti = grid.Ti(new_gs{:}); - end - - function ti = rotate(obj,a,rad) - gs = obj.gs; - - for i = 1:length(gs) - new_gs{i} = gs{i}.rotate(a,rad); - end - - ti = grid.Ti(new_gs{:}); - end - - function ti = rotate_edges(obj,n); - new_gs = cell(1,4); - for i = 0:3 - new_i = mod(i - n,4); - new_gs{new_i+1} = obj.gs{i+1}; - end - ti = grid.Ti(new_gs{:}); - end - end - - methods(Static) - function obj = points(p1, p2, p3, p4) - g1 = grid.Curve.line(p1,p2); - g2 = grid.Curve.line(p2,p3); - g3 = grid.Curve.line(p3,p4); - g4 = grid.Curve.line(p4,p1); - - obj = grid.Ti(g1,g2,g3,g4); - end - - function label(varargin) - if nargin == 2 && ischar(varargin{2}) - label_impl(varargin{:}); - else - for i = 1:length(varargin) - label_impl(varargin{i},inputname(i)); - end - end - - - function label_impl(ti,str) - S = ti.S; - - pc = S(0.5,0.5); - - margin = 0.1; - pw = S( margin, 0.5); - pe = S(1-margin, 0.5); - ps = S( 0.5, margin); - pn = S( 0.5, 1-margin); - - - ti.show(2,2); - grid.place_label(pc,str); - grid.place_label(pw,'w'); - grid.place_label(pe,'e'); - grid.place_label(ps,'s'); - grid.place_label(pn,'n'); - end - end - end -end \ No newline at end of file
--- a/+grid/equal_step_size.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -% Calculates M so that the stepsize m/M is as close to n/M as possible -function M = equal_step_size(n,N,m) - M = round((m*(N-1)+n)/n); -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equidistant.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,29 @@ +% Creates a cartesian grid of dimension length(m). +% over the doman xlim, ylim, ... +% Examples: +% g = grid.equidistant([mx, my], xlim, ylim) +% g = grid.equidistant([10, 15], {0,1}, {0,2}) +function g = equidistant(m, varargin) + if length(m) ~= length(varargin) + error('grid:equidistant:NonMatchingParameters','The number of provided dimensions do not match.') + end + + for i = 1:length(m) + if ~iscell(varargin{i}) || numel(varargin{i}) ~= 2 + error('grid:equidistant:InvalidLimits','The limits should be cell arrays with 2 elements.'); + end + + if varargin{i}{1} > varargin{i}{2} + error('grid:equidistant:InvalidLimits','The elements of the limit must be increasing.'); + end + end + + X = {}; + h = []; + for i = 1:length(m) + [X{i}, h(i)] = util.get_grid(varargin{i}{:},m(i)); + end + + g = grid.Cartesian(X{:}); + g.h = h; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equidistantCurvilinear.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,35 @@ +% Creates a curvilinear grid of dimension length(m). +% over the logical domain xi_lim, eta_lim, ... +% If all limits are ommited they are set to {0,1}. +% Examples: +% g = grid.equidistantCurvilinear(mapping, [m_xi, m_eta]) +% g = grid.equidistantCurvilinear(mapping, [m_xi, m_eta], xi_lim, eta_lim) +% g = grid.equidistantCurvilinear(mapping, [10, 15], {0,1}, {0,1}) +function g = equidistantCurvilinear(mapping, m, varargin) + if isempty(varargin) + varargin = repmat({{0,1}}, [1 length(m)]); + end + + if length(m) ~= length(varargin) + error('grid:equidistant:NonMatchingParameters','The number of provided dimensions do not match.') + end + + for i = 1:length(m) + if ~iscell(varargin{i}) || numel(varargin{i}) ~= 2 + error('grid:equidistant:InvalidLimits','The limits should be cell arrays with 2 elements.'); + end + + if varargin{i}{1} > varargin{i}{2} + error('grid:equidistant:InvalidLimits','The elements of the limit must be increasing.'); + end + end + + X = {}; + h = []; + for i = 1:length(m) + [X{i}, h(i)] = util.get_grid(varargin{i}{:},m(i)); + end + + g = grid.Curvilinear(mapping, X{:}); + g.h = h; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equidistantCurvilinearTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,7 @@ +function tests = equdistantCurvilinearTest() + tests = functiontests(localfunctions); +end + +function testNoTests(testCase) + testCase.verifyFail(); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equidistantTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,50 @@ +function tests = equidistantTest() + tests = functiontests(localfunctions); +end + + +function testErrorInvalidLimits(testCase) + in = { + {10,{1}}, + {10,[0,1]}, + {[10,10],{0,1},{1}}, + {[10,10],{1},{1,0}}, + {10,{1,0}}, + {[10, 5],{1,0}, {0,-1}}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.equidistant(in{i}{:}),'grid:equidistant:InvalidLimits',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorNonMatchingParam(testCase) + in = { + {[],{1}}, + {[],{1},{0,1}}, + {[5,5],{0,1},{0,1},{0,1}}, + {[5,5,4],{0,1},{0,1}}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.equidistant(in{i}{:}),'grid:equidistant:NonMatchingParameters',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + + +function testCompiles(testCase) + in = { + {5, {0,1}}, + {[3 3],{0,1},{0,2}}, + }; + + out = { + [[0; 0.25; 0.5; 0.75; 1]], + [[0; 0; 0; 0.5; 0.5; 0.5; 1; 1; 1;],[0; 1; 2; 0; 1; 2; 0; 1; 2;]], + }; + + for i = 1:length(in) + g = grid.equidistant(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/evalOn.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,39 @@ +% Takes a funciton and evaluates it on a grid to return a grid function in the +% form of a (n*k)x1 vector, where n is the number of grid points and k is the +% number of components of the function. +% g -- Grid to evaluate on. +% func -- Function to evaluate. May be a function handle or a constant. If +% it is a vector value it has to be provided as a column vector, +function gf = evalOn(g, func) + if ~isa(func, 'function_handle') + % We should have a constant. + if size(func,2) ~= 1 + error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') + end + + gf = repmat(func,[g.N, 1]); + return + end + % func should now be a function_handle + + % Get coordinates and convert to cell array for easier use as a parameter + x = g.points(); + X = {}; + for i = 1:size(x, 2) + X{i} = x(:,i); + end + + % Find the number of components + x0 = num2cell(x(1,:)); + f0 = func(x0{:}); + k = length(f0); + + if size(f0,2) ~= 1 + error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') + end + + gf = func(X{:}); + if k > 1 % Reorder so that componets sits together. + gf = reshape(reshape(gf, [g.N, k])', [g.N*k, 1]); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/evalOnTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,118 @@ +function tests = evalOnTest() + tests = functiontests(localfunctions); +end + +function testInputConstant(testCase) + in = { + 0, + 47, + 1, + [1; 2], + }; + + out = { + [0; 0; 0], + [47; 47; 47], + [1; 1; 1], + [1; 2; 1; 2; 1; 2], + }; + + g = getTestGrid('1d'); + + for i = 1:length(in) + gf = grid.evalOn(g,in{i}); + testCase.verifyEqual(gf, out{i}); + end +end + +function testInputScalarFunction1d(testCase) + in = { + @(x)1+x*0, + @(x)x, + @(x)x.*x, + }; + + out = { + [1; 1; 1], + [0; 1; 2], + [0; 1; 4], + }; + + g = getTestGrid('1d'); + + for i = 1:length(in) + gf = grid.evalOn(g,in{i}); + testCase.verifyEqual(gf, out{i}); + end +end + +function testInputScalarFunction2d(testCase) + in = { + @(x,y)1+x*0, + @(x,y)x-y, + @(x,y)x./(1+y), + }; + + out = { + [1; 1; 1; 1; 1; 1; 1; 1; 1], + [0; -1; -2; 1; 0; -1; 2; 1; 0], + [0; 0; 0; 1; 1/2; 1/3; 2; 1; 2/3], + }; + + g = getTestGrid('2d'); + + for i = 1:length(in) + gf = grid.evalOn(g, in{i}); + testCase.verifyEqual(gf, out{i}); + end +end + + +function testInputVectorFunction(testCase) + g = getTestGrid('1d'); + in = @(x)[x; -2*x]; + out = [0; 0; 1; -2; 2; -4]; + + gf = grid.evalOn(g,in); + testCase.verifyEqual(gf, out); + + g = getTestGrid('2d'); + in = @(x,y)[x.^2; -2*y]; + out = [ + 0; 0; + 0; -2; + 0; -4; + 1; 0; + 1; -2; + 1; -4; + 4; 0; + 4; -2; + 4; -4; + ]; + + gf = grid.evalOn(g,in); + testCase.verifyEqual(gf, out); +end + + +function testInputErrorVectorValued(testCase) + in = { + [1,2,3], + @(x,y)[x,-y]; + }; + + g = getTestGrid('2d'); + + for i = 1:length(in) + testCase.verifyError(@()grid.evalOn(g, in{i}),'grid:evalOn:VectorValuedWrongDim',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function g = getTestGrid(d) + switch d + case '1d' + g = grid.equidistant(3,{0,2}); + case '2d' + g = grid.equidistant([3,3],{0,2},{0,2}); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/funcToMatrix.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,5 @@ +% Converts a gridfunction to a matrix +% Takes a grid function and and a structured grid. +function F = funcToMatrix(g, gf) + reshapeRowMaj(gf, g.size()); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/funcToPlotMatrix.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,5 @@ +% Converts a gridfunction to a plot matrix +% Takes a grid function and and a structured grid. +function F = funcToPlotMatrix(g, gf) + reshapeToPlotMatrix(gf, g.size()); +end \ No newline at end of file
--- a/+grid/old/concat_curve.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% Concatenate two curves g1 and g2 intop -% g = concat_curve(g1,g2) -function g = concat_curve(g1,g2) - function v = g_fun(t) - if t < 1/2 - v = g1(2*t); - else - v = g2(2*t-1); - end - end - g = @g_fun; -end \ No newline at end of file
--- a/+grid/old/curve_discretise.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -% Discretises the curve g with the smallest number of points such that all segments -% are shorter than h. If do_plot is true the points of the discretisation and -% the normals of the curve in those points are plotted. -% -% [t,p,d] = curve_discretise(g,h,do_plot) -% -% t is a vector of input values to g. -% p is a cector of points. -% d are the length of the segments. -function [t,p,d] = curve_discretise(g,h,do_plot) - default_arg('do_plot',false) - - n = 10; - - [t,p,d] = curve_discretise_n(g,n); - - % ni = 0; - while any(d>h) - [t,p,d] = curve_discretise_n(g,n); - n = ceil(n*d(1)/h); - % ni = ni+1; - end - - % nj = 0; - while all(d<h) - [t,p,d] = curve_discretise_n(g,n); - n = n-1; - % nj = nj+1; - end - [t,p,d] = curve_discretise_n(g,n+1); - - % fprintf('ni = %d, nj = %d\n',ni,nj); - - if do_plot - fprintf('n:%d max: %f min: %f\n', n, max(d),min(d)); - p = grid.map_curve(g,t); - figure - show(g,t,h); - end - -end - -function [t,p,d] = curve_discretise_n(g,n) - t = linspace(0,1,n); - t = equalize_d(g,t); - d = D(g,t); - p = grid.map_curve(g,t); -end - -function d = D(g,t) - p = grid.map_curve(g,t); - - d = zeros(1,length(t)-1); - for i = 1:length(d) - d(i) = norm(p(:,i) - p(:,i+1)); - end -end - -function t = equalize_d(g,t) - d = D(g,t); - v = d-mean(d); - while any(abs(v)>0.01*mean(d)) - dt = t(2:end)-t(1:end-1); - t(2:end) = t(2:end) - cumsum(dt.*v./d); - - t = t/t(end); - d = D(g,t); - v = d-mean(d); - end -end - - -function show(g,t,hh) - p = grid.map_curve(g,t); - - - - h = grid.plot_curve(g); - h.LineWidth = 2; - axis equal - hold on - h = plot(p(1,:),p(2,:),'.'); - h.Color = [0.8500 0.3250 0.0980]; - h.MarkerSize = 24; - hold off - - n = grid.curve_normals(g,t); - hold on - for i = 1:length(t) - p0 = p(:,i); - p1 = p0 + hh*n(:,i); - l = [p0, p1]; - h = plot(l(1,:),l(2,:)); - h.Color = [0.8500 0.3250 0.0980]; - end - -end \ No newline at end of file
--- a/+grid/old/curve_interp.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -% Create a cubic spline from the points (x,y) using periodic conditions. -% g = curve_interp(x,y) -function g = curve_interp(x,y) - default_arg('x',[0 2 2 1 1 0]) - default_arg('y',[0 0 2 2 1 1]) - % solve for xp and yp - - % x(t) = at^4 + bt^2+ct+d - - % a = xp1 -2x1 + 2x0 + xp0 - % b = 3x1 -xp1 - 3x0 + 2xp0 - % c = xp0 - % d = x0 - - assert(length(x) == length(y)) - n = length(x); - A = spdiags(ones(n,1)*[2, 8, 2],-1:1,n,n); - A(n,1) = 2; - A(1,n) = 2; - - bx = zeros(n,1); - for i = 2:n-1 - bx(i) = -6*x(i-1)+6*x(i+1); - end - bx(1) = -6*x(n)+6*x(2); - bx(n) = -6*x(n-1)+6*x(1); - - by = zeros(n,1); - for i = 2:n-1 - by(i) = -6*y(i-1)+6*y(i+1); - end - by(1) = -6*y(n)+6*y(2); - by(n) = -6*y(n-1)+6*y(1); - - - xp = A\bx; - yp = A\by; - - x(end+1) = x(1); - y(end+1) = y(1); - - xp(end+1) = xp(1); - yp(end+1) = yp(1); - - function v = g_fun(t) - t = mod(t,1); - i = mod(floor(t*n),n) + 1; - t = t * n -(i-1); - X = (2*x(i)-2*x(i+1)+xp(i)+xp(i+1))*t.^3 + (-3*x(i)+3*x(i+1)-2*xp(i)-xp(i+1))*t.^2 + (xp(i))*t + x(i); - Y = (2*y(i)-2*y(i+1)+yp(i)+yp(i+1))*t.^3 + (-3*y(i)+3*y(i+1)-2*yp(i)-yp(i+1))*t.^2 + (yp(i))*t + y(i); - v = [X;Y]; - end - - g = @g_fun; -end - - -
--- a/+grid/old/max_h.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -function [d_max, i1_max, j1_max, i2_max, j2_max] = max_h(X,Y) - ni = size(X,1); - nj = size(X,2); - d_max = 0; - - i1_max = 0; - j1_max = 0; - i2_max = 0; - j2_max = 0; - - D = {[0,-1],[1,0],[0,1],[-1,0]}; - - for i = 1:ni - for j = 1:nj - p1 = [X(i,j); Y(i,j)]; - for k = 1:length(D) - i2 = i+D{k}(1); - j2 = j+D{k}(2); - if i2 >= 1 && i2 <= ni && j2 >= 1 && j2 <= nj - p2 = [X(i2,j2); Y(i2,j2)]; - d = norm(p2-p1); - if d > d_max; - d_max = d; - i1_max = i; - j1_max = j; - i2_max = i2; - j2_max = j2; - end - end - end - end - end -end
--- a/+grid/old/min_h.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +0,0 @@ -function [d_min, i1_min, j1_min, i2_min, j2_min] = min_h(X,Y) - ni = size(X,1); - nj = size(X,2); - d_min = norm([X(1,1);Y(1,1)] - [X(ni,nj);Y(ni,nj)]); - - i1_min = 0; - j1_min = 0; - i2_min = 0; - j2_min = 0; - - D = {[-1,-1],[0,-1],[1,-1],[1,0],[1,1],[0,1],[-1,1],[-1,0]}; - % D = {[0,-1],[1,0],[0,1],[-1,0]}; - - for i = 1:ni - for j = 1:nj - p1 = [X(i,j); Y(i,j)]; - for k = 1:length(D) - i2 = i+D{k}(1); - j2 = j+D{k}(2); - if i2 >= 1 && i2 <= ni && j2 >= 1 && j2 <= nj - p2 = [X(i2,j2); Y(i2,j2)]; - d = norm(p2-p1); - if d < d_min; - d_min = d; - i1_min = i; - j1_min = j; - i2_min = i2; - j2_min = j2; - end - end - end - end - end -end
--- a/+grid/old/plot_shape.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% Plot a shape using n points. Returns cell array of plot handles. -% hs = plot_shape(s,n) -function hs = plot_shape(s,n) - default_arg('n',100); - - hs = {}; - hold on - for i = 1:length(s) - hs{end+1} = grid.plot_curve(s{i},n); - end - hold off -end \ No newline at end of file
--- a/+grid/old/shape.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -% Creates a shape from a number of curves. A shape is a cell array of curves. -function s = shape(varargin); - s = varargin; -end \ No newline at end of file
--- a/+grid/old/shape_discretise.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -% Discretises a shape such that points on the curves are no further than h appart. -function p = shape_discretise(s,h) - p = []; - for i = 1:length(s) - [~,pt] = grid.curve_discretise(s{i},h); - p = [p, pt]; - end -end \ No newline at end of file
--- a/+grid/old/shape_linesegments.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -% Converts a shape into a cell array of linesegments shorter than h. -function l = shape_linesegments(s,h) - l = {}; - - for i = 1:length(s) - t = grid.curve_discretise(s{i},h); - l = [l, grid.curve_linesegments(s{i},t)]; - end -end \ No newline at end of file
--- a/+grid/old/triang_interp.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,132 +0,0 @@ -classdef triang_interp - properties - g1, g2 ,g3 % Curves encirling the tirangle in the positive direction. - A,B,C % The corners of the triangle - Sa, Sb, Sc % Mappings from square with different sides collapsed - end - - methods - function o = triang_interp(g1,g2,g3) - o.g1 = g1; - o.g2 = g2; - o.g3 = g3; - o.A = g1(0); - o.B = g2(0); - o.C = g3(0); - o.Sa = grid.triang_interp.square_to_triangle_interp(g2,g3,g1); - o.Sb = grid.triang_interp.square_to_triangle_interp(g3,g1,g2); - o.Sc = grid.triang_interp.square_to_triangle_interp(g1,g2,g3); - end - - - function show(o,N) - % Show the mapped meridians of the triangle. - % Might be used for the barycentric coordinates. - ma = @(t)o.Sa(1/2,1-t); - mb = @(t)o.Sb(1/2,1-t); - mc = @(t)o.Sc(1/2,1-t); - - na = @(t)o.Sa(t,1/2); - - ka = @(t)(o.g1(1-t)+o.g2(t))/2; - - h = grid.plot_curve(ma); - h.Color = Color.blue; - h = grid.plot_curve(mb); - h.Color = Color.blue; - h = grid.plot_curve(mc); - h.Color = Color.blue; - - h = grid.plot_curve(na); - h.Color = Color.red; - - h = grid.plot_curve(ka); - h.Color = Color.red; - - [a(1),a(2)] = ma(1/3); - [b(1),b(2)] = mb(1/3); - [c(1),c(2)] = mc(1/3); - - d = ka(1-1/3); - - - grid.label_pt(a,b,c,d); - - - % t = linspace(0,1,N); - % for i = 1:N - % sa = @(s)o.Sa(s,t(i)); - % sb = @(s)o.Sb(s,t(i)); - % sc = @(s)o.Sc(s,t(i)); - - % h = grid.plot_curve(sa); - % h.Color = Color.blue; - % h = grid.plot_curve(sb); - % h.Color = Color.blue; - % h = grid.plot_curve(sc); - % h.Color = Color.blue; - % end - - h = grid.plot_curve(o.g1); - h.LineWidth = 2; - h.Color = Color.red; - - h = grid.plot_curve(o.g2); - h.LineWidth = 2; - h.Color = Color.red; - - h = grid.plot_curve(o.g3); - h.LineWidth = 2; - h.Color = Color.red; - - end - - - end - - methods(Static) - % Makes a mapping from the unit square to a triangle by collapsing - % one of the sides of the squares to a corner on the triangle - % The collapsed side is mapped to the corner oposite to g1. - % This is done such that for S(s,t), S(s,1) = g1(s) - function S = square_to_triangle_interp(g1,g2,g3) - corner = grid.line_segment(g3(0),g3(0)); - S = grid.transfinite_interp(corner,g3,f(g1),f(g2)) - - % Function to flip a curve - function h = f(g) - h = @(t)g(1-t); - end - end - end - -end - -% % Return a mapping from u.v to x,y of the domain encircled by g1 g2 g3 in the the positive direction. created be using transfinite interpolation. -% function S = triang_interp(g1,g2,g3) -% A = g1(0) -% B = g2(0) -% C = g3(0) - -% function [x,y] = S_fun(u,v) -% w = sqrt((u-1)^2+v^2)/sqrt(2); % Parameter for g3 -% v = v*(1-u-v)*g1(u) + u*(1-u-v)*g2(v) + u*v*g3(w) ... -% +(1-u)*(1-v)*A+u*(1-v)*B + (1-u)*v*C; -% x = v(1); -% y = v(2); -% end -% S = @S_fun; -% end - - - -% function subsref(obj,S) -% if ~all(isnumeric(S.subs{:})) -% error('Only supports calling object with number') -% end -% if numel(S.subs{:}) > 1 -% disp('You''ve called the object with more than one argument'); -% else -% disp(['You called the object with argument = ',num2str(S.subs{:})]); -% end -% end \ No newline at end of file
--- a/+grid/old/triang_interp_pts.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% Creates a transfinite interpolation from connecting the four points wiht straight lines. -function [S, g1, g2, g3] = triang_interp_pts(p1,p2,p3) - if size(p1) ~= [2 1] - error('p1 is strange!'); - end - - g1 = @(t)(p1 + t*(p2-p1)); - g2 = @(t)(p2 + t*(p3-p2)); - g3 = @(t)(p3 + t*(p1-p3)); - - S = grid.triang_interp(g1,g2,g3); -end
--- a/+grid/old/triang_map.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,29 +0,0 @@ -% Creates a grid [X,Y] from the mapping function S at points in vectors u,v -function [X, Y] = traing_map(S,u,v) - error('not done') - if nargin == 2 - v = u; - end - - if isscalar(u) - u = linspace(0,1,u); - end - - if isscalar(v) - v = linspace(0,1,v); - end - - nu = length(u); - nv = length(v); - - X = zeros(nu,nv); - Y = zeros(nu,nv); - - for i = 1:nu - for j = 1:nv - [x,y] = S(u(i),v(j)); - X(i,j) = x; - Y(i,j) = y; - end - end -end
--- a/+grid/old/triang_plot_interp.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,114 +0,0 @@ -% Plots a transfinite interpolation in x,y space using nu and nv curves along u and v axes. - - - - - - -% Plots a interp of a triangle where one the interpolation is from a square -% with one side collapsed to -function h = triang_plot_interp_kindaworking(S,n) - u = linspace(0,1,n); - v = linspace(0,1,n); - - m = 100; - m = 20; - - Xl_curves = cell(n,1); - Xr_curves = cell(n,1); - Y_curves = cell(n,1); - - - function u = wierdness(v,d,N) - if N == 0 - u = 0; - else - u = N*d./(1-v); - end - end - - - %Y curves - t = linspace(0,1,m); - for i = 1:n - x = []; y = []; - for j = 1:length(t) - [x(j),y(j)] = S(t(j),v(i)); - end - Y_curves{i} = [x', y']; - end - - - % Right and left X curves - t = linspace(0,1,m); - d = u(2); - for i = 1:n - xl = []; yl = []; - xr = []; yr = []; - N = i-1; - t = linspace(0,1-N*d,m); - for j = 1:length(t) - w = wierdness(t(j),d,N); - [xr(j),yr(j)] = S(w,t(j)); - [xl(j),yl(j)] = S(1-w,t(j)); - end - Xl_curves{i} = [xl', yl']; - Xr_curves{i} = [xr', yr']; - end - - for i = 1:n-1 - line(Xl_curves{i}(:,1),Xl_curves{i}(:,2)) - line(Xr_curves{i}(:,1),Xr_curves{i}(:,2)) - line(Y_curves{i}(:,1),Y_curves{i}(:,2)) - end -end - - - - -function h = triang_plot_interp_nonworking(S,n) - - u = linspace(0,1,n); - v = linspace(0,1,n); - - m = 100; - - X_curves = cell(n-1,1); - Y_curves = cell(n-1,1); - K_curves = cell(n-1,1); - - - t = linspace(0,1,m); - for i = 1:n-1 - x = []; y = []; - for j = find(t+u(i) <= 1) - [x(j),y(j)] = S(u(i),t(j)); - end - X_curves{i} = [x', y']; - end - - for i = 1:n-1 - x = []; y = []; - for j = find(t+v(i) <= 1) - [x(j),y(j)] = S(t(j),v(i)); - end - Y_curves{i} = [x', y']; - end - - for i = 2:n - x = []; y = []; - for j = find(t<u(i)) - [x(j),y(j)] = S(t(j), u(i)-t(j)); - end - K_curves{i-1} = [x', y']; - end - - for i = 1:n-1 - line(X_curves{i}(:,1),X_curves{i}(:,2)) - line(Y_curves{i}(:,1),Y_curves{i}(:,2)) - line(K_curves{i}(:,1),K_curves{i}(:,2)) - end - - h = -1; - % h = plot(X_curves{:},Y_curves{:},K_curves{:}); -end
--- a/+grid/old/triang_show.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -% Show a grid as a matlab figure. -function triang_show(S,n) - - ih = ishold(); - - hold on - h_grid = grid.triang_plot_interp(S,n); - h_bord = grid.triang_plot_interp(S,2); - - set(h_grid,'Color',[0 0.4470 0.7410]); - set(h_bord,'Color',[0.8500 0.3250 0.0980]); - set(h_bord,'LineWidth',2); - - % axis auto - % axis equal - % axis square - - if ih - hold on - else - hold off - end -end
--- a/+grid/place_label.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% 'left' | 'center' | 'right' -% 'baseline' | 'top' | 'cap' | 'middle' | 'bottom' -function place_label(pt,str,horzAl,vertAl) - default_arg('horzAl','center'); - default_arg('vertAl', 'middle'); - - x = pt(1); - y = pt(2); - h = text(x,y,str); - h.HorizontalAlignment = horzAl; - h.VerticalAlignment = vertAl; -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/BoundaryGroup.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,51 @@ +% BoundaryGroup defines a boundary grouping in a multiblock grid. +classdef BoundaryGroup + properties + blockIDs + names + end + + methods + function obj = BoundaryGroup(varargin) + % Input arguemnts are arbitrary number or 1x2 cell arrays + % representing each boundary in the group. + % The 1st element of the cell array is an integer defining which grid it belongs to. + % The 2nd element of the cell array is the name of the boundary within the block. + % + % Ex: + % bg = multiblock.BoundaryGroup({1,'n'},{1,'s'},{2,'s'}) + + + obj.blockIDs = []; + obj.names = {}; + for i = 1:length(varargin) + if ~iscell(varargin{i}) || ~all(size(varargin{i}) == [1 2]) + error('multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', 'Inputs must be 1x2 cell arrays'); + end + obj.blockIDs(i) = varargin{i}{1}; + obj.names{i} = varargin{i}{2}; + end + end + + function display(obj, name) + + disp(' ') + disp([name, ' =']) + disp(' ') + + if length(obj.names) == 1 + fprintf(' {}\n\n') + return + end + + fprintf(' {') + + fprintf('%d:%s', obj.blockIDs(1), obj.names{1}) + for i = 2:length(obj.names) + fprintf(', %d:%s', obj.blockIDs(i), obj.names{i}); + end + + fprintf('}\n\n') + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/BoundaryGroupTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,30 @@ +function tests = BoundaryGroupTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + in = {{3,'n'},{2,'hoho'},{1,'s'}}; + + blockIDs = [3 2 1]; + names = {'n', 'hoho', 's'}; + + bg = multiblock.BoundaryGroup(in{:}); + testCase.verifyEqual(bg.blockIDs, blockIDs); + testCase.verifyEqual(bg.names, names); +end + +function testInputError(testCase) + in = { + {'n', 's'}, + {{3,'n'},{2,2,'hoho'},{1,'s'}}, + }; + + out = { + 'multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', + 'multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', + }; + + for i = 1:length(in) + testCase.verifyError(@()multiblock.BoundaryGroup(in{i}{:}), out{i}); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOp.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,144 @@ +classdef DiffOp < scheme.Scheme + properties + grid + order + diffOps + D + H + + blockmatrixDiv + end + + methods + function obj = DiffOp(doHand, grid, order, doParam) + % doHand -- may either be a function handle or a cell array of + % function handles for each grid. The function handle(s) + % should be on the form do = doHand(grid, order, ...) + % Additional parameters for each doHand may be provided in + % the doParam input. + % grid -- a multiblock grid + % order -- integer specifying the order of accuracy + % doParam -- may either be a cell array or a cell array of cell arrays + % for each block. If it is a cell array with length equal + % to the number of blocks then each element is sent to the + % corresponding function handle as extra parameters: + % doHand(..., doParam{i}{:}) Otherwise doParam is sent as + % extra parameters to all doHand: doHand(..., doParam{:}) + default_arg('doParam', []) + + [getHand, getParam] = parseInput(doHand, grid, doParam); + + nBlocks = grid.nBlocks(); + + obj.order = order; + + % Create the diffOps for each block + obj.diffOps = cell(1, nBlocks); + for i = 1:nBlocks + h = getHand(i); + p = getParam(i); + obj.diffOps{i} = h(grid.grids{i}, order, p{:}); + end + + + % Build the norm matrix + H = cell(nBlocks, nBlocks); + for i = 1:nBlocks + H{i,i} = obj.diffOps{i}.H; + end + obj.H = blockmatrix.toMatrix(H); + + + % Build the differentiation matrix + D = cell(nBlocks, nBlocks); + for i = 1:nBlocks + D{i,i} = obj.diffOps{i}.D; + end + + for i = 1:nBlocks + for j = i:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + D{i,i} = D{i,i} + ii; + D{i,j} = D{i,j} + ij; + + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = D{j,j} + jj; + D{j,i} = D{j,i} + ji; + end + end + obj.D = blockmatrix.toMatrix(D); + + obj.blockmatrixDiv = blockmatrix.getDivision(D); + + function [getHand, getParam] = parseInput(doHand, grid, doParam) + if ~isa(grid, 'multiblock.Grid') + error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); + end + + if iscell(doHand) && length(doHand) == grid.nBlocks() + getHand = @(i)doHand{i}; + elseif isa(doHand, 'function_handle') + getHand = @(i)doHand; + else + error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); + end + + if isempty(doParam) + getParam = @(i){}; + elseif iscell(doParam) && length(doParam) == grid.nBlocks() + getParam = @(i)doParam{i}; + else + getParam = @(i)doParam; + end + end + end + + function ops = splitOp(obj, op) + % Splits a matrix operator into a cell-matrix of matrix operators for + % each grid. + ops = sparse2cell(op, obj.NNN); + end + + % Creates the closere and penalty matrix for a given boundary condition, + % boundary -- the name of the boundary on the form [id][name] where + % id is the number of a block and name is the name of a + % boundary of that block example: 1s or 3w + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + I = boundary(1) + name = boundary(2:end); + + % Get the closure and penaly matrices + [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type, data); + + % Expand to matrix for full domain. + div = obj.blockmatrixDiv; + closure = blockmatrix.zero(div); + closure{I,I} = blockClosure; + + div{2} = 1; % Penalty is a column vector + penalty = blockmatrix.zero(div); + penalty{I} = blockPenalty; + + % Convert to matrix + closure = blockmatrix.toMatrix(closure); + penalty = blockmatrix.toMatrix(closure); + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + end + + % Size returns the number of degrees of freedom + function N = size(obj) + N = 0; + for i = 1:length(obj.diffOps) + N = N + obj.diffOps{i}.size(); + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOpTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,59 @@ +function tests = DiffOpTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + do = newMultiblockOp(); +end + +function testSplitOp(testCase) + testCase.verifyFail(); +end + +function testBoundary_condition(testCase) + testCase.verifyFail(); +end + +function testInterface(testCase) + testCase.verifyFail(); +end + +function testSize(testCase) + mbDo = newMultiblockOp(); + testCase.verifyEqual(mbDo.size(), 15) +end + + +function do = mockDiffOp(size, bc, interface) + do.H = 1; + do.D = 1; + + do.size = size; + do.boundary_condition = bc; + do.interface = interface; +end + + +function do = newMultiblockOp() + grids = { + grid.Cartesian([0 1 2], [3 4 5]); + grid.Cartesian([1 2 3], [10 20]); + }; + + conn = cell(2,2); + conn{1, 2} = {'s','n'}; + + mbGrid = multiblock.Grid(grids, conn); + + function [c, p] = boundary_condition(~,~,~,~) + c = 1; p = 1; + end + + function [c, p] = interface(~,~,~,~) + c = 1; p = 1; + end + + doHand = @(grid,~)mockDiffOp(@(~)prod(grid.size()), @boundary_condition, @interface); + + do = multiblock.DiffOp(doHand, mbGrid, 0); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Grid.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,111 @@ +classdef Grid < grid.Grid + properties + grids + connections + boundaryGroups + + nPoints + end + + % General multiblock grid + methods + + % grids -- cell array of N grids + % connections -- NxN upper triangular cell matrix. connections{i,j} + % specifies the connection between block i and j. If + % it's empty there is no connection otherwise it's a 2 + % -cell-vector with strings naming the boundaries to be + % connected. (inverted coupling?) + %% Should we have boundary groups at all? maybe it can be handled in a + %% cleaner way outside of the class. + function obj = Grid(grids, connections, boundaryGroups) + obj.grids = grids; + obj.connections = connections; + + obj.nPoints = 0; + for i = 1:length(grids) + obj.nPoints = obj.nPoints + grids{i}.N(); + end + + % if iscell(boundaryGroups) + end + + function n = size(obj) + n = length(obj.grids); + end + + % n returns the number of points in the grid + function o = N(obj) + o = obj.nPoints; + end + + function n = nBlocks(obj) + n = length(obj.grids); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.grids{1}.D(); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = []; + for i = 1:length(obj.grids) + X = [X; obj.grids{i}.points]; + end + end + + % Split a grid function on obj to a cell array of grid function on each block + function gfs = splitFunc(obj, gf) + nComponents = length(gf)/obj.nPoints; + nBlocks = length(obj.grids); + + % Collect number of points in each block + N = cell(1,nBlocks); + for i = 1:nBlocks + N{i} = obj.grids{i}.N(); + end + + gfs = mat2cell(gf, N, 1); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + gfs = obj.splitFunc(gf); + + for i = 1:length(obj.grids) + gfs{i} = obj.grids{i}.restrictFunc(gfs{i}, g.grids{i}); + end + + gf = cell2mat(gfs); + end + + % Projects the grid function gf on obj to the grid g. + function o = projectFunc(obj, gf, g) + error('not implemented') + + p = g.points(); + o = zeros(length(p),1); + for i = 1:length(p) + I = whatGrid(p(i)); + o(i) = obj.grids{I}.projectFunc(gf, p(i)); + end + + + function I = whatGrid(p) + % Find what grid a point lies on + end + + end + + function bs = getBoundaryNames(obj) + bs = []; + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + b = []; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/GridTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,37 @@ +function tests = GridTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + g = multiblock.Grid({},{}); +end + +function testMissing(testCase) + testCase.verifyFail(); +end + +function testGetBoundaryNames(testCase) + [grids, conn] = prepareAdjecentBlocks(); + + mbg = multiblock.Grid(grids, conn, multiblock.BoundaryGroup({1,'w'},{2,'w'}) ); + + testCase.verifyFail(); +end + +function testGetBoundary(testCase) + [grids, conn] = prepareAdjecentBlocks(); + + mbg = multiblock.Grid(grids, conn, multiblock.BoundaryGroup({1,'w'},{2,'w'}) ); + testCase.verifyFail(); +end + + +function [grids, conn] = prepareAdjecentBlocks() + grids = { + grid.Cartesian([0 1 2], [3 4 5]); + grid.Cartesian([1 2], [10 20]); + }; + + conn = cell(2,2); + conn{1, 2} = {'s','n'}; +end \ No newline at end of file
--- a/+multiblock/gridVector1d.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function x = gridVector1d(schms) - n = length(schms); - x = cell(n,1); - - for i = 1:n - x{i} = schms{i}.x; - end -end \ No newline at end of file
--- a/+multiblock/gridVector2d.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function [x,y] = gridVector2d(schms) - n = length(schms); - x = cell(n,1); - y = cell(n,1); - - for i = 1:n - x{i} = schms{i}.x; - y{i} = schms{i}.y; - end -end \ No newline at end of file
--- a/+multiblock/solutionVector2cell.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function v_cell = solutionVector2cell(schms,v,components) - if nargin == 2 - components = 1; - end - - N = length(schms); - v_cell = cell(N,1); - - i_start = 1; - for i =1:N - i_end = i_start + schms{i}.size()*components - 1; - v_cell{i} = v(i_start:i_end); - i_start = i_end+1; - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/Curve.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,388 @@ +classdef Curve + properties + g + gp + transformation + end + + methods + %TODO: + % Concatenation of curves + % Subsections of curves + % Stretching of curve parameter - done for arc length. + % Curve to cell array of linesegments + + % Returns a curve object. + % g -- curve parametrization for parameter between 0 and 1 + % gp -- parametrization of curve derivative + function obj = Curve(g,gp,transformation) + default_arg('gp',[]); + default_arg('transformation',[]); + p_test = g(0); + assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.'); + + if ~isempty(transformation) + transformation.base_g = g; + transformation.base_gp = gp; + [g,gp] = grid.Curve.transform_g(g,gp,transformation); + end + + obj.g = g; + obj.gp = gp; + obj.transformation = transformation; + + end + + function n = normal(obj,t) + assert(~isempty(obj.gp),'Curve has no derivative!'); + deriv = obj.gp(t); + normalization = sqrt(sum(deriv.^2,1)); + n = [-deriv(2,:)./normalization; deriv(1,:)./normalization]; + end + + + % Plots a curve g(t) for 0<t<1, using n points. Returns a handle h to the plotted curve. + % h = plot_curve(g,n) + function h = plot(obj,n) + default_arg('n',100); + + t = linspace(0,1,n); + p = obj.g(t); + h = line(p(1,:),p(2,:)); + end + + function h= plot_normals(obj,l,n,m) + default_arg('l',0.1); + default_arg('n',10); + default_arg('m',100); + t_n = linspace(0,1,n); + + normals = obj.normal(t_n)*l; + + n0 = obj.g(t_n); + n1 = n0 + normals; + + h = line([n0(1,:); n1(1,:)],[n0(2,:); n1(2,:)]); + set(h,'Color',Color.red); + obj.plot(m); + end + + function h= show(obj,name) + p = obj.g(1/2); + n = obj.normal(1/2); + p = p + n*0.1; + + % Add arrow + + h = text(p(1),p(2),name); + h.HorizontalAlignment = 'center'; + h.VerticalAlignment = 'middle'; + + obj.plot(); + end + % Shows curve with name and arrow for direction. + + % Length of arc from parameter t0 to t1 (which may be vectors). + % Computed using derivative. + function L = arcLength(obj,t0,t1) + assert(~isempty(obj.gp),'Curve has no derivative!'); + speed = @(t) sqrt(sum(obj.gp(t).^2,1)); + L = integral_vec(speed,t0,t1); + end + + % Creates the arc length parameterization of a curve. + % N -- number of points used to approximate the arclength function + function curve = arcLengthParametrization(obj,N) + default_arg('N',100); + assert(~isempty(obj.gp),'Curve has no derivative!'); + + % Construct arcLength function using splines + tvec = linspace(0,1,N); + arcVec = obj.arcLength(0,tvec); + tFunc = spline(arcVec,tvec); % t as a function of arcLength + L = obj.arcLength(0,1); + arcPar = @(s) tFunc(s*L); + + % New function and derivative + g_new = @(t)obj.g(arcPar(t)); + gp_new = @(t) normalize(obj.gp(arcPar(t))); + curve = grid.Curve(g_new,gp_new); + + end + + % how to make it work for methods without returns + function p = subsref(obj,S) + %Should i add error checking here? + %Maybe if you want performance you fetch obj.g and then use that + switch S(1).type + case '()' + p = obj.g(S.subs{1}); + % case '.' + + % p = obj.(S.subs); + otherwise + p = builtin('subsref',obj,S); + % error() + end + end + + + %% TRANSFORMATION OF A CURVE + function D = reverse(C) + % g = C.g; + % gp = C.gp; + % D = grid.Curve(@(t)g(1-t),@(t)-gp(1-t)); + D = C.transform([],[],-1); + end + + function D = transform(C,A,b,flip) + default_arg('A',[1 0; 0 1]); + default_arg('b',[0; 0]); + default_arg('flip',1); + if isempty(C.transformation) + g = C.g; + gp = C.gp; + transformation.A = A; + transformation.b = b; + transformation.flip = flip; + else + g = C.transformation.base_g; + gp = C.transformation.base_gp; + A_old = C.transformation.A; + b_old = C.transformation.b; + flip_old = C.transformation.flip; + + transformation.A = A*A_old; + transformation.b = A*b_old + b; + transformation.flip = flip*flip_old; + end + + D = grid.Curve(g,gp,transformation); + + end + + function D = translate(C,a) + g = C.g; + gp = C.gp; + + % function v = g_fun(t) + % x = g(t); + % v(1,:) = x(1,:)+a(1); + % v(2,:) = x(2,:)+a(2); + % end + + % D = grid.Curve(@g_fun,gp); + + D = C.transform([],a); + end + + function D = mirror(C, a, b) + assert_size(a,[2,1]); + assert_size(b,[2,1]); + + g = C.g; + gp = C.gp; + + l = b-a; + lx = l(1); + ly = l(2); + + + % fprintf('Singular?\n') + + A = [lx^2-ly^2 2*lx*ly; 2*lx*ly ly^2-lx^2]/(l'*l); + + % function v = g_fun(t) + % % v = a + A*(g(t)-a) + % x = g(t); + + % ax1 = x(1,:)-a(1); + % ax2 = x(2,:)-a(2); + % v(1,:) = a(1)+A(1,:)*[ax1;ax2]; + % v(2,:) = a(2)+A(2,:)*[ax1;ax2]; + % end + + % function v = gp_fun(t) + % v = A*gp(t); + % end + + % D = grid.Curve(@g_fun,@gp_fun); + + % g = A(g-a)+a = Ag - Aa + a; + b = - A*a + a; + D = C.transform(A,b); + + end + + function D = rotate(C,a,rad) + assert_size(a, [2,1]); + assert_size(rad, [1,1]); + g = C.g; + gp = C.gp; + + + A = [cos(rad) -sin(rad); sin(rad) cos(rad)]; + + % function v = g_fun(t) + % % v = a + A*(g(t)-a) + % x = g(t); + + % ax1 = x(1,:)-a(1); + % ax2 = x(2,:)-a(2); + % v(1,:) = a(1)+A(1,:)*[ax1;ax2]; + % v(2,:) = a(2)+A(2,:)*[ax1;ax2]; + % end + + % function v = gp_fun(t) + % v = A*gp(t); + % end + + % D = grid.Curve(@g_fun,@gp_fun); + + + % g = A(g-a)+a = Ag - Aa + a; + b = - A*a + a; + D = C.transform(A,b); + end + end + + methods (Static) + + % Computes the derivative of g: R -> R^2 using an operator D1 + function gp_out = numericalDerivative(g,D1) + m = length(D1); + t = linspace(0,1,m); + gVec = g(t)'; + gpVec = (D1*gVec)'; + + gp1_fun = spline(t,gpVec(1,:)); + gp2_fun = spline(t,gpVec(2,:)); + gp_out = @(t) [gp1_fun(t);gp2_fun(t)]; + end + + function obj = line(p1, p2) + + function v = g_fun(t) + v(1,:) = p1(1) + t.*(p2(1)-p1(1)); + v(2,:) = p1(2) + t.*(p2(2)-p1(2)); + end + g = @g_fun; + + obj = grid.Curve(g); + end + + function obj = circle(c,r,phi) + default_arg('phi',[0; 2*pi]) + default_arg('c',[0; 0]) + default_arg('r',1) + + function v = g_fun(t) + w = phi(1)+t*(phi(2)-phi(1)); + v(1,:) = c(1) + r*cos(w); + v(2,:) = c(2) + r*sin(w); + end + + function v = g_fun_deriv(t) + w = phi(1)+t*(phi(2)-phi(1)); + v(1,:) = -(phi(2)-phi(1))*r*sin(w); + v(2,:) = (phi(2)-phi(1))*r*cos(w); + end + + obj = grid.Curve(@g_fun,@g_fun_deriv); + end + + function obj = bezier(p0, p1, p2, p3) + function v = g_fun(t) + v(1,:) = (1-t).^3*p0(1) + 3*(1-t).^2.*t*p1(1) + 3*(1-t).*t.^2*p2(1) + t.^3*p3(1); + v(2,:) = (1-t).^3*p0(2) + 3*(1-t).^2.*t*p1(2) + 3*(1-t).*t.^2*p2(2) + t.^3*p3(2); + end + + function v = g_fun_deriv(t) + v(1,:) = 3*(1-t).^2*(p1(1)-p0(1)) + 6*(1-t).*t*(p2(1)-p1(1)) + 3*t.^2*(p3(1)-p2(1)); + v(2,:) = 3*(1-t).^2*(p1(2)-p0(2)) + 6*(1-t).*t*(p2(2)-p1(2)) + 3*t.^2*(p3(2)-p2(2)); + end + + obj = grid.Curve(@g_fun,@g_fun_deriv); + end + + + function [g_out,gp_out] = transform_g(g,gp,tr) + A = tr.A; + b = tr.b; + flip = tr.flip; + + function v = g_fun_noflip(t) + % v = A*g + b + x = g(t); + + v(1,:) = A(1,:)*x+b(1); + v(2,:) = A(2,:)*x+b(2); + end + + function v = g_fun_flip(t) + % v = A*g + b + x = g(1-t); + + v(1,:) = A(1,:)*x+b(1); + v(2,:) = A(2,:)*x+b(2); + end + + + switch flip + case 1 + g_out = @g_fun_noflip; + gp_out = @(t)A*gp(t); + case -1 + g_out = @g_fun_flip; + gp_out = @(t)-A*gp(1-t); + end + end + + end +end + + + +function g_norm = normalize(g0) + g1 = g0(1,:); + g2 = g0(2,:); + normalization = sqrt(sum(g0.^2,1)); + g_norm = [g1./normalization; g2./normalization]; +end + +function I = integral_vec(f,a,b) +% Wrapper around the built-in function integral that +% handles multiple limits. + + Na = length(a); + Nb = length(b); + assert(Na == 1 || Nb == 1 || Na==Nb,... + 'a and b must have same length, unless one is a scalar.'); + + if(Na>Nb); + I = zeros(size(a)); + for i = 1:Na + I(i) = integral(f,a(i),b); + end + elseif(Nb>Na) + I = zeros(size(b)); + for i = 1:Nb + I(i) = integral(f,a,b(i)); + end + else + I = zeros(size(b)); + for i = 1:Nb + I(i) = integral(f,a(i),b(i)); + end + end +end + +% Returns a function handle to the spline. +function f = spline(tval,fval,spline_order) + default_arg('spline_order',4); + [m,~] = size(tval); + assert(m==1,'Need row vectors.'); + + f_spline = spapi( optknt(tval,spline_order), tval, fval ); + f = @(t) fnval(f_spline,t); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/Ti.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,201 @@ +classdef Ti + properties + gs % {4}Curve + S % FunctionHandle(u,v) + end + + methods + % TODO function to label boundary names. + % function to find largest and smallest delta h in the grid. Maybe shouldnt live here + function obj = Ti(C1,C2,C3,C4) + obj.gs = {C1,C2,C3,C4}; + + g1 = C1.g; + g2 = C2.g; + g3 = C3.g; + g4 = C4.g; + + A = g1(0); + B = g2(0); + C = g3(0); + D = g4(0); + + function o = S_fun(u,v) + x1 = g1(u); + x2 = g2(v); + x3 = g3(1-u); + x4 = g4(1-v); + o1 = (1-v).*x1(1,:) + u.*x2(1,:) + v.*x3(1,:) + (1-u).*x4(1,:) ... + -((1-u)*(1-v).*A(1,:) + u*(1-v).*B(1,:) + u*v.*C(1,:) + (1-u)*v.*D(1,:)); + o2 = (1-v).*x1(2,:) + u.*x2(2,:) + v.*x3(2,:) + (1-u).*x4(2,:) ... + -((1-u)*(1-v).*A(2,:) + u*(1-v).*B(2,:) + u*v.*C(2,:) + (1-u)*v.*D(2,:)); + + o = [o1;o2]; + end + + obj.S = @S_fun; + end + + function [X,Y] = map(obj,u,v) + default_arg('v',u); + + if isscalar(u) + u = linspace(0,1,u); + end + + if isscalar(v) + v = linspace(0,1,v); + end + + S = obj.S; + + nu = length(u); + nv = length(v); + + X = zeros(nv,nu); + Y = zeros(nv,nu); + + u = rowVector(u); + v = rowVector(v); + + for i = 1:nv + p = S(u,v(i)); + X(i,:) = p(1,:); + Y(i,:) = p(2,:); + end + end + + function h = plot(obj,nu,nv) + S = obj.S; + + default_arg('nv',nu) + + u = linspace(0,1,nu); + v = linspace(0,1,nv); + + m = 100; + + X = zeros(nu+nv,m); + Y = zeros(nu+nv,m); + + + t = linspace(0,1,m); + for i = 1:nu + p = S(u(i),t); + X(i,:) = p(1,:); + Y(i,:) = p(2,:); + end + + for i = 1:nv + p = S(t,v(i)); + X(i+nu,:) = p(1,:); + Y(i+nu,:) = p(2,:); + end + + h = line(X',Y'); + end + + + function h = show(obj,nu,nv) + default_arg('nv',nu) + S = obj.S; + + if(nu>2 || nv>2) + h_grid = obj.plot(nu,nv); + set(h_grid,'Color',[0 0.4470 0.7410]); + end + + h_bord = obj.plot(2,2); + set(h_bord,'Color',[0.8500 0.3250 0.0980]); + set(h_bord,'LineWidth',2); + end + + + % TRANSFORMATIONS + function ti = translate(obj,a) + gs = obj.gs; + + for i = 1:length(gs) + new_gs{i} = gs{i}.translate(a); + end + + ti = grid.Ti(new_gs{:}); + end + + % Mirrors the Ti so that the resulting Ti is still left handed. + % (Corrected by reversing curves and switching e and w) + function ti = mirror(obj, a, b) + gs = obj.gs; + + new_gs = cell(1,4); + + new_gs{1} = gs{1}.mirror(a,b).reverse(); + new_gs{3} = gs{3}.mirror(a,b).reverse(); + new_gs{2} = gs{4}.mirror(a,b).reverse(); + new_gs{4} = gs{2}.mirror(a,b).reverse(); + + ti = grid.Ti(new_gs{:}); + end + + function ti = rotate(obj,a,rad) + gs = obj.gs; + + for i = 1:length(gs) + new_gs{i} = gs{i}.rotate(a,rad); + end + + ti = grid.Ti(new_gs{:}); + end + + function ti = rotate_edges(obj,n); + new_gs = cell(1,4); + for i = 0:3 + new_i = mod(i - n,4); + new_gs{new_i+1} = obj.gs{i+1}; + end + ti = grid.Ti(new_gs{:}); + end + end + + methods(Static) + function obj = points(p1, p2, p3, p4) + g1 = grid.Curve.line(p1,p2); + g2 = grid.Curve.line(p2,p3); + g3 = grid.Curve.line(p3,p4); + g4 = grid.Curve.line(p4,p1); + + obj = grid.Ti(g1,g2,g3,g4); + end + + function label(varargin) + if nargin == 2 && ischar(varargin{2}) + label_impl(varargin{:}); + else + for i = 1:length(varargin) + label_impl(varargin{i},inputname(i)); + end + end + + + function label_impl(ti,str) + S = ti.S; + + pc = S(0.5,0.5); + + margin = 0.1; + pw = S( margin, 0.5); + pe = S(1-margin, 0.5); + ps = S( 0.5, margin); + pn = S( 0.5, 1-margin); + + + ti.show(2,2); + grid.place_label(pc,str); + grid.place_label(pw,'w'); + grid.place_label(pe,'e'); + grid.place_label(ps,'s'); + grid.place_label(pn,'n'); + end + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/equal_step_size.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,4 @@ +% Calculates M so that the stepsize m/M is as close to n/M as possible +function M = equal_step_size(n,N,m) + M = round((m*(N-1)+n)/n); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/concat_curve.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,12 @@ +% Concatenate two curves g1 and g2 intop +% g = concat_curve(g1,g2) +function g = concat_curve(g1,g2) + function v = g_fun(t) + if t < 1/2 + v = g1(2*t); + else + v = g2(2*t-1); + end + end + g = @g_fun; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/curve_discretise.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,97 @@ +% Discretises the curve g with the smallest number of points such that all segments +% are shorter than h. If do_plot is true the points of the discretisation and +% the normals of the curve in those points are plotted. +% +% [t,p,d] = curve_discretise(g,h,do_plot) +% +% t is a vector of input values to g. +% p is a cector of points. +% d are the length of the segments. +function [t,p,d] = curve_discretise(g,h,do_plot) + default_arg('do_plot',false) + + n = 10; + + [t,p,d] = curve_discretise_n(g,n); + + % ni = 0; + while any(d>h) + [t,p,d] = curve_discretise_n(g,n); + n = ceil(n*d(1)/h); + % ni = ni+1; + end + + % nj = 0; + while all(d<h) + [t,p,d] = curve_discretise_n(g,n); + n = n-1; + % nj = nj+1; + end + [t,p,d] = curve_discretise_n(g,n+1); + + % fprintf('ni = %d, nj = %d\n',ni,nj); + + if do_plot + fprintf('n:%d max: %f min: %f\n', n, max(d),min(d)); + p = grid.map_curve(g,t); + figure + show(g,t,h); + end + +end + +function [t,p,d] = curve_discretise_n(g,n) + t = linspace(0,1,n); + t = equalize_d(g,t); + d = D(g,t); + p = grid.map_curve(g,t); +end + +function d = D(g,t) + p = grid.map_curve(g,t); + + d = zeros(1,length(t)-1); + for i = 1:length(d) + d(i) = norm(p(:,i) - p(:,i+1)); + end +end + +function t = equalize_d(g,t) + d = D(g,t); + v = d-mean(d); + while any(abs(v)>0.01*mean(d)) + dt = t(2:end)-t(1:end-1); + t(2:end) = t(2:end) - cumsum(dt.*v./d); + + t = t/t(end); + d = D(g,t); + v = d-mean(d); + end +end + + +function show(g,t,hh) + p = grid.map_curve(g,t); + + + + h = grid.plot_curve(g); + h.LineWidth = 2; + axis equal + hold on + h = plot(p(1,:),p(2,:),'.'); + h.Color = [0.8500 0.3250 0.0980]; + h.MarkerSize = 24; + hold off + + n = grid.curve_normals(g,t); + hold on + for i = 1:length(t) + p0 = p(:,i); + p1 = p0 + hh*n(:,i); + l = [p0, p1]; + h = plot(l(1,:),l(2,:)); + h.Color = [0.8500 0.3250 0.0980]; + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/curve_interp.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,58 @@ +% Create a cubic spline from the points (x,y) using periodic conditions. +% g = curve_interp(x,y) +function g = curve_interp(x,y) + default_arg('x',[0 2 2 1 1 0]) + default_arg('y',[0 0 2 2 1 1]) + % solve for xp and yp + + % x(t) = at^4 + bt^2+ct+d + + % a = xp1 -2x1 + 2x0 + xp0 + % b = 3x1 -xp1 - 3x0 + 2xp0 + % c = xp0 + % d = x0 + + assert(length(x) == length(y)) + n = length(x); + A = spdiags(ones(n,1)*[2, 8, 2],-1:1,n,n); + A(n,1) = 2; + A(1,n) = 2; + + bx = zeros(n,1); + for i = 2:n-1 + bx(i) = -6*x(i-1)+6*x(i+1); + end + bx(1) = -6*x(n)+6*x(2); + bx(n) = -6*x(n-1)+6*x(1); + + by = zeros(n,1); + for i = 2:n-1 + by(i) = -6*y(i-1)+6*y(i+1); + end + by(1) = -6*y(n)+6*y(2); + by(n) = -6*y(n-1)+6*y(1); + + + xp = A\bx; + yp = A\by; + + x(end+1) = x(1); + y(end+1) = y(1); + + xp(end+1) = xp(1); + yp(end+1) = yp(1); + + function v = g_fun(t) + t = mod(t,1); + i = mod(floor(t*n),n) + 1; + t = t * n -(i-1); + X = (2*x(i)-2*x(i+1)+xp(i)+xp(i+1))*t.^3 + (-3*x(i)+3*x(i+1)-2*xp(i)-xp(i+1))*t.^2 + (xp(i))*t + x(i); + Y = (2*y(i)-2*y(i+1)+yp(i)+yp(i+1))*t.^3 + (-3*y(i)+3*y(i+1)-2*yp(i)-yp(i+1))*t.^2 + (yp(i))*t + y(i); + v = [X;Y]; + end + + g = @g_fun; +end + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/max_h.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,33 @@ +function [d_max, i1_max, j1_max, i2_max, j2_max] = max_h(X,Y) + ni = size(X,1); + nj = size(X,2); + d_max = 0; + + i1_max = 0; + j1_max = 0; + i2_max = 0; + j2_max = 0; + + D = {[0,-1],[1,0],[0,1],[-1,0]}; + + for i = 1:ni + for j = 1:nj + p1 = [X(i,j); Y(i,j)]; + for k = 1:length(D) + i2 = i+D{k}(1); + j2 = j+D{k}(2); + if i2 >= 1 && i2 <= ni && j2 >= 1 && j2 <= nj + p2 = [X(i2,j2); Y(i2,j2)]; + d = norm(p2-p1); + if d > d_max; + d_max = d; + i1_max = i; + j1_max = j; + i2_max = i2; + j2_max = j2; + end + end + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/min_h.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,34 @@ +function [d_min, i1_min, j1_min, i2_min, j2_min] = min_h(X,Y) + ni = size(X,1); + nj = size(X,2); + d_min = norm([X(1,1);Y(1,1)] - [X(ni,nj);Y(ni,nj)]); + + i1_min = 0; + j1_min = 0; + i2_min = 0; + j2_min = 0; + + D = {[-1,-1],[0,-1],[1,-1],[1,0],[1,1],[0,1],[-1,1],[-1,0]}; + % D = {[0,-1],[1,0],[0,1],[-1,0]}; + + for i = 1:ni + for j = 1:nj + p1 = [X(i,j); Y(i,j)]; + for k = 1:length(D) + i2 = i+D{k}(1); + j2 = j+D{k}(2); + if i2 >= 1 && i2 <= ni && j2 >= 1 && j2 <= nj + p2 = [X(i2,j2); Y(i2,j2)]; + d = norm(p2-p1); + if d < d_min; + d_min = d; + i1_min = i; + j1_min = j; + i2_min = i2; + j2_min = j2; + end + end + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/plot_shape.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,12 @@ +% Plot a shape using n points. Returns cell array of plot handles. +% hs = plot_shape(s,n) +function hs = plot_shape(s,n) + default_arg('n',100); + + hs = {}; + hold on + for i = 1:length(s) + hs{end+1} = grid.plot_curve(s{i},n); + end + hold off +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/shape.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,4 @@ +% Creates a shape from a number of curves. A shape is a cell array of curves. +function s = shape(varargin); + s = varargin; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/shape_discretise.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,8 @@ +% Discretises a shape such that points on the curves are no further than h appart. +function p = shape_discretise(s,h) + p = []; + for i = 1:length(s) + [~,pt] = grid.curve_discretise(s{i},h); + p = [p, pt]; + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/shape_linesegments.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,9 @@ +% Converts a shape into a cell array of linesegments shorter than h. +function l = shape_linesegments(s,h) + l = {}; + + for i = 1:length(s) + t = grid.curve_discretise(s{i},h); + l = [l, grid.curve_linesegments(s{i},t)]; + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_interp.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,132 @@ +classdef triang_interp + properties + g1, g2 ,g3 % Curves encirling the tirangle in the positive direction. + A,B,C % The corners of the triangle + Sa, Sb, Sc % Mappings from square with different sides collapsed + end + + methods + function o = triang_interp(g1,g2,g3) + o.g1 = g1; + o.g2 = g2; + o.g3 = g3; + o.A = g1(0); + o.B = g2(0); + o.C = g3(0); + o.Sa = grid.triang_interp.square_to_triangle_interp(g2,g3,g1); + o.Sb = grid.triang_interp.square_to_triangle_interp(g3,g1,g2); + o.Sc = grid.triang_interp.square_to_triangle_interp(g1,g2,g3); + end + + + function show(o,N) + % Show the mapped meridians of the triangle. + % Might be used for the barycentric coordinates. + ma = @(t)o.Sa(1/2,1-t); + mb = @(t)o.Sb(1/2,1-t); + mc = @(t)o.Sc(1/2,1-t); + + na = @(t)o.Sa(t,1/2); + + ka = @(t)(o.g1(1-t)+o.g2(t))/2; + + h = grid.plot_curve(ma); + h.Color = Color.blue; + h = grid.plot_curve(mb); + h.Color = Color.blue; + h = grid.plot_curve(mc); + h.Color = Color.blue; + + h = grid.plot_curve(na); + h.Color = Color.red; + + h = grid.plot_curve(ka); + h.Color = Color.red; + + [a(1),a(2)] = ma(1/3); + [b(1),b(2)] = mb(1/3); + [c(1),c(2)] = mc(1/3); + + d = ka(1-1/3); + + + grid.label_pt(a,b,c,d); + + + % t = linspace(0,1,N); + % for i = 1:N + % sa = @(s)o.Sa(s,t(i)); + % sb = @(s)o.Sb(s,t(i)); + % sc = @(s)o.Sc(s,t(i)); + + % h = grid.plot_curve(sa); + % h.Color = Color.blue; + % h = grid.plot_curve(sb); + % h.Color = Color.blue; + % h = grid.plot_curve(sc); + % h.Color = Color.blue; + % end + + h = grid.plot_curve(o.g1); + h.LineWidth = 2; + h.Color = Color.red; + + h = grid.plot_curve(o.g2); + h.LineWidth = 2; + h.Color = Color.red; + + h = grid.plot_curve(o.g3); + h.LineWidth = 2; + h.Color = Color.red; + + end + + + end + + methods(Static) + % Makes a mapping from the unit square to a triangle by collapsing + % one of the sides of the squares to a corner on the triangle + % The collapsed side is mapped to the corner oposite to g1. + % This is done such that for S(s,t), S(s,1) = g1(s) + function S = square_to_triangle_interp(g1,g2,g3) + corner = grid.line_segment(g3(0),g3(0)); + S = grid.transfinite_interp(corner,g3,f(g1),f(g2)) + + % Function to flip a curve + function h = f(g) + h = @(t)g(1-t); + end + end + end + +end + +% % Return a mapping from u.v to x,y of the domain encircled by g1 g2 g3 in the the positive direction. created be using transfinite interpolation. +% function S = triang_interp(g1,g2,g3) +% A = g1(0) +% B = g2(0) +% C = g3(0) + +% function [x,y] = S_fun(u,v) +% w = sqrt((u-1)^2+v^2)/sqrt(2); % Parameter for g3 +% v = v*(1-u-v)*g1(u) + u*(1-u-v)*g2(v) + u*v*g3(w) ... +% +(1-u)*(1-v)*A+u*(1-v)*B + (1-u)*v*C; +% x = v(1); +% y = v(2); +% end +% S = @S_fun; +% end + + + +% function subsref(obj,S) +% if ~all(isnumeric(S.subs{:})) +% error('Only supports calling object with number') +% end +% if numel(S.subs{:}) > 1 +% disp('You''ve called the object with more than one argument'); +% else +% disp(['You called the object with argument = ',num2str(S.subs{:})]); +% end +% end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_interp_pts.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,12 @@ +% Creates a transfinite interpolation from connecting the four points wiht straight lines. +function [S, g1, g2, g3] = triang_interp_pts(p1,p2,p3) + if size(p1) ~= [2 1] + error('p1 is strange!'); + end + + g1 = @(t)(p1 + t*(p2-p1)); + g2 = @(t)(p2 + t*(p3-p2)); + g3 = @(t)(p3 + t*(p1-p3)); + + S = grid.triang_interp(g1,g2,g3); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_map.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,29 @@ +% Creates a grid [X,Y] from the mapping function S at points in vectors u,v +function [X, Y] = traing_map(S,u,v) + error('not done') + if nargin == 2 + v = u; + end + + if isscalar(u) + u = linspace(0,1,u); + end + + if isscalar(v) + v = linspace(0,1,v); + end + + nu = length(u); + nv = length(v); + + X = zeros(nu,nv); + Y = zeros(nu,nv); + + for i = 1:nu + for j = 1:nv + [x,y] = S(u(i),v(j)); + X(i,j) = x; + Y(i,j) = y; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_plot_interp.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,114 @@ +% Plots a transfinite interpolation in x,y space using nu and nv curves along u and v axes. + + + + + + +% Plots a interp of a triangle where one the interpolation is from a square +% with one side collapsed to +function h = triang_plot_interp_kindaworking(S,n) + u = linspace(0,1,n); + v = linspace(0,1,n); + + m = 100; + m = 20; + + Xl_curves = cell(n,1); + Xr_curves = cell(n,1); + Y_curves = cell(n,1); + + + function u = wierdness(v,d,N) + if N == 0 + u = 0; + else + u = N*d./(1-v); + end + end + + + %Y curves + t = linspace(0,1,m); + for i = 1:n + x = []; y = []; + for j = 1:length(t) + [x(j),y(j)] = S(t(j),v(i)); + end + Y_curves{i} = [x', y']; + end + + + % Right and left X curves + t = linspace(0,1,m); + d = u(2); + for i = 1:n + xl = []; yl = []; + xr = []; yr = []; + N = i-1; + t = linspace(0,1-N*d,m); + for j = 1:length(t) + w = wierdness(t(j),d,N); + [xr(j),yr(j)] = S(w,t(j)); + [xl(j),yl(j)] = S(1-w,t(j)); + end + Xl_curves{i} = [xl', yl']; + Xr_curves{i} = [xr', yr']; + end + + for i = 1:n-1 + line(Xl_curves{i}(:,1),Xl_curves{i}(:,2)) + line(Xr_curves{i}(:,1),Xr_curves{i}(:,2)) + line(Y_curves{i}(:,1),Y_curves{i}(:,2)) + end +end + + + + +function h = triang_plot_interp_nonworking(S,n) + + u = linspace(0,1,n); + v = linspace(0,1,n); + + m = 100; + + X_curves = cell(n-1,1); + Y_curves = cell(n-1,1); + K_curves = cell(n-1,1); + + + t = linspace(0,1,m); + for i = 1:n-1 + x = []; y = []; + for j = find(t+u(i) <= 1) + [x(j),y(j)] = S(u(i),t(j)); + end + X_curves{i} = [x', y']; + end + + for i = 1:n-1 + x = []; y = []; + for j = find(t+v(i) <= 1) + [x(j),y(j)] = S(t(j),v(i)); + end + Y_curves{i} = [x', y']; + end + + for i = 2:n + x = []; y = []; + for j = find(t<u(i)) + [x(j),y(j)] = S(t(j), u(i)-t(j)); + end + K_curves{i-1} = [x', y']; + end + + for i = 1:n-1 + line(X_curves{i}(:,1),X_curves{i}(:,2)) + line(Y_curves{i}(:,1),Y_curves{i}(:,2)) + line(K_curves{i}(:,1),K_curves{i}(:,2)) + end + + h = -1; + % h = plot(X_curves{:},Y_curves{:},K_curves{:}); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_show.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,23 @@ +% Show a grid as a matlab figure. +function triang_show(S,n) + + ih = ishold(); + + hold on + h_grid = grid.triang_plot_interp(S,n); + h_bord = grid.triang_plot_interp(S,2); + + set(h_grid,'Color',[0 0.4470 0.7410]); + set(h_bord,'Color',[0.8500 0.3250 0.0980]); + set(h_bord,'LineWidth',2); + + % axis auto + % axis equal + % axis square + + if ih + hold on + else + hold off + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/place_label.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,12 @@ +% 'left' | 'center' | 'right' +% 'baseline' | 'top' | 'cap' | 'middle' | 'bottom' +function place_label(pt,str,horzAl,vertAl) + default_arg('horzAl','center'); + default_arg('vertAl', 'middle'); + + x = pt(1); + y = pt(2); + h = text(x,y,str); + h.HorizontalAlignment = horzAl; + h.VerticalAlignment = vertAl; +end \ No newline at end of file
--- a/+scheme/Scheme.m Tue Feb 07 14:38:51 2017 +0000 +++ b/+scheme/Scheme.m Tue Feb 07 15:47:51 2017 +0100 @@ -1,39 +1,42 @@ -% Start with all matrix returns. When that works see how we should generalize to non-matrix stuff/nonlinear +% Start with all matrix returns. When that works see how we should generalize +% to non-matrix stuff/nonlinear classdef Scheme < handle properties (Abstract) - m % Number of points in each direction, possibly a vector order % Order accuracy for the approximation - % vectors u,v,w depending on dim that gives were gridpoints are in each dimension - % vectors x,y,z containing the x,y,z values corresponding to each grid point - % matrices X,Y,Z with point coordinates as multi dimensional vectors + grid D % non-stabalized scheme operator H % Discrete norm - - % Should also containg: - % the grid points used - % the grid spacing end methods (Abstract) - % Closure functions return the opertors applied to the own doamin to close the boundary - % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. - % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. - % data is a function returning the data that should be applied at the boundary. - % neighbour_scheme is an instance of Scheme that should be interfaced to. - % neighbour_boundary is a string specifying which boundary to interface to. - m = boundary_condition(obj,boundary,type,data) - m = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - N = size(obj) % Returns the number of degrees of freedom. + % Closure functions return the opertors applied to the own doamin to + % close the boundary Penalty functions return the opertors to force + % the solution. In the case of an interface it returns the operator + % applied to the other doamin. In some cases the penalty return value + % can be ommited and the closure function take care of both parts. + % boundary is a string specifying the boundary e.g. + % 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of + % boundary condition if there are several. + % data is a function returning the data that + % should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be + % interfaced to. + % neighbour_boundary is a string specifying which boundary to + % interface to. + [closure, penalty] = boundary_condition(obj,boundary,type,data) + [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % Returns the number of degrees of freedom. + N = size(obj) end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % Calculates the matrcis need for the inteface coupling between + % boundary bound_u of scheme schm_u and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
--- a/cell2sparse.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,54 +0,0 @@ -function A = cell2sparse(C) - - if isempty(C) - A = sparse([]); - return - end - - n = row_height(C); - m = col_width(C); - - N = sum(n); - M = sum(m); - - A = sparse(N,M); - - n_ind = [0 cumsum(n)]; - m_ind = [0 cumsum(m)]; - - for i = 1:size(C,1) - for j = 1:size(C,2) - if ~has_matrix(C{i,j}) - continue - end - A(n_ind(i)+1:n_ind(i+1),m_ind(j)+1:m_ind(j+1)) = C{i,j}; - end - end - -end - -function m = col_width(C) - for j = 1:size(C,2) - for i = 1:size(C,1) - if ~has_matrix(C{i,j}) - continue - end - m(j) = size(C{i,j},2); - end - end -end - -function n = row_height(C) - for i = 1:size(C,1) - for j = 1:size(C,2) - if ~has_matrix(C{i,j}) - continue - end - n(i) = size(C{i,j},1); - end - end -end - -function b = has_matrix(c) - b = ~(isempty(c) || (numel(c)==1 && c == 0)); -end \ No newline at end of file
--- a/cell2vector.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -% cell2vector accepts a column cell array of column vectors and returns a columnvector -% with the input concatenated. It also returns the number of elements in each vector. -% cv -- column cell array with column vectors -% v -- vector of the concatenated vectors -% n -- number of elements in each vector before concatenation. Can be used with vector2cell(). -function [v, n] = cell2vector(cv) - v = []; - n = zeros(length(cv),1); - - for i = 1:length(cv) - n(i) = length(cv{i}); - v = [v; cv{i}]; - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeRowMaj.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,12 @@ +% Reshapes a matrix as if it was stored in row major order. +function B = reshapeRowMaj(A, m) + D = length(m); + + if D == 1 + m = [m 1]; + D = 2; + end + + % Reshape and reverse order of indecies + B = permute(reshape(permute(A, ndims(A):-1:1), rot90(m,2)), D:-1:1); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeRowMajTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,113 @@ +function tests = reshapeRowMajTest() + tests = functiontests(localfunctions); +end + +function test1D(testCase) + in = { + {5,[1; 2; 3; 4; 5]}, + {5,[1 2 3 4 5]}, + }; + out = { + [1; 2; 3; 4; 5], + [1; 2; 3; 4; 5], + }; + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end + + +function testIdentity(testCase) + in = { + {[2,2], magic(2)}, + {[3,3], magic(3)}, + {[2,3], [1 2 3; 4 5 6]}, + }; + + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),in{i}{2}); + end +end + +function test2D(testCase) + in = { + {[2,2],[11; 12; 21; 22]}, + {[3,2],[1 2 3; 4 5 6]}, + {[6 1],[1 2 3; 4 5 6]}, + {[1 6],[1 2 3; 4 5 6]}, + }; + + out{1}(1,1) = 11; + out{1}(1,2) = 12; + out{1}(2,1) = 21; + out{1}(2,2) = 22; + + out{2} = [1 2; 3 4; 5 6]; + out{3} = [1; 2; 3; 4; 5; 6]; + out{4} = [1 2 3 4 5 6]; + + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end + +function test3D(testCase) + in = { + {[2, 2, 2], [111; 112; 121; 122; 211; 212; 221; 222]}, + {[8 1], cat(3,[1 2; 3 4],[5 6; 7 8])}, + {[1 8], cat(3,[1 2; 3 4],[5 6; 7 8])}, + {[2 4], cat(3,[1 2; 3 4],[5 6; 7 8])}, + {[4 2], cat(3,[1 2; 3 4],[5 6; 7 8])}, + }; + + out{1}(1,1,1) = 111; + out{1}(1,1,2) = 112; + out{1}(1,2,1) = 121; + out{1}(1,2,2) = 122; + out{1}(2,1,1) = 211; + out{1}(2,1,2) = 212; + out{1}(2,2,1) = 221; + out{1}(2,2,2) = 222; + + out{2} = [1; 5; 2; 6; 3; 7; 4; 8]; + out{3} = [1 5 2 6 3 7 4 8]; + out{4} = [1 5 2 6; 3 7 4 8]; + out{5} = [1 5; 2 6; 3 7; 4 8]; + + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end + +function testNonSquare(testCase) + in = { + {[2, 3, 4],[111; 112; 113; 114; 121; 122; 123; 124; 131; 132; 133; 134; 211; 212; 213; 214; 221; 222; 223; 224; 231; 232; 233; 234]}, + }; + out{1}(1,1,1) = 111; + out{1}(1,1,2) = 112; + out{1}(1,1,3) = 113; + out{1}(1,1,4) = 114; + out{1}(1,2,1) = 121; + out{1}(1,2,2) = 122; + out{1}(1,2,3) = 123; + out{1}(1,2,4) = 124; + out{1}(1,3,1) = 131; + out{1}(1,3,2) = 132; + out{1}(1,3,3) = 133; + out{1}(1,3,4) = 134; + out{1}(2,1,1) = 211; + out{1}(2,1,2) = 212; + out{1}(2,1,3) = 213; + out{1}(2,1,4) = 214; + out{1}(2,2,1) = 221; + out{1}(2,2,2) = 222; + out{1}(2,2,3) = 223; + out{1}(2,2,4) = 224; + out{1}(2,3,1) = 231; + out{1}(2,3,2) = 232; + out{1}(2,3,3) = 233; + out{1}(2,3,4) = 234; + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeToPlotMatrix.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,17 @@ +% Takes a grid function and reshapes it into a matrix of shape m for plotting. +function F = reshapeToPlotMatrix(gf, m) + D = length(m); + + switch D + case 1 + F = gf; + case 2 + F = reshape(gf, rot90(m,2)); + case 3 + % After the reshape the indecies will be M(z,y,x). Plot need them to be M(y,x,z) + p = [2 3 1]; % Permuation + F = permute(reshape(gf,rot90(m,2)), p); + otherwise + error('reshapeToPlotMatrix:NotImplemented','Grid function to matrix is not implemented for dimension = %d', length(m)); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeToPlotMatrixTest.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,48 @@ +function tests = reshapeToPlotMatrixTest() + tests = functiontests(localfunctions); +end + +function test1D(testCase) + inGf = [1 2 3 4 5]'; + inM = 5; + out = [1 2 3 4 5]'; + testCase.verifyEqual(reshapeToPlotMatrix(inGf, inM),out); +end + +function test2D(testCase) + x = 1:2; + y = 1:3; + + f = @(x,y) x + y*10; + + xx = [1; 1; 1; 2; 2; 2]; + yy = [1; 2; 3; 1; 2; 3]; + inGf = f(xx,yy); + + [X,Y] = meshgrid(x,y); + out = f(X,Y); + + inM = [2, 3]; + + testCase.verifyEqual(reshapeToPlotMatrix(inGf, inM),out); +end + +function test3D(testCase) + x = 1:2; + y = 1:3; + z = 1:4; + + f = @(x,y,z) x + y*10 + z*100; + + xx = [repmat(1, [12, 1]); repmat(2, [12, 1])]; + yy = repmat([1; 1; 1; 1; 2; 2; 2; 2; 3; 3; 3; 3], [2, 1]); + zz = repmat([1; 2; 3; 4], [6, 1]); + inGf = f(xx,yy,zz); + + [X,Y,Z] = meshgrid(x,y,z); + out = f(X,Y,Z); + + inM = [2, 3, 4]; + + testCase.verifyEqual(reshapeToPlotMatrix(inGf, inM),out); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/runtestsPackage.m Tue Feb 07 15:47:51 2017 +0100 @@ -0,0 +1,4 @@ +function res = runtestsPackage(pkgName) + ts = matlab.unittest.TestSuite.fromPackage(pkgName); + res = ts.run(); +end \ No newline at end of file
--- a/vector2cell.m Tue Feb 07 14:38:51 2017 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -% Splits column vector v into segments of length n and returns the result as a column cell array. -% v -- column vector to be split -% n -- number of elements in each part -% -% cv -- cell array of vectors with lenght n(i) -function cv = vector2cell(v,n) - cv = cell(length(n),1); - - ind = [0; cumsum(n)]; - for i = 1:length(n) - ind_i = (ind(i)+1):ind(i+1); - cv{i} = v(ind_i); - end -end \ No newline at end of file