From 21a93bca5f57a2c22a729ada4ddc82e20c9480fe Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Thu, 26 Mar 2026 21:17:57 -0700 Subject: [PATCH 01/10] test --- test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 test.txt diff --git a/test.txt b/test.txt new file mode 100644 index 000000000..e69de29bb From d44483c7bdc518599b8cc513d1e16c3559b9983f Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Thu, 26 Mar 2026 21:35:54 -0700 Subject: [PATCH 02/10] jacobian --- src/matlab_octave/jacobian2D.m | 107 +++++++++++++++++++-------- src/matlab_octave/jacobian2DLegacy.m | 37 +++++++++ test.txt | 0 3 files changed, 113 insertions(+), 31 deletions(-) create mode 100644 src/matlab_octave/jacobian2DLegacy.m delete mode 100644 test.txt diff --git a/src/matlab_octave/jacobian2D.m b/src/matlab_octave/jacobian2D.m index b0d0a2be8..81860b12d 100644 --- a/src/matlab_octave/jacobian2D.m +++ b/src/matlab_octave/jacobian2D.m @@ -1,44 +1,89 @@ -function [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y) +function [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, m, dx, n, dy, dc, nc) % PURPOSE -% +% 2D Jacobian Metrics for Curvilinear Operators +% % DESCRIPTION % Returns: -% J : Determinant of the Jacobian (XeYn - XnYe) -% Xe : dx/de metric -% Xn : dx/dn metric -% Ye : dy/de metric -% Yn : dy/dn metric +% J : Determinant of the Jacobian (XeYn - XnYe) on the +% centers if optional arguments are specified, else nodes +% Xe : dx/de metric on the centers if optional arguments are +% specified, else nodes +% Xn : dx/dn metric on the centers if optional arguments are +% specified, else nodes +% Ye : dy/de metric on the centers if optional arguments are +% specified, else nodes +% Yn : dy/dn metric on the centers if optional arguments are +% specified, else nodes % % Parameters: % k : Order of accuracy -% X : x-coordinates (physical) of meshgrid -% Y : y-coordinates (physical) of meshgrid -% -% SYNTAX -% [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y) -% +% X : x-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% Y : y-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% (optional) m : Number of cells in xi direction +% (optional) dx : Step size in xi direction +% (optional) n : Number of cells in eta direction +% (optional) dy : Step size in eta direction +% (optional) dc : a0 (4x1 vector for left, right, bottom, top +% boundaries, resp.) +% (optional) nc : b0 (4x1 vector for left, right, bottom, top +% boundaries, resp.) % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - [n, m] = size(X); - - X = reshape(X', [], 1); - Y = reshape(Y', [], 1); - - N = nodal2D(k, m, 1, n, 1); - - X = N*X; - Y = N*Y; - - mn = m*n; - - Xe = X(1:mn); - Xn = X(mn+1:end); - Ye = Y(1:mn); - Yn = Y(mn+1:end); - - J = Xe.*Yn-Xn.*Ye; + if nargin ~= 3 && nargin ~= 9 + error("jacobian2D:InvalidNumArgs", "jacobian2D expects 3 or 9 arguments") + elseif nargin == 3 + [J,Xe,Xn,Ye,Yn] = jacobian2DLegacy(k,X,Y); + return; + end + + assert(all(size(dc) == [4 1]), "dc is a 4x1 vector") + assert(all(size(nc) == [4 1]), "nc is a 4x1 vector") + + % Periodic Handling + if isempty(find(dc(1:2).^2 + nc(1:2).^2,1)) + Ge = gradPeriodic(k,m,dx); + IFCx = interpolFacesToCentersG1DPeriodic(k,m); + Im = speye(m); + else + Ge = gradNonPeriodic(k,m,dx); + IFCx = interpolFacesToCentersG1D(k,m); + Im = speye(m+2); + end + if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) + Gn = gradPeriodic(k,n,dy); + IFCy = interpolFacesToCentersG1DPeriodic(k,n); + In = speye(n); + else + Gn = gradNonPeriodic(k,n,dy); + IFCy = interpolFacesToCentersG1D(k,n); + In = speye(n+2); + end + + numC = numel(X); + + % Get metrics on centers + % We don't augment the identity matrices + % so that we don't lose information + % around the boundaries + Ge = kron(In,IFCx) * kron(In,Ge); + Gn = kron(IFCy,Im) * kron(Gn,Im); + G = [Ge; Gn]; + + X = reshape(X',[],1); + Y = reshape(Y',[],1); + metrics = G * [X Y]; + + Xe = metrics(1:numC,1); + Xn = metrics(1+numC:end,1); + Ye = metrics(1:numC,2); + Yn = metrics(1+numC:end,2); + + J = Xe .* Yn - Xn .* Ye; + end \ No newline at end of file diff --git a/src/matlab_octave/jacobian2DLegacy.m b/src/matlab_octave/jacobian2DLegacy.m new file mode 100644 index 000000000..ddafd3959 --- /dev/null +++ b/src/matlab_octave/jacobian2DLegacy.m @@ -0,0 +1,37 @@ +function [J, Xe, Xn, Ye, Yn] = jacobian2DLegacy(k, X, Y) +% Returns: +% J : Determinant of the Jacobian (XeYn - XnYe) on the nodes +% Xe : dx/de metric on the nodes +% Xn : dx/dn metric on the nodes +% Ye : dy/de metric on the nodes +% Yn : dy/dn metric on the nodes +% +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid nodes +% Y : y-coordinates (physical) of meshgrid nodes +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- + + [n, m] = size(X); + + X = reshape(X', [], 1); + Y = reshape(Y', [], 1); + + N = nodal2D(k, m, 1, n, 1); + + X = N*X; + Y = N*Y; + + mn = m*n; + + Xe = X(1:mn); + Xn = X(mn+1:end); + Ye = Y(1:mn); + Yn = Y(mn+1:end); + + J = Xe.*Yn-Xn.*Ye; +end \ No newline at end of file diff --git a/test.txt b/test.txt deleted file mode 100644 index e69de29bb..000000000 From 300eab2cbd82afd0c512004bac677109078bb107 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Thu, 26 Mar 2026 21:44:03 -0700 Subject: [PATCH 03/10] gradient --- src/matlab_octave/grad2DCurv.m | 121 ++++++++++++++++----------- src/matlab_octave/grad2DCurvLegacy.m | 58 +++++++++++++ src/matlab_octave/jacobian2D.m | 3 +- 3 files changed, 130 insertions(+), 52 deletions(-) create mode 100644 src/matlab_octave/grad2DCurvLegacy.m diff --git a/src/matlab_octave/grad2DCurv.m b/src/matlab_octave/grad2DCurv.m index d765d5376..7aabd3263 100644 --- a/src/matlab_octave/grad2DCurv.m +++ b/src/matlab_octave/grad2DCurv.m @@ -1,65 +1,84 @@ -function G = grad2DCurv(k, X, Y) +function G = grad2DCurv(k, X, Y, m, dx, n, dy, dc, nc) % PURPOSE -% Returns a 2D curvilinear mimetic gradient +% Returns a 2D curvilinear mimetic gradient operator % % DESCRIPTION +% Returns: +% G : 2D curvilinear mimetic gradient operator. It outpus +% to the extended faces (normal faces plus boundaries) +% +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% Y : y-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% (optional) m : Number of cells in xi direction +% (optional) dx : Step size in xi direction +% (optional) n : Number of cells in eta direction +% (optional) dy : Step size in eta direction +% (optional) dc : a0 (4x1 vector for left, right, bottom, top +% boundaries, resp.) +% (optional) nc : b0 (4x1 vector for left, right, bottom, top +% boundaries, resp.) % % SYNTAX % G = grad2DCurv(k, X, Y) +% G = grad2DCurv(k, X, Y, m, dx, n, dy, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - % Get the determinant of the jacobian and the metrics - [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y); - - % Dimensions of nodal grid - [n, m] = size(X); - - % Make them surfaces so they can be interpolated - J = reshape(J, m, n)'; - Xe = reshape(Xe, m, n)'; - Xn = reshape(Xn, m, n)'; - Ye = reshape(Ye, m, n)'; - Yn = reshape(Yn, m, n)'; - - % Logical grids - [Xl, Yl] = meshgrid(1:m, 1:n); - - % Interpolate the metrics on the logical grid for positions u and v - Ju = interp2(Xl, Yl, J, (Xl(1:end-1, :)+Xl(2:end, :))/2,... - (Yl(1:end-1, :)+Yl(2:end, :))/2); - Jv = interp2(Xl, Yl, J, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,... - (Yl(:, 1:end-1)+Yl(:, 2:end))/2); - Xev = interp2(Xl, Yl, Xe, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,... - (Yl(:, 1:end-1)+Yl(:, 2:end))/2); - Xnv = interp2(Xl, Yl, Xn, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,... - (Yl(:, 1:end-1)+Yl(:, 2:end))/2); - Yeu = interp2(Xl, Yl, Ye, (Xl(1:end-1, :)+Xl(2:end, :))/2,... - (Yl(1:end-1, :)+Yl(2:end, :))/2); - Ynu = interp2(Xl, Yl, Yn, (Xl(1:end-1, :)+Xl(2:end, :))/2,... - (Yl(1:end-1, :)+Yl(2:end, :))/2); - - % Convert metrics to diagonal matrices so they can be multiplied by the - % logical operators - Ju = spdiags(1./reshape(Ju', [], 1), 0, numel(Ju), numel(Ju)); - Jv = spdiags(1./reshape(Jv', [], 1), 0, numel(Jv), numel(Jv)); - Xev = spdiags(reshape(Xev', [], 1), 0, numel(Xev), numel(Xev)); - Xnv = spdiags(reshape(Xnv', [], 1), 0, numel(Xnv), numel(Xnv)); - Yeu = spdiags(reshape(Yeu', [], 1), 0, numel(Yeu), numel(Yeu)); - Ynu = spdiags(reshape(Ynu', [], 1), 0, numel(Ynu), numel(Ynu)); - - % Construct 2D uniform mimetic gradient operator (d/de, d/dn) - G = grad2D(k, m-1, 1, n-1, 1); - Ge = G(1:m*(n-1), :); - Gn = G(m*(n-1)+1:end, :); - - % Apply transformation - Gx = Ju*(Ynu*Ge-Yeu*GI2(Gn, m-1, n-1, 'Gn')); - Gy = Jv*(-Xnv*GI2(Ge, m-1, n-1, 'Ge')+Xev*Gn); - - % Final 2D curvilinear mimetic gradient operator (d/dx, d/dy) + + if nargin ~=3 && nargin ~= 9 + error("grad2DCurv:InvalidNumArgs", ... + "grad2DCurv expects 3 or 9 arguments") + elseif nargin == 3 + G = grad2DCurvLegacy(k,X,Y); + return; + end + + % Periodic Handling + if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) + Ge = gradPeriodic(k,m,dx); + ICFx = interpolCentersToFacesD1DPeriodic(k,m); + IFCx = interpolFacesToCentersG1DPeriodic(k,m); + Im = speye(m); + else + Ge = gradNonPeriodic(k,m,dx); + ICFx = interpolCentersToFacesD1D(k,m); + IFCx = interpolFacesToCentersG1D(k,m); + Im = speye(m+2); + end + if isempty(find(dc(3:4).^2 + nc(3:4).^2, 1)) + Gn = gradPeriodic(k,n,dy); + ICFy = interpolCentersToFacesD1DPeriodic(k,n); + IFCy = interpolFacesToCentersG1DPeriodic(k,n); + In = speye(n); + else + Gn = gradNonPeriodic(k,n,dy); + ICFy = interpolCentersToFacesD1D(k,n); + IFCy = interpolFacesToCentersG1D(k,n); + In = speye(n+2); + end + + % Make Ge and Gn act on and output to the centers + % This allows them to be added together + Ge = kron(In,IFCx) * kron(In,Ge); + Gn = kron(IFCy,Im) * kron(Gn,Im); + + % Apply Metrics + [J,Xe,Xn,Ye,Yn] = jacobian2D(k,X,Y,m,dx,n,dy,dc,nc); + + Gx = (Yn ./ J) .* Ge - (Ye ./ J) .* Gn; + Gy = (Xe ./ J) .* Gn - (Xn ./ J) .* Ge; + + % Now have them output to the extended faces + Gx = kron(In,ICFx) * Gx; + Gy = kron(ICFy,Im) * Gy; + G = [Gx; Gy]; + end diff --git a/src/matlab_octave/grad2DCurvLegacy.m b/src/matlab_octave/grad2DCurvLegacy.m new file mode 100644 index 000000000..9b0b1f359 --- /dev/null +++ b/src/matlab_octave/grad2DCurvLegacy.m @@ -0,0 +1,58 @@ +function G = grad2DCurvLegacy(k, X, Y) +% Returns a 2D curvilinear mimetic gradient +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- + % Get the determinant of the jacobian and the metrics + [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y); + + % Dimensions of nodal grid + [n, m] = size(X); + + % Make them surfaces so they can be interpolated + J = reshape(J, m, n)'; + Xe = reshape(Xe, m, n)'; + Xn = reshape(Xn, m, n)'; + Ye = reshape(Ye, m, n)'; + Yn = reshape(Yn, m, n)'; + + % Logical grids + [Xl, Yl] = meshgrid(1:m, 1:n); + + % Interpolate the metrics on the logical grid for positions u and v + Ju = interp2(Xl, Yl, J, (Xl(1:end-1, :)+Xl(2:end, :))/2,... + (Yl(1:end-1, :)+Yl(2:end, :))/2); + Jv = interp2(Xl, Yl, J, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,... + (Yl(:, 1:end-1)+Yl(:, 2:end))/2); + Xev = interp2(Xl, Yl, Xe, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,... + (Yl(:, 1:end-1)+Yl(:, 2:end))/2); + Xnv = interp2(Xl, Yl, Xn, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,... + (Yl(:, 1:end-1)+Yl(:, 2:end))/2); + Yeu = interp2(Xl, Yl, Ye, (Xl(1:end-1, :)+Xl(2:end, :))/2,... + (Yl(1:end-1, :)+Yl(2:end, :))/2); + Ynu = interp2(Xl, Yl, Yn, (Xl(1:end-1, :)+Xl(2:end, :))/2,... + (Yl(1:end-1, :)+Yl(2:end, :))/2); + + % Convert metrics to diagonal matrices so they can be multiplied by the + % logical operators + Ju = spdiags(1./reshape(Ju', [], 1), 0, numel(Ju), numel(Ju)); + Jv = spdiags(1./reshape(Jv', [], 1), 0, numel(Jv), numel(Jv)); + Xev = spdiags(reshape(Xev', [], 1), 0, numel(Xev), numel(Xev)); + Xnv = spdiags(reshape(Xnv', [], 1), 0, numel(Xnv), numel(Xnv)); + Yeu = spdiags(reshape(Yeu', [], 1), 0, numel(Yeu), numel(Yeu)); + Ynu = spdiags(reshape(Ynu', [], 1), 0, numel(Ynu), numel(Ynu)); + + % Construct 2D uniform mimetic gradient operator (d/de, d/dn) + G = grad2D(k, m-1, 1, n-1, 1); + Ge = G(1:m*(n-1), :); + Gn = G(m*(n-1)+1:end, :); + + % Apply transformation + Gx = Ju*(Ynu*Ge-Yeu*GI2(Gn, m-1, n-1, 'Gn')); + Gy = Jv*(-Xnv*GI2(Ge, m-1, n-1, 'Ge')+Xev*Gn); + + % Final 2D curvilinear mimetic gradient operator (d/dx, d/dy) + G = [Gx; Gy]; +end diff --git a/src/matlab_octave/jacobian2D.m b/src/matlab_octave/jacobian2D.m index 81860b12d..5e064335a 100644 --- a/src/matlab_octave/jacobian2D.m +++ b/src/matlab_octave/jacobian2D.m @@ -36,7 +36,8 @@ % ---------------------------------------------------------------------------- if nargin ~= 3 && nargin ~= 9 - error("jacobian2D:InvalidNumArgs", "jacobian2D expects 3 or 9 arguments") + error("jacobian2D:InvalidNumArgs", ... + "jacobian2D expects 3 or 9 arguments") elseif nargin == 3 [J,Xe,Xn,Ye,Yn] = jacobian2DLegacy(k,X,Y); return; From 29b5863c0b578a6e6d79e089a2c115325ef047b6 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Thu, 26 Mar 2026 21:53:49 -0700 Subject: [PATCH 04/10] divergence --- src/matlab_octave/div2DCurv.m | 122 ++++++++++++++++++---------- src/matlab_octave/div2DCurvLegacy.m | 54 ++++++++++++ src/matlab_octave/grad2DCurv.m | 8 +- 3 files changed, 137 insertions(+), 47 deletions(-) create mode 100644 src/matlab_octave/div2DCurvLegacy.m diff --git a/src/matlab_octave/div2DCurv.m b/src/matlab_octave/div2DCurv.m index 039161245..805b0b0bb 100644 --- a/src/matlab_octave/div2DCurv.m +++ b/src/matlab_octave/div2DCurv.m @@ -1,11 +1,31 @@ -function D = div2DCurv(k, X, Y) +function D = div2DCurv(k, X, Y, m, dx, n, dy, dc, nc) % PURPOSE -% Returns a 2D curvilinear mimetic divergence +% Returns a 2D curvilinear mimetic divergence operator % % DESCRIPTION +% Returns: +% D : 2D curvilinear mimetic divergence operator. If optional +% arguments are specified, it acts on the extended faces +% (normal faces plus boundaries) +% +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% Y : y-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% (optional) m : Number of cells in xi direction +% (optional) dx : Step size in xi direction +% (optional) n : Number of cells in eta direction +% (optional) dy : Step size in eta direction +% (optional) dc : a0 (4x1 vector for left, right, bottom, top +% boundaries, resp.) +% (optional) nc : b0 (4x1 vector for left, right, bottom, top +% boundaries, resp.) % % SYNTAX % D = div2DCurv(k, X, Y) +% D = div2DCurv(k, X, Y, m, dx, n, dy, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -13,48 +33,60 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - % Get the determinant of the jacobian and the metrics - [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y); - - % Dimensions of nodal grid - [n, m] = size(X); - - % Make them surfaces so they can be interpolated - J = reshape(J, m, n)'; - Xe = reshape(Xe, m, n)'; - Xn = reshape(Xn, m, n)'; - Ye = reshape(Ye, m, n)'; - Yn = reshape(Yn, m, n)'; - - % Logical grids - [Xl, Yl] = meshgrid(1:m, 1:n); - % Staggered logical grid - [Xs, Ys] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n]); - - % Interpolate the metrics on the logical grid for positions (Xs, Ys) - J = interp2(Xl, Yl, J, Xs, Ys); - Xe = interp2(Xl, Yl, Xe, Xs, Ys); - Xn = interp2(Xl, Yl, Xn, Xs, Ys); - Ye = interp2(Xl, Yl, Ye, Xs, Ys); - Yn = interp2(Xl, Yl, Yn, Xs, Ys); - - % Convert metrics to diagonal matrices so they can be multiplied by the - % logical operators - J = spdiags(1./reshape(J', [], 1), 0, numel(J), numel(J)); - Xe = spdiags(reshape(Xe', [], 1), 0, numel(Xe), numel(Xe)); - Xn = spdiags(reshape(Xn', [], 1), 0, numel(Xn), numel(Xn)); - Ye = spdiags(reshape(Ye', [], 1), 0, numel(Ye), numel(Ye)); - Yn = spdiags(reshape(Yn', [], 1), 0, numel(Yn), numel(Yn)); - - % Construct 2D uniform mimetic divergence operator (d/de, d/dn) - D = div2D(k, m-1, 1, n-1, 1); - De = D(:, 1:m*(n-1)); - Dn = D(:, m*(n-1)+1:end); - - % Apply transformation - Dx = J*(Yn*De-Ye*DI2(m-1, n-1, 'Dn')); - Dy = J*(-Xn*DI2(m-1, n-1, 'De')+Xe*Dn); - - % Final 2D curvilinear mimetic divergence operator (d/dx, d/dy) + if nargin ~= 3 && nargin ~= 9 + error("div2DCurv:InvalidNumArgs", ... + "div2DCurv expects 3 or 9 arguments") + elseif nargin == 3 + D = div2DCurvLegacy(k,X,Y); + return; + end + + assert(all(size(dc) == [4 1]), "dc is a 4x1 vector") + assert(all(size(nc) == [4 1]), "nc is a 4x1 vector") + + % Periodic Handling + if isempty(find(dc(1:2).^2 + nc(1:2).^2,1)) + De = divPeriodic(k,m,dx); + IFCx = interpolFacesToCentersG1DPeriodic(k,m); + ICFx = interpolCentersToFacesD1DPeriodic(k,m); + Im = speye(m); + else + De = divNonPeriodic(k,m,dx); + IFCx = interpolFacesToCentersG1D(k,m); + ICFx = interpolCentersToFacesD1D(k,m); + Im = speye(m+2); + end + if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) + Dn = divPeriodic(k,n,dy); + IFCy = interpolFacesToCentersG1DPeriodic(k,n); + ICFy = interpolCentersToFacesD1DPeriodic(k,n); + In = speye(n); + else + Dn = divNonPeriodic(k,n,dy); + IFCy = interpolFacesToCentersG1D(k,n); + ICFy = interpolCentersToFacesD1D(k,n); + In = speye(n+2); + end + + % Make De and Dn act on and output to the centers + % This allows them to be added together + De = kron(In,De) * kron(In,ICFx); + Dn = kron(Dn,Im) * kron(ICFy,Im); + + % Apply metrics + [J,Xe,Xn,Ye,Yn] = jacobian2D(k,X,Y,m,dx,n,dy,dc,nc); + + Dx = (Yn ./ J) .* De - (Ye ./ J) .* Dn; + Dy = (Xe ./ J) .* Dn - (Xn ./ J) .* De; + + % Now have them act on the extended faces + Dx = Dx * kron(In,IFCx); + Dy = Dy * kron(IFCy,Im); + D = [Dx Dy]; + + % Ensure no output on boundary -- Probably a cleaner way to do this + bdry = find(sum(spones(div2D(2,m,1,n,1,dc,nc)), 2) == 0); + D(bdry,:) = sparse(size(bdry,1),size(D,2)); + end diff --git a/src/matlab_octave/div2DCurvLegacy.m b/src/matlab_octave/div2DCurvLegacy.m new file mode 100644 index 000000000..f0026a84a --- /dev/null +++ b/src/matlab_octave/div2DCurvLegacy.m @@ -0,0 +1,54 @@ +function D = div2DCurvLegacy(k, X, Y) +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- + +% Returns a 2D curvilinear mimetic divergence + + % Get the determinant of the jacobian and the metrics + [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y); + + % Dimensions of nodal grid + [n, m] = size(X); + + % Make them surfaces so they can be interpolated + J = reshape(J, m, n)'; + Xe = reshape(Xe, m, n)'; + Xn = reshape(Xn, m, n)'; + Ye = reshape(Ye, m, n)'; + Yn = reshape(Yn, m, n)'; + + % Logical grids + [Xl, Yl] = meshgrid(1:m, 1:n); + % Staggered logical grid + [Xs, Ys] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n]); + + % Interpolate the metrics on the logical grid for positions (Xs, Ys) + J = interp2(Xl, Yl, J, Xs, Ys); + Xe = interp2(Xl, Yl, Xe, Xs, Ys); + Xn = interp2(Xl, Yl, Xn, Xs, Ys); + Ye = interp2(Xl, Yl, Ye, Xs, Ys); + Yn = interp2(Xl, Yl, Yn, Xs, Ys); + + % Convert metrics to diagonal matrices so they can be multiplied by the + % logical operators + J = spdiags(1./reshape(J', [], 1), 0, numel(J), numel(J)); + Xe = spdiags(reshape(Xe', [], 1), 0, numel(Xe), numel(Xe)); + Xn = spdiags(reshape(Xn', [], 1), 0, numel(Xn), numel(Xn)); + Ye = spdiags(reshape(Ye', [], 1), 0, numel(Ye), numel(Ye)); + Yn = spdiags(reshape(Yn', [], 1), 0, numel(Yn), numel(Yn)); + + % Construct 2D uniform mimetic divergence operator (d/de, d/dn) + D = div2D(k, m-1, 1, n-1, 1); + De = D(:, 1:m*(n-1)); + Dn = D(:, m*(n-1)+1:end); + + % Apply transformation + Dx = J*(Yn*De-Ye*DI2(m-1, n-1, 'Dn')); + Dy = J*(-Xn*DI2(m-1, n-1, 'De')+Xe*Dn); + + % Final 2D curvilinear mimetic divergence operator (d/dx, d/dy) + D = [Dx Dy]; +end diff --git a/src/matlab_octave/grad2DCurv.m b/src/matlab_octave/grad2DCurv.m index 7aabd3263..d329200b9 100644 --- a/src/matlab_octave/grad2DCurv.m +++ b/src/matlab_octave/grad2DCurv.m @@ -4,8 +4,9 @@ % % DESCRIPTION % Returns: -% G : 2D curvilinear mimetic gradient operator. It outpus -% to the extended faces (normal faces plus boundaries) +% G : 2D curvilinear mimetic gradient operator. If optional +% arguments are specified, it outputs to the extended +% faces (normal faces plus boundaries) % % Parameters: % k : Order of accuracy @@ -40,6 +41,9 @@ return; end + assert(all(size(dc) == [4 1]), "dc is a 4x1 vector") + assert(all(size(nc) == [4 1]), "nc is a 4x1 vector") + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) Ge = gradPeriodic(k,m,dx); From f1010faa4b4a5e92a607756b5208b6a41ac26ae2 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Thu, 26 Mar 2026 21:54:08 -0700 Subject: [PATCH 05/10] Comparison Example --- examples/matlab_octave/curvConvergenceTest.m | 222 +++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 examples/matlab_octave/curvConvergenceTest.m diff --git a/examples/matlab_octave/curvConvergenceTest.m b/examples/matlab_octave/curvConvergenceTest.m new file mode 100644 index 000000000..564bb0d2c --- /dev/null +++ b/examples/matlab_octave/curvConvergenceTest.m @@ -0,0 +1,222 @@ +%% Convegence Test for Cuvilinear Operators +% +% In 2D space, MOLE operators can only calculate the xi and eta derivatives +% On a curvilinear domain, the x and xi, and y and eta directions are +% not always the same. We can assume a transformation of x = x(xi, eta) +% and y = y(xi, eta) with Jacobian J = x_xi y_eta - x_eta y_xi. +% +% If we have a function u(x,y), via chain rule we get: +% u_xi = u_x x_xi + u_y y_xi +% u_eta = u_x x_eta + u_y y_eta +% +% Assuming a non-zero Jacobian, we can invert the mapping to obtain: +% u_x = (u_xi y_eta - u_eta y_xi) / J +% u_y = (u_eta x_xi - u_xi x_eta) / J +% +% Everything on the right is easily calculable using MOLE's operators +% + +%% Comparison between current implementation and new implementation +% +% Current Implementation: +% - Nodal Jacobian derivatives +% - MATLAB interpolation (2nd order) +% - Hard to understand implementation (DI2 and GI2 functions: It is +% obvious what they are doing, but not how they are doing it.) +% - Cannot handle periodic domains +% +% New Implementation +% - MOLE gradient used to calculate Jacobian derivatives +% - MOLE interpolators +% - Easy to understand implementation +% - Can handle periodic domains +% + +%% Problem +% +% ∆u = 0 +% +% 1 < R < 2 +% 0 < θ < pi +% +% x = R cos(θ) +% y = R sin(θ) +% +% Dirichlet Boundary Conditions +% +% Exact Solution: (many options) +% u(x,y) = sin(x) sinh(y) +% u(x,y) = x^2 - y^2 +% u(x,y) = ln(x^2 + y^2) +% u(x,y) = x^3 - 3 x y^2 +% +close all; clear; clc; + + +%% Fixed mesh, varying k Comparison +% Parameters +k = [2 4 6 8]; +m = 100; +n = 25; +numCenters = (m+2) * (n+2); +dx = pi / m; +dy = 1 / n; + +% ue = @(X,Y) sin(X) .* sinh(Y); +% ue = @(X,Y) X.^2 - Y.^2; +% ue = @(X,Y) log(X.^2 + Y.^2); +ue = @(X,Y) X.^3 - 3 * X .* Y.^2; + +errorCur = zeros(size(k)); +errorNew = zeros(size(k)); + +dc = [1;1;1;1]; nc = [0;0;0;0]; + +for i = 1:numel(k) + + fprintf("Starting k = %d\n", k(i)) + % Nodes to make current operators + rsN = 1:dy:2; tsN = pi:-dx:0; + [TSN,RSN] = meshgrid(tsN,rsN); + xn = RSN .* cos(TSN); + yn = RSN .* sin(TSN); + + % Centers to make new operators + rsC = [1 1+dy/2:dy:2-dy/2 2]; + tsC = [pi pi-dx/2:-dx:dx/2 0]; + [TSC,RSC] = meshgrid(tsC,rsC); + xc = RSC .* cos(TSC); xc = reshape(xc',[],1); + yc = RSC .* sin(TSC); yc = reshape(yc',[],1); + + % Build operators + curG = grad2DCurv(k(i),xn,yn); + curD = div2DCurv(k(i),xn,yn); + curL = curD * curG; + + newG = grad2DCurv(k(i),xc,yc,m,dx,n,dy,dc,nc); + newD = div2DCurv(k(i),xc,yc,m,dx,n,dy,dc,nc); + newL = newD * newG; + + % Boundary Conditions + X = reshape(xc,m+2,n+2)'; % Reshape for plotting and easier BC + Y = reshape(yc,m+2,n+2)'; + u = ue(X,Y); + + l = u(:,1); r = u(:,end); + b = u(1,:)'; t = u(end,:)'; + v = {l(2:end-1);r(2:end-1);b;t}; + B = zeros(numCenters,1); + [curL,B] = addScalarBC2D(curL,B,k(i),m,dx,n,dy,dc,nc,v); + [newL,B] = addScalarBC2D(newL,B,k(i),m,dx,n,dy,dc,nc,v); + + curU = curL \ B; + newU = newL \ B; + + curU = reshape(curU,m+2,n+2)'; + newU = reshape(newU,m+2,n+2)'; + + % Maximum Absolute Error + errorCur(i) = max(max(abs(u-curU))); + errorNew(i) = max(max(abs(u-newU))); + +end + +% Plot results +figure +hold on +plot(k, errorCur, 'LineWidth', 2) +plot(k, errorNew, 'LineWidth', 2) +hold off +yscale log +xlim([k(1), k(end)]) +xlabel('k') +ylabel('Maximum Absolute Error') +legend('Current Implementation', 'New Implementation','Location','southwest') +grid on +title("Maximum Absolute Error vs k") +subtitle("∆u = 0, 1 < R < 2, 0 < theta < π, m = 100, n = 25") + + +%% Fixed k, varying mesh Comparison +k = 4; +minCells = 2*k+1; +m = [4 8 16 32 64] * minCells; +n = [1 2 4 8 16] * minCells; +numCenters = (m+2) .* (n+2); +dx = pi ./ m; +dy = 1 ./ n; + +errorCur = zeros(size(numCenters)); +errorNew = zeros(size(numCenters)); + +dc = [1;1;1;1]; nc = [0;0;0;0]; + +for i = 1:numel(numCenters) + + fprintf("Starting m = %d, n = %d\n", m(i),n(i)) + % Nodes for current operators + rsN = 1:dy(i):2; tsN = pi:-dx(i):0; + [TSN,RSN] = meshgrid(tsN,rsN); + xn = RSN .* cos(TSN); + yn = RSN .* sin(TSN); + + % Centers for new operators + rsC = [1 1+dy(i)/2:dy(i):2-dy(i)/2 2]; + tsC = [pi pi-dx(i)/2:-dx(i):dx(i)/2 0]; + [TSC,RSC] = meshgrid(tsC,rsC); + xc = RSC .* cos(TSC); xc = reshape(xc',[],1); + yc = RSC .* sin(TSC); yc = reshape(yc',[],1); + + % Build operators + curG = grad2DCurv(k,xn,yn); + curD = div2DCurv(k,xn,yn); + curL = curD * curG; + + newG = grad2DCurv(k,xc,yc,m(i),dx(i),n(i),dy(i),dc,nc); + newD = div2DCurv(k,xc,yc,m(i),dx(i),n(i),dy(i),dc,nc); + newL = newD * newG; + + % Boundary Conditions + X = reshape(xc,m(i)+2,n(i)+2)'; % Reshape for plotting and easier BC + Y = reshape(yc,m(i)+2,n(i)+2)'; + u = ue(X,Y); + + l = u(:,1); r = u(:,end); + b = u(1,:)'; t = u(end,:)'; + v = {l(2:end-1);r(2:end-1);b;t}; + B = zeros(numCenters(i),1); + [curL,B] = addScalarBC2D(curL,B,k,m(i),dx(i),n(i),dy(i),dc,nc,v); + [newL,B] = addScalarBC2D(newL,B,k,m(i),dx(i),n(i),dy(i),dc,nc,v); + + curU = curL \ B; + newU = newL \ B; + + curU = reshape(curU,m(i)+2,n(i)+2)'; + newU = reshape(newU,m(i)+2,n(i)+2)'; + + % Maximum Absolute Error + errorCur(i) = max(max(abs(u-curU))); + errorNew(i) = max(max(abs(u-newU))); + +end + +% Plot results +h = numCenters.^-0.5; +refCur = errorCur(1) * (h / h(1)).^k; +refNew = errorNew(1) * (h / h(1)).^k; +figure +hold on +plot(h, errorCur, 'LineWidth', 2) +plot(h, errorNew, 'LineWidth', 2) +plot(h, refCur, 'LineWidth', 1.5, 'Color', 'k', 'LineStyle', '--') +plot(h, refNew, 'LineWidth', 1.5, 'Color', 'k', 'LineStyle', '--') +hold off +xscale log +yscale log +xlim([min(h) max(h)]) +xlabel('1 / sqrt(Number of Centers)') +ylabel('Maximum Absolute Error') +legend('Current Implementation', 'New Implementation',"k = " + k + " Reference Line",'Location','southeast') +grid on +title("Maximum Absolute Error vs 1 / sqrt(Number of Centers)") +subtitle("∆u = 0, 1 < R < 2, 0 < theta < π, k = " + k) \ No newline at end of file From 1b44e1b5c7b738d3ef99560d1a2bceeda09b2e9b Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Thu, 9 Apr 2026 16:12:56 -0700 Subject: [PATCH 06/10] jacobian From df5143259250d13383ce2b160d2cf1328e2da99a Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Sun, 26 Apr 2026 19:08:38 -0700 Subject: [PATCH 07/10] 3d and some tweaks --- .../matlab_octave/elliptic2DCurvPeriodic.m | 93 +++++++++ src/matlab_octave/div2DCurv.m | 13 +- src/matlab_octave/div2DCurvLegacy.m | 5 + src/matlab_octave/div3DCurv.m | 174 +++++++++++------ src/matlab_octave/div3DCurvLegacy.m | 75 ++++++++ src/matlab_octave/grad2DCurvLegacy.m | 5 + src/matlab_octave/grad3DCurv.m | 179 ++++++++++-------- src/matlab_octave/grad3DCurvLegacy.m | 97 ++++++++++ src/matlab_octave/jacobian2DLegacy.m | 5 + src/matlab_octave/jacobian3D.m | 148 +++++++++++---- src/matlab_octave/jacobian3DLegacy.m | 55 ++++++ 11 files changed, 661 insertions(+), 188 deletions(-) create mode 100644 examples/matlab_octave/elliptic2DCurvPeriodic.m create mode 100644 src/matlab_octave/div3DCurvLegacy.m create mode 100644 src/matlab_octave/grad3DCurvLegacy.m create mode 100644 src/matlab_octave/jacobian3DLegacy.m diff --git a/examples/matlab_octave/elliptic2DCurvPeriodic.m b/examples/matlab_octave/elliptic2DCurvPeriodic.m new file mode 100644 index 000000000..2060a473f --- /dev/null +++ b/examples/matlab_octave/elliptic2DCurvPeriodic.m @@ -0,0 +1,93 @@ +clear; clc; close all; +addpath("../../src/matlab_octave/") +% +% ∆u = 0 +% +% 1 < r < 2 +% 0 < θ < 2 pi +% +% x = r cos(θ) +% y = r sin(θ) +% +% Dirichlet Boundary Conditions +% +% Exact solution +% u(x, y) = x^2 - y^2 +% + +% Parameters +k = 4; +m = 100; +n = 25; +dx = 2 * pi / m; +dy = 1 / n; +dc = [0;0;1;1]; +nc = [0;0;0;0]; + +% Exact solution +ue = @(X, Y) X.^2 - Y.^2; + +% Create Mesh +rs = [1 (1 + dy / 2) : dy : (2 - dy / 2) 2]; +ts = 0 : dx : (2*pi - dx); +[T,R] = meshgrid(ts, rs); +xc = R .* cos(T); +yc = R .* sin(T); +xc = reshape(xc', [], 1); +yc = reshape(yc', [], 1); + +% Build Operators +G = grad2DCurv(k, xc, yc, m, dx, n, dy, dc, nc); +D = div2DCurv(k, xc, yc, m, dx, n, dy, dc, nc); +L = D * G; + +% Boundary Conditions +X = reshape(xc, m, n + 2)'; % Reshape for plotting and easier BC +Y = reshape(yc, m, n + 2)'; + +u = ue(X, Y); + +l = 0; % Left and right boundaries don't exist +r = 0; +b = u(1, :)'; +t = u(end, :)'; +v = {l; r; b; t}; +B = zeros(size(xc)); +[L0, B0] = addScalarBC2D(L, B, k, m, dx, n, dy, dc, nc, v); + +ua = L0 \ B0; +ua = reshape(ua, m, n + 2)'; + +% Plot results +% Join left and right edges for nicer plotting +X = [X X(:,1)]; +Y = [Y Y(:,1)]; +u = [u u(:,1)]; +ua = [ua ua(:,1)]; + +figure +surf(X, Y, ua, "EdgeColor", "none") +title("Approximate Solution") +xlabel("X") +ylabel("Y") +axis equal +view([0 90]) +cb = colorbar; +cb.Label.String = "u(X,Y)"; +colormap("jet") + +figure +surf(X, Y, u, "EdgeColor", "none") +title("Exact Solution") +xlabel("X") +ylabel("Y") +axis equal +view([0 90]) +cb = colorbar; +cb.Label.String = "u(X,Y)"; +colormap("jet") + +max_err = max(max(abs(ua - u))); +disp("Maximum Absolute Error: " + max_err) +rel_err = 100 * max_err / (max(max(u)) - min(min(u))); +disp("Maximum Relative Error: " + rel_err) \ No newline at end of file diff --git a/src/matlab_octave/div2DCurv.m b/src/matlab_octave/div2DCurv.m index 805b0b0bb..e3aaab8ab 100644 --- a/src/matlab_octave/div2DCurv.m +++ b/src/matlab_octave/div2DCurv.m @@ -50,22 +50,30 @@ IFCx = interpolFacesToCentersG1DPeriodic(k,m); ICFx = interpolCentersToFacesD1DPeriodic(k,m); Im = speye(m); + Bx = sparse(m, m); else De = divNonPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1D(k,m); ICFx = interpolCentersToFacesD1D(k,m); Im = speye(m+2); + Bx = sparse(m + 2, m + 1); + Bx(1, 1) = 1; + Bx(end, end) = 1; end if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) Dn = divPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1DPeriodic(k,n); ICFy = interpolCentersToFacesD1DPeriodic(k,n); In = speye(n); + By = sparse(n, n); else Dn = divNonPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1D(k,n); ICFy = interpolCentersToFacesD1D(k,n); In = speye(n+2); + By = sparse(n + 2, n + 1); + By(1, 1) = 1; + By(end, end) = 1; end % Make De and Dn act on and output to the centers @@ -85,8 +93,9 @@ D = [Dx Dy]; - % Ensure no output on boundary -- Probably a cleaner way to do this - bdry = find(sum(spones(div2D(2,m,1,n,1,dc,nc)), 2) == 0); + % Ensure no output on boundary + B = [kron(By, Im) kron(In, Bx)]; + bdry = find(sum(B, 2) ~= 0); D(bdry,:) = sparse(size(bdry,1),size(D,2)); end diff --git a/src/matlab_octave/div2DCurvLegacy.m b/src/matlab_octave/div2DCurvLegacy.m index f0026a84a..c90cab807 100644 --- a/src/matlab_octave/div2DCurvLegacy.m +++ b/src/matlab_octave/div2DCurvLegacy.m @@ -1,4 +1,9 @@ function D = div2DCurvLegacy(k, X, Y) +% +% ---------------------------------------------------------------------------- +% !!! WARNING: DEPRICATED BY div2DCurv.m !!! +% ---------------------------------------------------------------------------- +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/div3DCurv.m b/src/matlab_octave/div3DCurv.m index ec260fc9d..58727f01b 100644 --- a/src/matlab_octave/div3DCurv.m +++ b/src/matlab_octave/div3DCurv.m @@ -1,11 +1,33 @@ -function D = div3DCurv(k, X, Y, Z) +function D = div3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % PURPOSE -% Returns a 3D curvilinear mimetic divergence +% Returns a 3D curvilinear mimetic divergence operator % % DESCRIPTION +% Returns: +% D : 3D curvilinear mimetic divergence operator. If optional +% arguments are specified, it acts on the extended faces +% (normal faces plus boundaries) +% +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% Y : y-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% (optional) m : Number of cells in xi direction +% (optional) dx : Step size in xi direction +% (optional) n : Number of cells in eta direction +% (optional) dy : Step size in eta direction +% (optional) o : Number of cells in kappa direction +% (optional) dz : Step size in kappa direction +% (optional) dc : a0 (6x1 vector for left, right, bottom, top +% boundaries, resp.) +% (optional) nc : b0 (6x1 vector for left, right, bottom, top +% boundaries, resp.) % % SYNTAX % D = div3DCurv(k, X, Y, Z) +% D = div3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -13,65 +35,93 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - % Get the determinant of the jacobian and the metrics - [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z); - - % Dimensions of nodal grid - [n, m, o] = size(X); - - % Make them volumes so they can be interpolated - J = permute(reshape(J, m, n, o), [2, 1, 3]); - A = permute(reshape(Yn.*Zc-Zn.*Yc, m, n, o), [2, 1, 3]); - B = permute(reshape(Zn.*Xc-Xn.*Zc, m, n, o), [2, 1, 3]); - C = permute(reshape(Xn.*Yc-Yn.*Xc, m, n, o), [2, 1, 3]); - D = permute(reshape(Ze.*Yc-Ye.*Zc, m, n, o), [2, 1, 3]); - E = permute(reshape(Xe.*Zc-Ze.*Xc, m, n, o), [2, 1, 3]); - F = permute(reshape(Ye.*Xc-Xe.*Yc, m, n, o), [2, 1, 3]); - G = permute(reshape(Ye.*Zn-Ze.*Yn, m, n, o), [2, 1, 3]); - H = permute(reshape(Ze.*Xn-Xe.*Zn, m, n, o), [2, 1, 3]); - I = permute(reshape(Xe.*Yn-Ye.*Xn, m, n, o), [2, 1, 3]); - - % Logical grids - [Xl, Yl, Zl] = meshgrid(1:m, 1:n, 1:o); - % Staggered logical grid - [Xs, Ys, Zs] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n], [1 1.5 : 1 : o-0.5 o]); - - % Interpolate the metrics on the logical grid for positions (Xs, Ys, Zs) - J = interp3(Xl, Yl, Zl, J, Xs, Ys, Zs); - A = interp3(Xl, Yl, Zl, A, Xs, Ys, Zs); - B = interp3(Xl, Yl, Zl, B, Xs, Ys, Zs); - C = interp3(Xl, Yl, Zl, C, Xs, Ys, Zs); - D = interp3(Xl, Yl, Zl, D, Xs, Ys, Zs); - E = interp3(Xl, Yl, Zl, E, Xs, Ys, Zs); - F = interp3(Xl, Yl, Zl, F, Xs, Ys, Zs); - G = interp3(Xl, Yl, Zl, G, Xs, Ys, Zs); - H = interp3(Xl, Yl, Zl, H, Xs, Ys, Zs); - I = interp3(Xl, Yl, Zl, I, Xs, Ys, Zs); - - % Convert metrics to diagonal matrices so they can be multiplied by the - % logical operators - J = spdiags(1./reshape(permute(J, [2, 1, 3]), [], 1), 0, numel(J), numel(J)); - A = spdiags(reshape(permute(A, [2, 1, 3]), [], 1), 0, numel(A), numel(A)); - B = spdiags(reshape(permute(B, [2, 1, 3]), [], 1), 0, numel(B), numel(B)); - C = spdiags(reshape(permute(C, [2, 1, 3]), [], 1), 0, numel(C), numel(C)); - D = spdiags(reshape(permute(D, [2, 1, 3]), [], 1), 0, numel(D), numel(D)); - E = spdiags(reshape(permute(E, [2, 1, 3]), [], 1), 0, numel(E), numel(E)); - F = spdiags(reshape(permute(F, [2, 1, 3]), [], 1), 0, numel(F), numel(F)); - G = spdiags(reshape(permute(G, [2, 1, 3]), [], 1), 0, numel(G), numel(G)); - H = spdiags(reshape(permute(H, [2, 1, 3]), [], 1), 0, numel(H), numel(H)); - I = spdiags(reshape(permute(I, [2, 1, 3]), [], 1), 0, numel(I), numel(I)); - - % Construct 3D uniform mimetic divergence operator (d/de + d/dn + d/dc) - Div = div3D(k, m-1, 1, n-1, 1, o-1, 1); - De = Div(:, 1:m*(n-1)*(o-1)); - Dn = Div(:, m*(n-1)*(o-1)+1:m*(n-1)*(o-1)+(m-1)*n*(o-1)); - Dc = Div(:, m*(n-1)*(o-1)+(m-1)*n*(o-1)+1:end); - - % Apply transformation - Dx = J*(A*De+D*DI3(m-1, n-1, o-1, 'Dn')+G*DI3(m-1, n-1, o-1, 'Dc')); - Dy = J*(B*DI3(m-1, n-1, o-1, 'De')+E*Dn+H*DI3(m-1, n-1, o-1, 'Dcc')); - Dz = J*(C*DI3(m-1, n-1, o-1, 'Dee')+F*DI3(m-1, n-1, o-1, 'Dnn')+I*Dc); - - % Final 3D curvilinear mimetic divergence operator (d/dx + d/dy + d/dz) + if nargin ~= 4 && nargin ~= 12 + error("div2DCurv:InvalidNumArgs", ... + "div2DCurv expects 4 or 12 arguments") + elseif nargin == 4 + D = div3DCurvLegacy(k, X, Y, Z); + return; + end + + assert(all(size(dc) == [6 1]), "dc is a 6x1 vector") + assert(all(size(nc) == [6 1]), "nc is a 6x1 vector") + + % Periodic Handling + if isempty(find(dc(1:2).^2 + nc(1:2).^2,1)) + De = divPeriodic(k,m,dx); + IFCx = interpolFacesToCentersG1DPeriodic(k,m); + ICFx = interpolCentersToFacesD1DPeriodic(k,m); + Im = speye(m); + Bx = sparse(m, m); + else + De = divNonPeriodic(k,m,dx); + IFCx = interpolFacesToCentersG1D(k,m); + ICFx = interpolCentersToFacesD1D(k,m); + Im = speye(m+2); + Bx = sparse(m + 2, m + 1); + Bx(1, 1) = 1; + Bx(end, end) = 1; + end + if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) + Dn = divPeriodic(k,n,dy); + IFCy = interpolFacesToCentersG1DPeriodic(k,n); + ICFy = interpolCentersToFacesD1DPeriodic(k,n); + In = speye(n); + By = sparse(n, n); + else + Dn = divNonPeriodic(k,n,dy); + IFCy = interpolFacesToCentersG1D(k,n); + ICFy = interpolCentersToFacesD1D(k,n); + In = speye(n+2); + By = sparse(n + 2, n + 1); + By(1, 1) = 1; + By(end, end) = 1; + end + if isempty(find(dc(5:6).^2 + nc(5:6).^2, 1)) + Dk = divPeriodic(k, o, dz); + IFCz = interpolFacesToCentersG1DPeriodic(k, o); + ICFz = interpolCentersToFacesD1DPeriodic(k, o); + Io = speye(o); + Bz = sparse(o, o); + else + Dk = divNonPeriodic(k, o, dz); + IFCz = interpolFacesToCentersG1D(k, o); + ICFz = interpolCentersToFacesD1D(k, o); + Io = speye(o + 2); + Bz = sparse(o + 2, o + 1); + Bz(1, 1) = 1; + Bz(end, end) = 1; + end + + % Make De, Dn, and Dk act on and output to the centers + % This allows them to be added together + De = kron(kron(Io, In), De) * kron(kron(Io, In), ICFx); + Dn = kron(kron(Io, Dn), Im) * kron(kron(Io, ICFy), Im); + Dk = kron(kron(Dk, In), Im) * kron(kron(ICFz, In), Im); + + % Apply metrics + [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc); + + Dx = ((Yn .* Zk - Zn .* Yk) ./ J) .* De ... + + ((Ze .* Yk - Ye .* Zk) ./ J) .* Dn ... + + ((Ye .* Zn - Ze .* Yn) ./ J) .* Dk; + Dy = ((Zn .* Xk - Xn .* Zk) ./ J) .* De ... + + ((Xe .* Zk - Ze .* Xk) ./ J) .* Dn ... + + ((Ze .* Xn - Xe .* Zn) ./ J) .* Dk; + Dz = ((Xn .* Yk - Yn .* Xk) ./ J) .* De ... + + ((Ye .* Xk - Xe .* Yk) ./ J) .* Dn ... + + ((Xe .* Yn - Ye .* Xn) ./ J) .* Dk; + + % Now have them act on the extended faces + Dx = Dx * kron(kron(Io, In), IFCx); + Dy = Dy * kron(kron(Io, IFCy), Im); + Dz = Dz * kron(kron(IFCz, In), Im); + D = [Dx Dy Dz]; + + % Ensure no output on boundary + B = [kron(kron(Io, In), Bx) kron(kron(Io, By), Im) kron(kron(Bz, In), Im)]; + bdry = find(sum(B, 2) ~= 0); + D(bdry, :) = sparse(size(bdry, 1), size(D, 2)); + end diff --git a/src/matlab_octave/div3DCurvLegacy.m b/src/matlab_octave/div3DCurvLegacy.m new file mode 100644 index 000000000..6acc77be9 --- /dev/null +++ b/src/matlab_octave/div3DCurvLegacy.m @@ -0,0 +1,75 @@ +function D = div3DCurvLegacy(k, X, Y, Z) +% +% ---------------------------------------------------------------------------- +% !!! WARNING: DEPRICATED BY div2DCurv.m !!! +% ---------------------------------------------------------------------------- +% +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- +% Returns a 3D curvilinear mimetic divergence + + % Get the determinant of the jacobian and the metrics + [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z); + + % Dimensions of nodal grid + [n, m, o] = size(X); + + % Make them volumes so they can be interpolated + J = permute(reshape(J, m, n, o), [2, 1, 3]); + A = permute(reshape(Yn.*Zc-Zn.*Yc, m, n, o), [2, 1, 3]); + B = permute(reshape(Zn.*Xc-Xn.*Zc, m, n, o), [2, 1, 3]); + C = permute(reshape(Xn.*Yc-Yn.*Xc, m, n, o), [2, 1, 3]); + D = permute(reshape(Ze.*Yc-Ye.*Zc, m, n, o), [2, 1, 3]); + E = permute(reshape(Xe.*Zc-Ze.*Xc, m, n, o), [2, 1, 3]); + F = permute(reshape(Ye.*Xc-Xe.*Yc, m, n, o), [2, 1, 3]); + G = permute(reshape(Ye.*Zn-Ze.*Yn, m, n, o), [2, 1, 3]); + H = permute(reshape(Ze.*Xn-Xe.*Zn, m, n, o), [2, 1, 3]); + I = permute(reshape(Xe.*Yn-Ye.*Xn, m, n, o), [2, 1, 3]); + + % Logical grids + [Xl, Yl, Zl] = meshgrid(1:m, 1:n, 1:o); + % Staggered logical grid + [Xs, Ys, Zs] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n], [1 1.5 : 1 : o-0.5 o]); + + % Interpolate the metrics on the logical grid for positions (Xs, Ys, Zs) + J = interp3(Xl, Yl, Zl, J, Xs, Ys, Zs); + A = interp3(Xl, Yl, Zl, A, Xs, Ys, Zs); + B = interp3(Xl, Yl, Zl, B, Xs, Ys, Zs); + C = interp3(Xl, Yl, Zl, C, Xs, Ys, Zs); + D = interp3(Xl, Yl, Zl, D, Xs, Ys, Zs); + E = interp3(Xl, Yl, Zl, E, Xs, Ys, Zs); + F = interp3(Xl, Yl, Zl, F, Xs, Ys, Zs); + G = interp3(Xl, Yl, Zl, G, Xs, Ys, Zs); + H = interp3(Xl, Yl, Zl, H, Xs, Ys, Zs); + I = interp3(Xl, Yl, Zl, I, Xs, Ys, Zs); + + % Convert metrics to diagonal matrices so they can be multiplied by the + % logical operators + J = spdiags(1./reshape(permute(J, [2, 1, 3]), [], 1), 0, numel(J), numel(J)); + A = spdiags(reshape(permute(A, [2, 1, 3]), [], 1), 0, numel(A), numel(A)); + B = spdiags(reshape(permute(B, [2, 1, 3]), [], 1), 0, numel(B), numel(B)); + C = spdiags(reshape(permute(C, [2, 1, 3]), [], 1), 0, numel(C), numel(C)); + D = spdiags(reshape(permute(D, [2, 1, 3]), [], 1), 0, numel(D), numel(D)); + E = spdiags(reshape(permute(E, [2, 1, 3]), [], 1), 0, numel(E), numel(E)); + F = spdiags(reshape(permute(F, [2, 1, 3]), [], 1), 0, numel(F), numel(F)); + G = spdiags(reshape(permute(G, [2, 1, 3]), [], 1), 0, numel(G), numel(G)); + H = spdiags(reshape(permute(H, [2, 1, 3]), [], 1), 0, numel(H), numel(H)); + I = spdiags(reshape(permute(I, [2, 1, 3]), [], 1), 0, numel(I), numel(I)); + + % Construct 3D uniform mimetic divergence operator (d/de + d/dn + d/dc) + Div = div3D(k, m-1, 1, n-1, 1, o-1, 1); + De = Div(:, 1:m*(n-1)*(o-1)); + Dn = Div(:, m*(n-1)*(o-1)+1:m*(n-1)*(o-1)+(m-1)*n*(o-1)); + Dc = Div(:, m*(n-1)*(o-1)+(m-1)*n*(o-1)+1:end); + + % Apply transformation + Dx = J*(A*De+D*DI3(m-1, n-1, o-1, 'Dn')+G*DI3(m-1, n-1, o-1, 'Dc')); + Dy = J*(B*DI3(m-1, n-1, o-1, 'De')+E*Dn+H*DI3(m-1, n-1, o-1, 'Dcc')); + Dz = J*(C*DI3(m-1, n-1, o-1, 'Dee')+F*DI3(m-1, n-1, o-1, 'Dnn')+I*Dc); + + % Final 3D curvilinear mimetic divergence operator (d/dx + d/dy + d/dz) + D = [Dx Dy Dz]; +end diff --git a/src/matlab_octave/grad2DCurvLegacy.m b/src/matlab_octave/grad2DCurvLegacy.m index 9b0b1f359..c30cfd744 100644 --- a/src/matlab_octave/grad2DCurvLegacy.m +++ b/src/matlab_octave/grad2DCurvLegacy.m @@ -1,4 +1,9 @@ function G = grad2DCurvLegacy(k, X, Y) +% +% ---------------------------------------------------------------------------- +% !!! WARNING: DEPRICATED BY grad2DCurv.m !!! +% ---------------------------------------------------------------------------- +% % Returns a 2D curvilinear mimetic gradient % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later diff --git a/src/matlab_octave/grad3DCurv.m b/src/matlab_octave/grad3DCurv.m index 7540f7c6c..c463d467a 100644 --- a/src/matlab_octave/grad3DCurv.m +++ b/src/matlab_octave/grad3DCurv.m @@ -1,11 +1,33 @@ -function G = grad3DCurv(k, X, Y, Z) +function G = grad3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % PURPOSE -% Returns a 3D curvilinear mimetic gradient +% Returns a 3D curvilinear mimetic gradient operator % % DESCRIPTION +% Returns: +% G : 3D curvilinear mimetic gradient operator. If optional +% arguments are specified, it outputs to the extended +% faces (normal faces plus boundaries) +% +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% Y : y-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes +% (optional) m : Number of cells in xi direction +% (optional) dx : Step size in xi direction +% (optional) n : Number of cells in eta direction +% (optional) dy : Step size in eta direction +% (optional) o : Number of cells in kappa direction +% (optional) dz : Step size in kappa direction +% (optional) dc : a0 (6x1 vector for left, right, bottom, top +% boundaries, resp.) +% (optional) nc : b0 (6x1 vector for left, right, bottom, top +% boundaries, resp.) % % SYNTAX % G = grad3DCurv(k, X, Y, Z) +% G = grad3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -13,87 +35,76 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - % Get the determinant of the jacobian and the metrics - [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z); - - % Dimensions of nodal grid - [n, m, o] = size(X); - - % Make them volumes so they can be interpolated - J = permute(reshape(J, m, n, o), [2, 1, 3]); - A = permute(reshape(Yn.*Zc-Zn.*Yc, m, n, o), [2, 1, 3]); - B = permute(reshape(Zn.*Xc-Xn.*Zc, m, n, o), [2, 1, 3]); - C = permute(reshape(Xn.*Yc-Yn.*Xc, m, n, o), [2, 1, 3]); - D = permute(reshape(Ze.*Yc-Ye.*Zc, m, n, o), [2, 1, 3]); - E = permute(reshape(Xe.*Zc-Ze.*Xc, m, n, o), [2, 1, 3]); - F = permute(reshape(Ye.*Xc-Xe.*Yc, m, n, o), [2, 1, 3]); - G = permute(reshape(Ye.*Zn-Ze.*Yn, m, n, o), [2, 1, 3]); - H = permute(reshape(Ze.*Xn-Xe.*Zn, m, n, o), [2, 1, 3]); - I = permute(reshape(Xe.*Yn-Ye.*Xn, m, n, o), [2, 1, 3]); - - % Logical grids - [Xl, Yl, Zl] = meshgrid(1:m, 1:n, 1:o); - - % Interpolate the metrics on the logical grid for positions u, v and w - Xl_ = (Xl(1:end-1, :, :)+Xl(2:end, :, :))/2; - Xl_ = (Xl_(:, :, 1:end-1)+Xl_(:, :, 2:end))/2; - Yl_ = (Yl(1:end-1, :, :)+Yl(2:end, :, :))/2; - Yl_ = (Yl_(:, :, 1:end-1)+Yl_(:, :, 2:end))/2; - Zl_ = (Zl(1:end-1, :, :)+Zl(2:end, :, :))/2; - Zl_ = (Zl_(:, :, 1:end-1)+Zl_(:, :, 2:end))/2; - Ju = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_); - A = interp3(Xl, Yl, Zl, A, Xl_, Yl_, Zl_); - D = interp3(Xl, Yl, Zl, D, Xl_, Yl_, Zl_); - G = interp3(Xl, Yl, Zl, G, Xl_, Yl_, Zl_); - - Xl_ = (Xl(:, 1:end-1, :)+Xl(:, 2:end, :))/2; - Xl_ = (Xl_(:, :, 1:end-1)+Xl_(:, :, 2:end))/2; - Yl_ = (Yl(:, 1:end-1, :)+Yl(:, 2:end, :))/2; - Yl_ = (Yl_(:, :, 1:end-1)+Yl_(:, :, 2:end))/2; - Zl_ = (Zl(:, 1:end-1, :)+Zl(:, 2:end, :))/2; - Zl_ = (Zl_(:, :, 1:end-1)+Zl_(:, :, 2:end))/2; - Jv = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_); - B = interp3(Xl, Yl, Zl, B, Xl_, Yl_, Zl_); - E = interp3(Xl, Yl, Zl, E, Xl_, Yl_, Zl_); - H = interp3(Xl, Yl, Zl, H, Xl_, Yl_, Zl_); - - Xl_ = (Xl(1:end-1, :, :)+Xl(2:end, :, :))/2; - Xl_ = (Xl_(:, 1:end-1, :)+Xl_(:, 2:end, :))/2; - Yl_ = (Yl(1:end-1, :, :)+Yl(2:end, :, :))/2; - Yl_ = (Yl_(:, 1:end-1, :)+Yl_(:, 2:end, :))/2; - Zl_ = (Zl(1:end-1, :, :)+Zl(2:end, :, :))/2; - Zl_ = (Zl_(:, 1:end-1, :)+Zl_(:, 2:end, :))/2; - Jw = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_); - C = interp3(Xl, Yl, Zl, C, Xl_, Yl_, Zl_); - F = interp3(Xl, Yl, Zl, F, Xl_, Yl_, Zl_); - I = interp3(Xl, Yl, Zl, I, Xl_, Yl_, Zl_); - - % Convert metrics to diagonal matrices so they can be multiplied by the - % logical operators - Ju = spdiags(1./reshape(permute(Ju, [2, 1, 3]), [], 1), 0, numel(Ju), numel(Ju)); - Jv = spdiags(1./reshape(permute(Jv, [2, 1, 3]), [], 1), 0, numel(Jv), numel(Jv)); - Jw = spdiags(1./reshape(permute(Jw, [2, 1, 3]), [], 1), 0, numel(Jw), numel(Jw)); - A = spdiags(reshape(permute(A, [2, 1, 3]), [], 1), 0, numel(A), numel(A)); - B = spdiags(reshape(permute(B, [2, 1, 3]), [], 1), 0, numel(B), numel(B)); - C = spdiags(reshape(permute(C, [2, 1, 3]), [], 1), 0, numel(C), numel(C)); - D = spdiags(reshape(permute(D, [2, 1, 3]), [], 1), 0, numel(D), numel(D)); - E = spdiags(reshape(permute(E, [2, 1, 3]), [], 1), 0, numel(E), numel(E)); - F = spdiags(reshape(permute(F, [2, 1, 3]), [], 1), 0, numel(F), numel(F)); - G = spdiags(reshape(permute(G, [2, 1, 3]), [], 1), 0, numel(G), numel(G)); - H = spdiags(reshape(permute(H, [2, 1, 3]), [], 1), 0, numel(H), numel(H)); - I = spdiags(reshape(permute(I, [2, 1, 3]), [], 1), 0, numel(I), numel(I)); - - % Construct 3D uniform mimetic gradient operator (d/de, d/dn, d/dc) - Grad = grad3D(k, m-1, 1, n-1, 1, o-1, 1); - Ge = Grad(1:m*(n-1)*(o-1), :); - Gn = Grad(m*(n-1)*(o-1)+1:m*(n-1)*(o-1)+(m-1)*n*(o-1), :); - Gc = Grad(m*(n-1)*(o-1)+(m-1)*n*(o-1)+1:end, :); - - % Apply transformation - Gx = Ju*(A*Ge+D*GI13(Gn, m-1, n-1, o-1, 'Gn')+G*GI13(Gc, m-1, n-1, o-1, 'Gc')); - Gy = Jv*(B*GI13(Ge, m-1, n-1, o-1, 'Ge')+E*Gn+H*GI13(Gc, m-1, n-1, o-1, 'Gcy')); - Gz = Jw*(C*GI13(Ge, m-1, n-1, o-1, 'Gee')+F*GI13(Gn, m-1, n-1, o-1, 'Gnn')+I*Gc); - - % Final 3D curvilinear mimetic gradient operator (d/dx, d/dy, d/dz) + if nargin ~= 4 && nargin ~= 12 + error("grad3DCurv:InvalidNumArgs", ... + "grad3DCurv expects 4 or 12 arguments") + elseif nargin == 4 + G = grad3DCurvLegacy(k, X, Y, Z); + return; + end + + assert(all(size(dc) == [6 1]), "dc is a 6x1 vector") + assert(all(size(nc) == [6 1]), "nc is a 6x1 vector") + + % Periodic Handling + if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) + Ge = gradPeriodic(k,m,dx); + ICFx = interpolCentersToFacesD1DPeriodic(k,m); + IFCx = interpolFacesToCentersG1DPeriodic(k,m); + Im = speye(m); + else + Ge = gradNonPeriodic(k,m,dx); + ICFx = interpolCentersToFacesD1D(k,m); + IFCx = interpolFacesToCentersG1D(k,m); + Im = speye(m+2); + end + if isempty(find(dc(3:4).^2 + nc(3:4).^2, 1)) + Gn = gradPeriodic(k,n,dy); + ICFy = interpolCentersToFacesD1DPeriodic(k,n); + IFCy = interpolFacesToCentersG1DPeriodic(k,n); + In = speye(n); + else + Gn = gradNonPeriodic(k,n,dy); + ICFy = interpolCentersToFacesD1D(k,n); + IFCy = interpolFacesToCentersG1D(k,n); + In = speye(n+2); + end + if isempty(find(dc(5:6).^2 + nc(5:6).^2, 1)) + Gk = gradPeriodic(k, o, dz); + ICFz = interpolCentersToFacesD1DPeriodic(k,o); + IFCz = interpolFacesToCentersG1DPeriodic(k,o); + Io = speye(o); + else + Gk = gradNonPeriodic(k,o,dz); + ICFz = interpolCentersToFacesD1D(k,o); + IFCz = interpolFacesToCentersG1D(k,o); + Io = speye(o+2); + end + + % Make Ge, Gn, and Gk act on and output to the centers + % This allows them to be added together + Ge = kron(kron(Io, In), IFCx) * kron(kron(Io, In), Ge); + Gn = kron(kron(Io, IFCy), Im) * kron(kron(Io, Gn), Im); + Gk = kron(kron(IFCz, In), Im) * kron(kron(Gk, In), Im); + + % Apply metrics + [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc); + + Gx = ((Yn .* Zk - Zn .* Yk) ./ J) .* Ge ... + + ((Ze .* Yk - Ye .* Zk) ./ J) .* Gn ... + + ((Ye .* Zn - Ze .* Yn) ./ J) .* Gk; + Gy = ((Zn .* Xk - Xn .* Zk) ./ J) .* Ge ... + + ((Xe .* Zk - Ze .* Xk) ./ J) .* Gn ... + + ((Ze .* Xn - Xe .* Zn) ./ J) .* Gk; + Gz = ((Xn .* Yk - Yn .* Xk) ./ J) .* Ge ... + + ((Ye .* Xk - Xe .* Yk) ./ J) .* Gn ... + + ((Xe .* Yn - Ye .* Xn) ./ J) .* Gk; + + % Now have them output to the extended faces + Gx = kron(kron(Io, In), ICFx) * Gx; + Gy = kron(kron(Io, ICFy), Im) * Gy; + Gz = kron(kron(ICFz, In), Im) * Gz; + G = [Gx; Gy; Gz]; + end diff --git a/src/matlab_octave/grad3DCurvLegacy.m b/src/matlab_octave/grad3DCurvLegacy.m new file mode 100644 index 000000000..dad96f4c1 --- /dev/null +++ b/src/matlab_octave/grad3DCurvLegacy.m @@ -0,0 +1,97 @@ +function G = grad3DCurvLegacy(k, X, Y, Z) +% +% ---------------------------------------------------------------------------- +% !!! WARNING: DEPRICATED BY grad2DCurv.m !!! +% ---------------------------------------------------------------------------- +% +% Returns a 3D curvilinear mimetic gradient +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- + + % Get the determinant of the jacobian and the metrics + [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z); + + % Dimensions of nodal grid + [n, m, o] = size(X); + + % Make them volumes so they can be interpolated + J = permute(reshape(J, m, n, o), [2, 1, 3]); + A = permute(reshape(Yn.*Zc-Zn.*Yc, m, n, o), [2, 1, 3]); + B = permute(reshape(Zn.*Xc-Xn.*Zc, m, n, o), [2, 1, 3]); + C = permute(reshape(Xn.*Yc-Yn.*Xc, m, n, o), [2, 1, 3]); + D = permute(reshape(Ze.*Yc-Ye.*Zc, m, n, o), [2, 1, 3]); + E = permute(reshape(Xe.*Zc-Ze.*Xc, m, n, o), [2, 1, 3]); + F = permute(reshape(Ye.*Xc-Xe.*Yc, m, n, o), [2, 1, 3]); + G = permute(reshape(Ye.*Zn-Ze.*Yn, m, n, o), [2, 1, 3]); + H = permute(reshape(Ze.*Xn-Xe.*Zn, m, n, o), [2, 1, 3]); + I = permute(reshape(Xe.*Yn-Ye.*Xn, m, n, o), [2, 1, 3]); + + % Logical grids + [Xl, Yl, Zl] = meshgrid(1:m, 1:n, 1:o); + + % Interpolate the metrics on the logical grid for positions u, v and w + Xl_ = (Xl(1:end-1, :, :)+Xl(2:end, :, :))/2; + Xl_ = (Xl_(:, :, 1:end-1)+Xl_(:, :, 2:end))/2; + Yl_ = (Yl(1:end-1, :, :)+Yl(2:end, :, :))/2; + Yl_ = (Yl_(:, :, 1:end-1)+Yl_(:, :, 2:end))/2; + Zl_ = (Zl(1:end-1, :, :)+Zl(2:end, :, :))/2; + Zl_ = (Zl_(:, :, 1:end-1)+Zl_(:, :, 2:end))/2; + Ju = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_); + A = interp3(Xl, Yl, Zl, A, Xl_, Yl_, Zl_); + D = interp3(Xl, Yl, Zl, D, Xl_, Yl_, Zl_); + G = interp3(Xl, Yl, Zl, G, Xl_, Yl_, Zl_); + + Xl_ = (Xl(:, 1:end-1, :)+Xl(:, 2:end, :))/2; + Xl_ = (Xl_(:, :, 1:end-1)+Xl_(:, :, 2:end))/2; + Yl_ = (Yl(:, 1:end-1, :)+Yl(:, 2:end, :))/2; + Yl_ = (Yl_(:, :, 1:end-1)+Yl_(:, :, 2:end))/2; + Zl_ = (Zl(:, 1:end-1, :)+Zl(:, 2:end, :))/2; + Zl_ = (Zl_(:, :, 1:end-1)+Zl_(:, :, 2:end))/2; + Jv = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_); + B = interp3(Xl, Yl, Zl, B, Xl_, Yl_, Zl_); + E = interp3(Xl, Yl, Zl, E, Xl_, Yl_, Zl_); + H = interp3(Xl, Yl, Zl, H, Xl_, Yl_, Zl_); + + Xl_ = (Xl(1:end-1, :, :)+Xl(2:end, :, :))/2; + Xl_ = (Xl_(:, 1:end-1, :)+Xl_(:, 2:end, :))/2; + Yl_ = (Yl(1:end-1, :, :)+Yl(2:end, :, :))/2; + Yl_ = (Yl_(:, 1:end-1, :)+Yl_(:, 2:end, :))/2; + Zl_ = (Zl(1:end-1, :, :)+Zl(2:end, :, :))/2; + Zl_ = (Zl_(:, 1:end-1, :)+Zl_(:, 2:end, :))/2; + Jw = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_); + C = interp3(Xl, Yl, Zl, C, Xl_, Yl_, Zl_); + F = interp3(Xl, Yl, Zl, F, Xl_, Yl_, Zl_); + I = interp3(Xl, Yl, Zl, I, Xl_, Yl_, Zl_); + + % Convert metrics to diagonal matrices so they can be multiplied by the + % logical operators + Ju = spdiags(1./reshape(permute(Ju, [2, 1, 3]), [], 1), 0, numel(Ju), numel(Ju)); + Jv = spdiags(1./reshape(permute(Jv, [2, 1, 3]), [], 1), 0, numel(Jv), numel(Jv)); + Jw = spdiags(1./reshape(permute(Jw, [2, 1, 3]), [], 1), 0, numel(Jw), numel(Jw)); + A = spdiags(reshape(permute(A, [2, 1, 3]), [], 1), 0, numel(A), numel(A)); + B = spdiags(reshape(permute(B, [2, 1, 3]), [], 1), 0, numel(B), numel(B)); + C = spdiags(reshape(permute(C, [2, 1, 3]), [], 1), 0, numel(C), numel(C)); + D = spdiags(reshape(permute(D, [2, 1, 3]), [], 1), 0, numel(D), numel(D)); + E = spdiags(reshape(permute(E, [2, 1, 3]), [], 1), 0, numel(E), numel(E)); + F = spdiags(reshape(permute(F, [2, 1, 3]), [], 1), 0, numel(F), numel(F)); + G = spdiags(reshape(permute(G, [2, 1, 3]), [], 1), 0, numel(G), numel(G)); + H = spdiags(reshape(permute(H, [2, 1, 3]), [], 1), 0, numel(H), numel(H)); + I = spdiags(reshape(permute(I, [2, 1, 3]), [], 1), 0, numel(I), numel(I)); + + % Construct 3D uniform mimetic gradient operator (d/de, d/dn, d/dc) + Grad = grad3D(k, m-1, 1, n-1, 1, o-1, 1); + Ge = Grad(1:m*(n-1)*(o-1), :); + Gn = Grad(m*(n-1)*(o-1)+1:m*(n-1)*(o-1)+(m-1)*n*(o-1), :); + Gc = Grad(m*(n-1)*(o-1)+(m-1)*n*(o-1)+1:end, :); + + % Apply transformation + Gx = Ju*(A*Ge+D*GI13(Gn, m-1, n-1, o-1, 'Gn')+G*GI13(Gc, m-1, n-1, o-1, 'Gc')); + Gy = Jv*(B*GI13(Ge, m-1, n-1, o-1, 'Ge')+E*Gn+H*GI13(Gc, m-1, n-1, o-1, 'Gcy')); + Gz = Jw*(C*GI13(Ge, m-1, n-1, o-1, 'Gee')+F*GI13(Gn, m-1, n-1, o-1, 'Gnn')+I*Gc); + + % Final 3D curvilinear mimetic gradient operator (d/dx, d/dy, d/dz) + G = [Gx; Gy; Gz]; +end diff --git a/src/matlab_octave/jacobian2DLegacy.m b/src/matlab_octave/jacobian2DLegacy.m index ddafd3959..25bcc7b39 100644 --- a/src/matlab_octave/jacobian2DLegacy.m +++ b/src/matlab_octave/jacobian2DLegacy.m @@ -1,4 +1,9 @@ function [J, Xe, Xn, Ye, Yn] = jacobian2DLegacy(k, X, Y) +% +% ---------------------------------------------------------------------------- +% !!! WARNING: DEPRICATED BY jacobian2D.m !!! +% ---------------------------------------------------------------------------- +% % Returns: % J : Determinant of the Jacobian (XeYn - XnYe) on the nodes % Xe : dx/de metric on the nodes diff --git a/src/matlab_octave/jacobian3D.m b/src/matlab_octave/jacobian3D.m index 4e74f1608..4a85da851 100644 --- a/src/matlab_octave/jacobian3D.m +++ b/src/matlab_octave/jacobian3D.m @@ -1,27 +1,52 @@ -function [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z) +function [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % PURPOSE +% 3D Jacobian metrics for curvilinear operators % % DESCRIPTION % Returns: -% J : Determinant of the Jacobian -% Xe : dx/de metric -% Xn : dx/dn metric -% Xc : dx/dc metric -% Ye : dy/de metric -% Yn : dy/dn metric -% Yc : dy/dc metric -% Ze : dz/de metric -% Zn : dz/dn metric -% Zc : dz/dc metric +% J : Determinant of the Jacobian on the centers if optional +% arguments are specified, else nodes +% Xe : dx/de metric on the centers if optional arguments are +% specified, else nodes +% Xn : dx/dn metric on the centers if optional arguments are +% specified, else nodes +% Xk : dx/dk metric on the centers if optional arguments are +% specified, else nodes +% Ye : dy/de metric on the centers if optional arguments are +% specified, else nodes +% Yn : dy/dn metric on the centers if optional arguments are +% specified, else nodes +% Yk : dy/dk metric on the centers if optional arguments are +% specified, else nodes +% Ze : dz/de metric on the centers if optional arguments are +% specified, else nodes +% Zn : dz/dn metric on the centers if optional arguments are +% specified, else nodes +% Zk : dz/dk metric on the centers if optional arguments are +% specified, else nodes % % Parameters: % k : Order of accuracy -% X : x-coordinates (physical) of meshgrid -% Y : y-coordinates (physical) of meshgrid -% Z : z-coordinates (physical) of meshgrid +% X : x-coordinates (physical) of meshgrid centers if optional +% arguments are specified, else nodes +% Y : y-coordinates (physical) of meshgrid centers if optional +% arguments are specified, else nodes +% Z : z-coordinates (physical) of meshgrid centers if optional +% arguments are specified, else nodes +% (optional) m : Number of cells in xi direction +% (optional) dx : Step size in xi direction +% (optional) n : Number of cells in eta direction +% (optional) dy : Step size in eta direction +% (optional) o : Number of cells in kappa direction +% (optional) dz : Step size in kappa direction +% (optional) dc : a0 (6x1 vector for left, right, bottom, top +% boundaries, resp.) +% (optional) nc : b0 (6x1 vector for left, right, bottom, top +% boundaries, resp.) % % SYNTAX -% [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z) +% [J, Xe, Xn, Xv, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z) +% [J, Xe, Xn, Xv, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -29,29 +54,72 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - [n, m, o] = size(X); - - X = reshape(permute(X, [2, 1, 3]), [], 1); - Y = reshape(permute(Y, [2, 1, 3]), [], 1); - Z = reshape(permute(Z, [2, 1, 3]), [], 1); - - N = nodal3D(k, m, 1, n, 1, o, 1); - - X = N*X; - Y = N*Y; - Z = N*Z; - - mno = m*n*o; - - Xe = X(1:mno); - Xn = X(mno+1:2*mno); - Xc = X(2*mno+1:end); - Ye = Y(1:mno); - Yn = Y(mno+1:2*mno); - Yc = Y(2*mno+1:end); - Ze = Z(1:mno); - Zn = Z(mno+1:2*mno); - Zc = Z(2*mno+1:end); - - J = Xe.*(Yn.*Zc-Yc.*Zn)-Ye.*(Xn.*Zc-Xc.*Zn)+Ze.*(Xn.*Yc-Xc.*Yn); + if nargin ~= 4 && nargin ~= 12 + error("jacobian3D:InvalidNumArgs", ... + "jacobian3D expects 4 or 9 arguments") + elseif nargin == 4 + [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3DLegacy(k, X, Y, Z); + return; + end + + assert(all(size(dc) == [6 1]), "dc is a 6x1 vector") + assert(all(size(nc) == [6 1]), "nc is a 6x1 vector") + + % Periodic Handling + if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) + Ge = gradPeriodic(k, m, dx); + IFCx = interpolFacesToCentersG1DPeriodic(k, m); + Im = speye(m); + else + Ge = gradNonPeriodic(k, m, dx); + IFCx = interpolFacesToCentersG1D(k, m); + Im = speye(m + 2); + end + if isempty(find(dc(3:4).^2 + nc(3:4).^2, 1)) + Gn = gradPeriodic(k, n, dy); + IFCy = interpolFacesToCentersG1DPeriodic(k, n); + In = speye(n); + else + Gn = gradNonPeriodic(k, n, dy); + IFCy = interpolFacesToCentersG1D(k, n); + In = speye(n + 2); + end + if isempty(find(dc(5:6).^2 + nc(5:6).^2, 1)) + Gk = gradPeriodic(k, o, dz); + IFCz = interpolFacesToCentersG1DPeriodic(k, o); + Io = speye(o); + else + Gk = gradNonPeriodic(k, o, dz); + IFCz = interpolFacesToCentersG1D(k, o); + Io = speye(o + 2); + end + + numC = numel(X); + + % Get metrics on centers + % We don't augment the identity matrices + % so that we don't lose information + % around the boundaries + Ge = kron(kron(Io, In), IFCx) * kron(kron(Io, In), Ge); + Gn = kron(kron(Io, IFCy), Im) * kron(kron(Io, Gn), Im); + Gk = kron(kron(IFCz, In), Im) * kron(kron(Gk, In), Im); + G = [Ge; Gn; Gk]; + + X = reshape(X', [], 1); + Y = reshape(Y', [], 1); + Z = reshape(Z', [], 1); + metrics = G * [X Y Z]; + + Xe = metrics(1:numC, 1); + Xn = metrics(1+numC:2*numC, 1); + Xk = metrics(1+2*numC:end, 1); + Ye = metrics(1:numC, 2); + Yn = metrics(1+numC:2*numC, 2); + Yk = metrics(1+2*numC:end, 2); + Ze = metrics(1:numC, 3); + Zn = metrics(1+numC:2*numC, 3); + Zk = metrics(1+2*numC:end, 3); + + J = Xe .* (Yn .* Zk - Yk .* Zn) - Xn .* (Ye .* Zk - Yk .* Ze) + Xk .* (Ye .* Zn - Yn .* Ze); + end \ No newline at end of file diff --git a/src/matlab_octave/jacobian3DLegacy.m b/src/matlab_octave/jacobian3DLegacy.m new file mode 100644 index 000000000..562b8392f --- /dev/null +++ b/src/matlab_octave/jacobian3DLegacy.m @@ -0,0 +1,55 @@ +function [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3DLegacy(k, X, Y, Z) +% +% ---------------------------------------------------------------------------- +% !!! WARNING: DEPRICATED BY jacobian3D.m !!! +% ---------------------------------------------------------------------------- +% +% Returns: +% J : Determinant of the Jacobian +% Xe : dx/de metric +% Xn : dx/dn metric +% Xc : dx/dc metric +% Ye : dy/de metric +% Yn : dy/dn metric +% Yc : dy/dc metric +% Ze : dz/de metric +% Zn : dz/dn metric +% Zc : dz/dc metric +% +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid +% Y : y-coordinates (physical) of meshgrid +% Z : z-coordinates (physical) of meshgrid +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- + + [n, m, o] = size(X); + + X = reshape(permute(X, [2, 1, 3]), [], 1); + Y = reshape(permute(Y, [2, 1, 3]), [], 1); + Z = reshape(permute(Z, [2, 1, 3]), [], 1); + + N = nodal3D(k, m, 1, n, 1, o, 1); + + X = N*X; + Y = N*Y; + Z = N*Z; + + mno = m*n*o; + + Xe = X(1:mno); + Xn = X(mno+1:2*mno); + Xc = X(2*mno+1:end); + Ye = Y(1:mno); + Yn = Y(mno+1:2*mno); + Yc = Y(2*mno+1:end); + Ze = Z(1:mno); + Zn = Z(mno+1:2*mno); + Zc = Z(2*mno+1:end); + + J = Xe.*(Yn.*Zc-Yc.*Zn)-Ye.*(Xn.*Zc-Xc.*Zn)+Ze.*(Xn.*Yc-Xc.*Yn); +end \ No newline at end of file From a3ce9dc76e0358d12cbf0687823c0931afbc233c Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Wed, 29 Apr 2026 22:26:26 -0700 Subject: [PATCH 08/10] examples --- examples/matlab_octave/test_div2DCurv.m | 151 +++++++++-------- examples/matlab_octave/test_div3DCurv.m | 137 +++++----------- examples/matlab_octave/test_grad2DCurv.m | 90 +++++++++-- examples/matlab_octave/test_grad3DCurv.m | 152 +++++++++++++----- examples/matlab_octave/test_lap3DCurv.m | 37 +++-- examples/matlab_octave/test_lap3DCurv_case2.m | 44 ++--- src/matlab_octave/jacobian2D.m | 2 - src/matlab_octave/jacobian3D.m | 3 - 8 files changed, 354 insertions(+), 262 deletions(-) diff --git a/examples/matlab_octave/test_div2DCurv.m b/examples/matlab_octave/test_div2DCurv.m index 3ea388cd5..cbd5f01fc 100644 --- a/examples/matlab_octave/test_div2DCurv.m +++ b/examples/matlab_octave/test_div2DCurv.m @@ -1,103 +1,93 @@ % Tests the 2D curvilinear divergence +% on an annulus of 1 < r < 2 clc close all +clear addpath('../../src/matlab_octave') % Parameters k = 2; -m = 20; -n = 30; - -% Grid -r1 = 1; % Inner radius -r2 = 2; % Outer radius -nR = linspace(r1, r2, m) ; -nT = linspace(0, 2*pi, n) ; -[R, T] = meshgrid(nR, nT) ; -% Convert grid to cartesian coordinates -X = R.*cos(T); -Y = R.*sin(T); - -% Test on another grid -% [X, Y] = genCurvGrid(n, m); - -mesh(X, Y, zeros(n, m), 'Marker', '.', 'MarkerSize', 10) -% Az El -view([0 90]) +m = 30; +n = 20; +dx = 1 / m; +dy = 1 / n; +dc = [0; 0; 1; 1]; +nc = [0; 0; 0; 0]; + +% Nodal Grid Generation +r1 = 1; +r2 = 2; +nR = linspace(r1, r2, n + 1); +nT = linspace(2 * pi, 0, m + 1); +nT = nT(1:end-1); % No duplicate points +[T, R] = meshgrid(nT, nR); +X = R .* cos(T); +Y = R .* sin(T); + +xn = reshape(X', [], 1); +yn = reshape(Y', [], 1); + +% Interpolators +NtoC = interpolNodesToCenters2D(k, m, n, dc, nc); +CtoU = interpolCentersToFacesD1DPeriodic(k, m); +CtoU = kron(speye(n + 2), CtoU); % From the centers to the extended faces +CtoV = interpolCentersToFacesD1D(k, n); +CtoV = kron(CtoV, speye(m)); % From the centers to the extended faces + +% Plot Mesh +mesh([X X(:, 1)], [Y Y(:, 1)], zeros(n + 1, m + 1), 'Marker', '.', 'MarkerSize', 10) +view(0, 90) axis equal hold on -[n, m] = size(X); -n = n-1; -m = m-1; - -Ux = (X(1:end-1, :) + X(2:end, :))/2; -Uy = (Y(1:end-1, :) + Y(2:end, :))/2; -scatter3(Ux(:), Uy(:), zeros(n*(m+1), 1), '+', 'MarkerEdgeColor', 'k') - -Vx = (X(:, 1:end-1) + X(:, 2:end))/2; -Vy = (Y(:, 1:end-1) + Y(:, 2:end))/2; -scatter3(Vx(:), Vy(:), zeros((n+1)*m, 1), '*', 'MarkerEdgeColor', 'k') - -Cx = (Vx(1:end-1, :) + Vx(2:end, :))/2; -Cy = (Uy(:, 1:end-1) + Uy(:, 2:end))/2; -scatter3(Cx(:), Cy(:), zeros(n*m, 1), '.', 'MarkerEdgeColor', 'r') - -%% Interpolators -% Some useful numbers -num_nodes = (m+1)*(n+1); -num_centers = (m+2)*(n+2); -num_u = (m+1)*n; -num_v = m*(n+1); - -% MOLE Interpolator Matrices of 'k' order -NtoC = interpolNodesToCenters2D(k, n, m); -CtoF = interpolCentersToFacesD2D(k, n, m); % V is the top part!! - -% Center to specific face, u or v -CtoV = CtoF(1:(n+1)*m, 1:num_centers); -CtoU = CtoF((n+1)*m+1:end, num_centers+1:end); - -% Node to X interpolator is pre-multiplied -NtoU = CtoU * NtoC; -NtoV = CtoV * NtoC; +Ux = CtoU * NtoC * xn; +Uy = CtoU * NtoC * yn; +scatter3(Ux, Uy, zeros(size(Ux)), '+', 'MarkerEdgeColor', 'k') -% Interpolate U values -Ugiven = sin(X); -U = NtoU * Ugiven(:); -U = reshape(U,n,m+1); +Vx = CtoV * NtoC * xn; +Vy = CtoV * NtoC * yn; +scatter3(Vx, Vy, zeros(size(Vx)), '*', 'MarkerEdgeColor', 'k') -% Interpolate V values -Vgiven = cos(Y); -V = NtoV * Vgiven(:); -V = reshape(V,n+1,m); +Cx = NtoC * xn; +Cy = NtoC * yn; +scatter3(Cx, Cy, zeros(size(Cx)), 'o', 'MarkerEdgeColor', 'r') +legend("Nodal Points", "u", "v", "Centers") +hold off -Cgiven = cos(X)-sin(Y); -C = NtoC * Cgiven(:); -C = reshape(C,n+2,m+2); +tic +D = div2DCurv(k, Cx, Cy, m, dx, n, dy, dc, nc); +toc -% Interpolate Nodal grid to centered grid. -Cx = NtoC * X(:); -Cx = reshape(Cx,n+2,m+2); +% Vector Field +U = sin(Ux); +V = cos(Vy); -Cy = NtoC * Y(:); -Cy = reshape(Cy, n+2,m+2); +% Exact Solution +C = cos(Cx) - sin(Cy); +C = reshape(C, m, n + 2)'; -scatter3(Cx(:), Cy(:), zeros((m+2)*(n+2), 1), 'o', 'MarkerEdgeColor', 'r') -legend('Nodal points', 'u', 'v', 'Centers', 'All centers') -hold off +% Approximate Solution +Ccomp = D * [U; V]; +Ccomp = reshape(Ccomp, m, n + 2)'; -tic -D = div2DCurv(k, X, Y); -toc +% Remove boundary points (divergence returns 0 on boundaries) +Cx = reshape(Cx, m, n + 2)'; +Cy = reshape(Cy, m, n + 2)'; +Cx = Cx(2:end-1, :); +Cy = Cy(2:end-1, :); +C = C(2:end-1, :); +Ccomp = Ccomp(2:end-1, :); -Ccomp = D*[reshape(U', [], 1); reshape(V', [], 1)]; -Ccomp = reshape(Ccomp, m+2, n+2); +% Join left and right sides for plotting +Cx = [Cx Cx(:, 1)]; +Cy = [Cy Cy(:, 1)]; +C = [C C(:, 1)]; +Ccomp = [Ccomp Ccomp(:, 1)]; figure subplot(2, 1, 1) -surf(Cx(2:end-1, 2:end-1), Cy(2:end-1, 2:end-1), C(2:end-1, 2:end-1), 'EdgeColor', 'none'); +surf(Cx, Cy, C, 'EdgeColor', 'none'); view([0 90]) colorbar title('Exact') @@ -106,7 +96,7 @@ axis equal shading interp subplot(2, 1, 2) -surf(Cx(2:end-1, 2:end-1), Cy(2:end-1, 2:end-1), Ccomp(2:end-1, 2:end-1)', 'EdgeColor', 'none') +surf(Cx, Cy, Ccomp, 'EdgeColor', 'none') view([0 90]) colorbar title('Approx') @@ -114,3 +104,6 @@ ylabel('y') axis equal shading interp + +l2_norm = norm(C(:) - Ccomp(:)); +disp("L2 norm: " + l2_norm) \ No newline at end of file diff --git a/examples/matlab_octave/test_div3DCurv.m b/examples/matlab_octave/test_div3DCurv.m index d09243d1b..2186f3b3c 100644 --- a/examples/matlab_octave/test_div3DCurv.m +++ b/examples/matlab_octave/test_div3DCurv.m @@ -1,15 +1,20 @@ % Tests the 3D curvilinear divergence clc close all -clear all +clear addpath('../../src/matlab_octave') % Parameters k = 2; % Order of accuracy -m = 20; % Number of nodes along x-axis -n = 30; % Number of nodes along y-axis -o = 40; % Number of nodes along z-axis +m = 20; % Number of nodes along xi-axis +n = 30; % Number of nodes along eta-axis +o = 40; % Number of nodes along kappa-axis +dx = 1 / (m - 1); +dy = 1 / (n - 1); +dz = 1 / (o - 1); +dc = [1; 1; 1; 1; 1; 1]; +nc = [0; 0; 0; 0; 0; 0]; xlength = linspace(1,m,m); ylength = linspace(1,n,n); @@ -22,85 +27,31 @@ Y = repmat(Y, [1 1 o]); [~, ~, Z] = meshgrid(1:m, 1:n, 1:o); -Ux = (X(1:end-1, : ,1:end-1)+X(2:end, : ,2:end))*0.50; -Uy = (Y(1:end-1, : ,1:end-1)+Y(2:end, : ,2:end))*0.50; -Uz = (Z(1:end-1, : ,1:end-1)+Z(2:end, : ,2:end))*0.50; -% -Vx = (X( : ,1:end-1,1:end-1)+X( : ,2:end,2:end))*0.50; -Vy = (Y( : ,1:end-1,1:end-1)+Y( : ,2:end,2:end))*0.50; -Vz = (Z( : ,1:end-1,1:end-1)+Z( : ,2:end,2:end))*0.50; -% -Wx = (X(1:end-1,1:end-1,:)+X(2:end,2:end,:))*0.50; -Wy = (Y(1:end-1,1:end-1,:)+Y(2:end,2:end,:))*0.50; -Wz = (Z(1:end-1,1:end-1,:)+Z(2:end,2:end,:))*0.50; - -% MOLE Operators want matrices to be in X,Y,Z -% for row, col, slice. Meshgrid creates a matrix -% col, row, slice. So, the permutation is needed. -X = permute( X, [2 1 3] ); -Y = permute( Y, [2 1 3] ); -Z = permute( Z, [2 1 3] ); - -Ux = permute( Ux, [2 1 3] ); -Uy = permute( Uy, [2 1 3] ); -Uz = permute( Uz, [2 1 3] ); - -Vx = permute( Vx, [2 1 3] ); -Vy = permute( Vy, [2 1 3] ); -Vz = permute( Vz, [2 1 3] ); - -Wx = permute( Wx, [2 1 3] ); -Wy = permute( Wy, [2 1 3] ); -Wz = permute( Wz, [2 1 3] ); - -%% Interpolators -% Some useful numbers -NtoC = interpolNodesToCenters3D(k+2, m-1, n-1, o-1); -CtoF3D = interpolCentersToFacesD3D(k+2,m-1,n-1,o-1);% v,u,w - -size_u = (m)*(n-1)*(o-1); -size_v = (m-1)*(n)*(o-1); -size_w = (m-1)*(n-1)*(o); -size_c = (m+1)*(n+1)*(o+1); -size_n = (m)*(n)*(o); - -CtoU = CtoF3D(1:size_u, 1:size_c); -CtoV = CtoF3D(size_u+1:size_u+size_v, size_c+1:2*size_c); -CtoW = CtoF3D(size_u+size_v+1:end, 2*size_c+1:end); - -NtoU = CtoU * NtoC; -NtoV = CtoV * NtoC; -NtoW = CtoW * NtoC; -% -% Interpolate U values -Ugiven = sin(X); %X.^2; -U = NtoU * Ugiven(:); -U = reshape(U, m, n-1 , o-1); - -% Interpolate V values -Vgiven = zeros(size(Y)); %Y.^2; -V = NtoV * Vgiven(:); -V = reshape(V, m-1, n, o-1); - -% Interpolate W values -Wgiven = zeros(size(Z)); %Z.^2; -W = NtoW * Wgiven(:); -W = reshape(W, m-1, n-1, o); - - -U = reshape(U, [], 1); -V = reshape(V, [], 1); -W = reshape(W, [], 1); - -% Get 3D curvilinear mimetic divergence -% The operator want everything in X,Y,Z, BUT!!!! -% uses the original meshgrid X,Y,Z in (y,x,z) coordinates!! -% This needs to be changed!!! This is insane !!! -Xog = permute(X, [2 1 3]); -Yog = permute(Y, [2 1 3]); -Zog = permute(Z, [2 1 3]); - -D = div3DCurv(k, Xog, Yog, Zog); +X = permute(X, [2 1 3]); +Y = permute(Y, [2 1 3]); +Z = permute(Z, [2 1 3]); + +% Interpolators +NtoC = interpolNodesToCenters3D(k, m - 1, n - 1, o - 1, dc, nc); + +CtoU = interpolCentersToFacesD1D(k, m - 1); +CtoU = kron(kron(speye(o + 1), speye(n + 1)), CtoU); + +CtoV = interpolCentersToFacesD1D(k, n - 1); +CtoV = kron(kron(speye(o + 1), CtoV), speye(m + 1)); + +CtoW = interpolCentersToFacesD1D(k, o - 1); +CtoW = kron(kron(CtoW, speye(n + 1)), speye(m + 1)); + +xc = NtoC * reshape(X, [], 1); +yc = NtoC * reshape(Y, [], 1); +zc = NtoC * reshape(Z, [], 1); + +U = CtoU * sin(xc); +V = CtoV * zeros(size(yc)); +W = CtoW * zeros(size(zc)); + +D = div3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); % Apply the operator to the field Ccomp = D*[U; V; W]; @@ -109,16 +60,13 @@ Ccomp = reshape(Ccomp, m+1, n+1, o+1); Ccomp = Ccomp(2:end-1, 2:end-1, 2:end-1); -% Compute centroids? -X = (X(1:end-1, :, :)+X(2:end, :, :))/2; -X = (X(:, 1:end-1, :)+X(:, 2:end, :))/2; -X = (X(:, :, 1:end-1)+X(:, :, 2:end))/2; -Y = (Y(1:end-1, :, :)+Y(2:end, :, :))/2; -Y = (Y(:, 1:end-1, :)+Y(:, 2:end, :))/2; -Y = (Y(:, :, 1:end-1)+Y(:, :, 2:end))/2; -Z = (Z(1:end-1, :, :)+Z(2:end, :, :))/2; -Z = (Z(:, 1:end-1, :)+Z(:, 2:end, :))/2; -Z = (Z(:, :, 1:end-1)+Z(:, :, 2:end))/2; +X = reshape(xc, m + 1, n + 1, o + 1); +Y = reshape(yc, m + 1, n + 1, o + 1); +Z = reshape(zc, m + 1, n + 1, o + 1); + +X = X(2:end-1, 2:end-1, 2:end-1); +Y = Y(2:end-1, 2:end-1, 2:end-1); +Z = Z(2:end-1, 2:end-1, 2:end-1); figure scatter3(X(:), Y(:), Z(:), 50, Ccomp(:), 'Filled'); @@ -162,3 +110,6 @@ colorbar view([0 90]) shading interp + +l2_norm = norm(div(:) - Ccomp(:)); +disp("L2 norm: " + l2_norm) \ No newline at end of file diff --git a/examples/matlab_octave/test_grad2DCurv.m b/examples/matlab_octave/test_grad2DCurv.m index ffb071dc5..454cc55be 100644 --- a/examples/matlab_octave/test_grad2DCurv.m +++ b/examples/matlab_octave/test_grad2DCurv.m @@ -1,3 +1,10 @@ +% Tests the 2D Curvilinear Gradient +% on a curvilinear mesh +% +% u = x^2 + y^2 +% Gu = <2x, 2y> +% + clc close all @@ -7,6 +14,10 @@ k = 2; % Order of accuracy m = 40; % Number of nodes along x-axis n = 40; % Number of nodes along y-axis +dx = 1 / (m - 1); +dy = 1 / (n - 1); +dc = [1; 1; 1; 1]; +nc = [0; 0; 0; 0]; [X, Y] = genCurvGrid(n, m); % [X, Y] = meshgrid(1:m, 1:n); @@ -17,53 +28,104 @@ axis equal set(gcf, 'Color', 'w') -C = X.^2+Y.^2; % Given scalar field (on a nodal grid) % Staggered logical grid [Xs, Ys] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n]); -% Interpolate scalar field from nodal logical to staggered logical -Cs = interp2(C, Xs, Ys); +Cs = Xs.^2 + Ys.^2; % Reshape the field so it can be multiplied by the operator later on C_ = reshape(Cs', [], 1); +Xs = reshape(Xs', [], 1); +Ys = reshape(Ys', [], 1); % Get 2D curvilinear mimetic gradient -G = grad2DCurv(k, X, Y); +G = grad2DCurv(k, Xs, Ys, m - 1, dx, n - 1, dy, dc, nc); % Apply the operator to the field TMP = G*C_; -Gx = TMP(1:m*(n-1)); -Gy = TMP(m*(n-1)+1:end); +Gx = TMP(1:(m+1)*(n)); +Gy = TMP((m+1)*(n)+1:end); % Reshape for visualization -Gx = reshape(Gx, m, n-1)'; -Gy = reshape(Gy, m-1, n)'; +Gx = reshape(Gx, m+1, n)'; +Gy = reshape(Gy, m, n+1)'; + +CtoU = interpolCentersToFacesD1D(k, m - 1); +CtoU = kron(speye(n + 1), CtoU); +CtoV = interpolCentersToFacesD1D(k, n - 1); +CtoV = kron(CtoV, speye(m + 1)); + +Ux = CtoU * Xs; +Uy = CtoU * Ys; +Vx = CtoV * Xs; +Vy = CtoV * Ys; + +Ux = reshape(Ux, m + 1, n)'; +Uy = reshape(Uy, m + 1, n)'; +Vx = reshape(Vx, m, n + 1)'; +Vy = reshape(Vy, m, n + 1)'; % Plot results figure set(gcf, 'Color', 'w') subplot(3, 1, 1) -surf(X, Y, C, 'EdgeColor', 'none'); +surf(reshape(Xs, m + 1, n + 1)', reshape(Ys, m + 1, n + 1)', Cs, 'EdgeColor', 'none'); +colorbar +xlabel('x') +ylabel('y') +title('C Approximate') +axis equal +view([0 90]) +shading interp +subplot(3, 1, 2) +surf(Ux, Uy, Gx, 'EdgeColor', 'none'); +colorbar +xlabel('x') +ylabel('y') +title('U Approximate') +axis equal +view([0 90]) +shading interp +subplot(3, 1, 3) +surf(Vx, Vy, Gy, 'EdgeColor', 'none'); +colorbar +xlabel('x') +ylabel('y') +title('V Approximate') +axis equal +view([0 90]) +shading interp + +figure +set(gcf, 'Color', 'w') +subplot(3, 1, 1) +surf(reshape(Xs, m + 1, n + 1)', reshape(Ys, m + 1, n + 1)', Cs, 'EdgeColor', 'none'); colorbar xlabel('x') ylabel('y') -title('C') +title('C Analytical') axis equal view([0 90]) shading interp subplot(3, 1, 2) -surf((X(1:end-1, :)+X(2:end, :))/2, (Y(1:end-1, :)+Y(2:end, :))/2, Gx, 'EdgeColor', 'none'); +surf(Ux, Uy, 2 * Ux, 'EdgeColor', 'none'); colorbar xlabel('x') ylabel('y') -title('U') +title('U Analytical') axis equal view([0 90]) shading interp subplot(3, 1, 3) -surf((X(:, 1:end-1)+X(:, 2:end))/2, (Y(:, 1:end-1)+Y(:, 2:end))/2, Gy, 'EdgeColor', 'none'); +surf(Vx, Vy, 2 * Vy, 'EdgeColor', 'none'); colorbar xlabel('x') ylabel('y') -title('V') +title('V Analytical') axis equal view([0 90]) shading interp + +l2_norm_u = norm(Gx - 2 * Ux); +l2_norm_v = norm(Gy - 2 * Vy); + +disp("L2 norm of U: " + l2_norm_u) +disp("L2 norm of V: " + l2_norm_v) \ No newline at end of file diff --git a/examples/matlab_octave/test_grad3DCurv.m b/examples/matlab_octave/test_grad3DCurv.m index cf38d9495..aab7afd2b 100644 --- a/examples/matlab_octave/test_grad3DCurv.m +++ b/examples/matlab_octave/test_grad3DCurv.m @@ -1,58 +1,102 @@ clc close all +clear addpath('../../src/matlab_octave') % Parameters -k = 2; % Order of accuracy -m = 20; % Number of nodes along x-axis -n = 20; % Number of nodes along y-axis -o = 20; % Number of nodes along z-axis +k = 2; % Order of accuracy +m = 20; % Number of nodes along x-axis +n = 20; % Number of nodes along y-axis +o = 20; % Number of nodes along z-axis +dx = 1 / (m - 1); % Step size along x-axis +dy = 1 / (n - 1); % Step size along y-axis +dz = 1 / (o - 1); % Step size along z-axis +dc = [1; 1; 1; 1; 1; 1]; % Dirichlet Coefficients +nc = [0; 0; 0; 0; 0; 0]; % Neumann Coefficients [X, Y] = genCurvGrid(n, m); X = repmat(X, [1 1 o]); Y = repmat(Y, [1 1 o]); [~, ~, Z] = meshgrid(1:m, 1:n, 1:o); -C = X.^2+Y.^2+Z.^2; % Given scalar field (on a nodal grid) +% Reshape to vectors +xn = reshape(permute(X, [2, 1, 3]), [], 1); +yn = reshape(permute(Y, [2, 1, 3]), [], 1); +zn = reshape(permute(Z, [2, 1, 3]), [], 1); -% Plot the physical grid -scatter3(X(:), Y(:), Z(:), 50, C(:), 'Filled'); +% Interpolators +INC = interpolNodesToCenters3D(k, m - 1, n - 1, o - 1); +ICFx = interpolCentersToFacesD1D(k, m - 1); +ICFy = interpolCentersToFacesD1D(k, n - 1); +ICFz = interpolCentersToFacesD1D(k, o - 1); +Im = speye(m + 1); +In = speye(n + 1); +Io = speye(o + 1); +ICFx = kron(kron(Io, In), ICFx); +ICFy = kron(kron(Io, ICFy), Im); +ICFz = kron(kron(ICFz, In), Im); + +% Staggered Grid +xc = INC * xn; +yc = INC * yn; +zc = INC * zn; +C = xc.^2 + yc.^2 + zc.^2; + +Ux = ICFx * xc; +Uy = ICFy * yc; +Uz = ICFz * zc; + +Vx = ICFx * xc; +Vy = ICFy * yc; +Vz = ICFz * zc; + +Wx = ICFx * xc; +Wy = ICFy * yc; +Wz = ICFz * zc; + +% Plot the centers +scatter3(xc(:), yc(:), zc(:), 50, C(:), 'Filled'); title('Given scalar field') xlabel('x') ylabel('y') zlabel('z') axis equal -% Staggered logical grid -[Xs, Ys, Zs] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n], [1 1.5 : 1 : o-0.5 o]); -% Interpolate scalar field from nodal logical to staggered logical -Cs = interp3(C, Xs, Ys, Zs); -% Reshape the field so it can be multiplied by the operator later on -C_ = reshape(permute(Cs, [2, 1, 3]), [], 1); - % Get 3D curvilinear mimetic gradient -G = grad3DCurv(k, X, Y, Z); +G = grad3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); % Apply the operator to the field -TMP = G*C_; -Gx = TMP(1:m*(n-1)*(o-1)); -Gy = TMP(m*(n-1)*(o-1)+1:m*(n-1)*(o-1)+(m-1)*n*(o-1)); -Gz = TMP(m*(n-1)*(o-1)+(m-1)*n*(o-1)+1:end); +TMP = G*C; +Gx = TMP(1 : m * (n + 1) * (o + 1)); +Gy = TMP(m * (n + 1) * (o + 1) + 1 : m * (n + 1) * (o + 1) + (m + 1) * n * (o + 1)); +Gz = TMP(m * (n + 1) * (o + 1) + (m + 1) * n * (o + 1) + 1 : end); % Reshape for visualization -Gx = Gx(1:m*(n-1)); -Gx = reshape(Gx, m, n-1)'; -Gy = Gy(1:(m-1)*n); -Gy = reshape(Gy, m-1, n)'; -Gz = reshape(Gz, m-1, n-1, o); -Gz = squeeze(Gz(:, 1, :))'; +Uxp = Ux(1 : m * (n + 1)); +Uyp = Uy(1 : m * (n + 1)); +Gxp = Gx(1 : m * (n + 1)); +Uxp = reshape(Uxp, m, n + 1)'; +Uyp = reshape(Uyp, m, n + 1)'; +Gxp = reshape(Gxp, m, n + 1)'; + +Vxp = Vx(1 : (m + 1) * n); +Vyp = Vy(1 : (m + 1) * n); +Gyp = Gy(1 : (m + 1) * n); +Vxp = reshape(Vxp, m + 1, n)'; +Vyp = reshape(Vyp, m + 1, n)'; +Gyp = reshape(Gyp, m + 1, n)'; + +Wxp = reshape(Wx, m + 1, n + 1, o); +Wzp = reshape(Wz, m + 1, n + 1, o); +Gzp = reshape(Gz, m + 1, n + 1, o); +Wxp = squeeze(Wxp(:, 1, :)); +Wzp = squeeze(Wzp(:, 1, :)); +Gzp = squeeze(Gzp(:, 1, :)); figure -Xu = (X(1:end-1, :, 1)+X(2:end, :, 1))/2; -Yu = (Y(1:end-1, :, 1)+Y(2:end, :, 1))/2; -surf(Xu, Yu, Gx, 'EdgeColor', 'none') -title('Gx') +surf(Uxp, Uyp, Gxp, 'EdgeColor', 'none') +title('Gx Approximate') xlabel('x') ylabel('y') zlabel('z') @@ -61,10 +105,8 @@ shading interp figure -Xv = (X(:, 1:end-1, 1)+X(:, 2:end, 1))/2; -Yv = (Y(:, 1:end-1, 1)+Y(:, 2:end, 1))/2; -surf(Xv, Yv, Gy, 'EdgeColor', 'none') -title('Gy') +surf(Vxp, Vyp, Gyp, 'EdgeColor', 'none') +title('Gy Approximate') xlabel('x') ylabel('y') zlabel('z') @@ -73,13 +115,49 @@ shading interp figure -Xw = squeeze((X(1, 1:end-1, :)+X(1, 2:end, :))/2); -Zw = squeeze((Z(1, 1:end-1, :)+Z(1, 2:end, :))/2); -surf(Xw, Zw, Gz', 'EdgeColor', 'none') -title('Gz') +surf(Wxp, Wzp, Gzp, 'EdgeColor', 'none') +title('Gz Approximate') xlabel('x') ylabel('z') zlabel('y') axis equal view([0 90]) shading interp + +figure +surf(Uxp, Uyp, 2 * Uxp, 'EdgeColor', 'none') +title('Gx Analytical') +xlabel('x') +ylabel('y') +zlabel('Gx') +axis equal +view([0 90]) +shading interp + +figure +surf(Vxp, Vyp, 2 * Vyp, 'EdgeColor', 'none') +title('Gy Analytical') +xlabel('x') +ylabel('y') +zlabel('Gy') +axis equal +view([0 90]) +shading interp + +figure +surf(Wxp, Wzp, 2 * Wzp, 'EdgeColor', 'none') +title('Gz Analytical') +xlabel('x') +ylabel('z') +zlabel('Gz') +axis equal +view([0 90]) +shading interp + +l2_norm_u = norm(Gx - 2 * Ux); +l2_norm_v = norm(Gy - 2 * Vy); +l2_norm_w = norm(Gz - 2 * Wz); + +disp("L2 norm of U: " + l2_norm_u) +disp("L2 norm of V: " + l2_norm_v) +disp("L2 norm of W: " + l2_norm_w) \ No newline at end of file diff --git a/examples/matlab_octave/test_lap3DCurv.m b/examples/matlab_octave/test_lap3DCurv.m index 4de8fa0f9..e00775b08 100644 --- a/examples/matlab_octave/test_lap3DCurv.m +++ b/examples/matlab_octave/test_lap3DCurv.m @@ -1,5 +1,5 @@ % Tests the 3D curvilinear laplacian - +clear clc close all @@ -10,25 +10,31 @@ m = 25; % Number of nodes along x-axis n = 20; % Number of nodes along y-axis o = 15; % Number of nodes along z-axis +dx = 1 / (m - 1); +dy = 1 / (n - 1); +dz = 1 / (o - 1); +dc = [1; 1; 1; 1; 1; 1]; +nc = [0; 0; 0; 0; 0; 0]; [X, Y] = genCurvGrid(n, m); X = repmat(X, [1 1 o]); Y = repmat(Y, [1 1 o]); [~, ~, Z] = meshgrid(1:m, 1:n, 1:o); -F = X.^2+Y.^2+Z.^2; % Unknown function +% Interpolate from Nodes to Centers +INC = interpolNodesToCenters3D(k, m - 1, n - 1, o - 1, dc, nc); -% Staggered logical grid -[Xs, Ys, Zs] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n], [1 1.5 : 1 : o-0.5 o]); -% Interpolate scalar field from nodal logical to staggered logical -RHS = interp3(F, Xs, Ys, Zs); +xc = INC * reshape(permute(X, [2, 1, 3]), [], 1); +yc = INC * reshape(permute(Y, [2, 1, 3]), [], 1); +zc = INC * reshape(permute(Z, [2, 1, 3]), [], 1); -RHS = reshape(permute(RHS, [2, 1, 3]), [], 1); +ue = xc.^2 + yc.^2 + zc.^2; +RHS = ue; % Get 3D curvilinear mimetic divergence -D = div3DCurv(k, X, Y, Z); +D = div3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); % Get 3D curvilinear mimetic gradient -G = grad3DCurv(k, X, Y, Z); +G = grad3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); % Dirichlet BCs BC = robinBC3D(k, m-1, 1, n-1, 1, o-1, 1, 1, 0); % Laplacian operator with BCs @@ -38,13 +44,12 @@ RHS(idx) = 6; % RHS = f''(x, y, z) in the inner domain % Solve the system of linear equations -SOL = L\RHS; - -SOL = permute(reshape(SOL, m+1, n+1, o+1), [2, 1, 3]); +SOL = L \ RHS; % Plot the exact solution +figure subplot(2, 1, 1) -scatter3(X(:), Y(:), Z(:), 100, F(:), 'Filled'); +scatter3(xc, yc, zc, 100, ue, 'Filled'); title('Exact') xlabel('x') ylabel('y') @@ -54,11 +59,13 @@ % Plot the approximation subplot(2, 1, 2) -SOL = SOL(2:end, 2:end, 2:end); % So we can plot with X, Y and Z -scatter3(X(:), Y(:), Z(:), 100, SOL(:), 'Filled'); +scatter3(xc, yc, zc, 100, SOL, 'Filled'); title('Approximation') xlabel('x') ylabel('y') zlabel('z') axis equal colorbar + +l2_norm = norm(ue - SOL); +disp("L2 norm of error: " + l2_norm) \ No newline at end of file diff --git a/examples/matlab_octave/test_lap3DCurv_case2.m b/examples/matlab_octave/test_lap3DCurv_case2.m index 5b8e3c948..8c62bc3bf 100644 --- a/examples/matlab_octave/test_lap3DCurv_case2.m +++ b/examples/matlab_octave/test_lap3DCurv_case2.m @@ -1,15 +1,20 @@ % Tests the 3D curvilinear laplacian - +clear clc close all addpath('../../src/matlab_octave') % Parameters -k = 2; % Order of accuracy -m = 20; % Number of nodes along x-axis -n = 20; % Number of nodes along y-axis -o = 10; % Number of nodes along z-axis +k = 2; % Order of accuracy +m = 20; % Number of nodes along x-axis +n = 20; % Number of nodes along y-axis +o = 10; % Number of nodes along z-axis +dx = 1 / (m - 1); % Step size along x-axis +dy = 1 / (n - 1); % Step size along y-axis +dz = 1 / (o - 1); % Step size along z-axis +dc = [1; 1; 1; 1; 1; 1]; % Dirichlet Coefficients +nc = [0; 0; 0; 0; 0; 0]; % Neumann Coefficients a = -pi; % xlim(0) b = pi; % xlim(1) @@ -20,19 +25,20 @@ [X, Y, Z] = meshgrid(linspace(a, b, m), linspace(c, d, n), linspace(e, f, o)); X = X + sin(Y); -F = -(X.^2+Y.^2+Z.^2); % Unknown function +% Interpolate Nodes to Centers +INC = interpolNodesToCenters3D(k, m - 1, n - 1, o - 1, dc, nc); -% Staggered logical grid -[Xs, Ys, Zs] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n], [1 1.5 : 1 : o-0.5 o]); -% Interpolate scalar field from nodal logical to staggered logical -RHS = interp3(F, Xs, Ys, Zs); +xc = INC * reshape(permute(X, [2, 1, 3]), [], 1); +yc = INC * reshape(permute(Y, [2, 1, 3]), [], 1); +zc = INC * reshape(permute(Z, [2, 1, 3]), [], 1); -RHS = reshape(permute(RHS, [2, 1, 3]), [], 1); +ue = -(xc.^2 + yc.^2 + zc.^2); % Unknown function +RHS = ue; % Get 3D curvilinear mimetic divergence -D = div3DCurv(k, X, Y, Z); +D = div3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); % Get 3D curvilinear mimetic gradient -G = grad3DCurv(k, X, Y, Z); +G = grad3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); % Dirichlet BCs BC = robinBC3D(k, m-1, 1, n-1, 1, o-1, 1, 1, 0); % Laplacian operator with BCs @@ -42,13 +48,11 @@ RHS(idx) = -6; % RHS = f''(x, y, z) in the inner domain % Solve the system of linear equations -SOL = L\RHS; - -SOL = permute(reshape(SOL, m+1, n+1, o+1), [2, 1, 3]); +SOL = L \ RHS; % Plot the exact solution subplot(2, 1, 1) -scatter3(X(:), Y(:), Z(:), 100, F(:), 'Filled'); +scatter3(xc, yc, zc, 100, ue, 'Filled'); title('Exact') xlabel('x') ylabel('y') @@ -58,11 +62,13 @@ % Plot the approximation subplot(2, 1, 2) -SOL = SOL(2:end, 2:end, 2:end); % So we can plot with X, Y and Z -scatter3(X(:), Y(:), Z(:), 100, SOL(:), 'Filled'); +scatter3(xc, yc, zc, 100, SOL, 'Filled'); title('Approximation') xlabel('x') ylabel('y') zlabel('z') axis equal colorbar + +l2_norm = norm(ue - SOL); +disp("L2 norm of error: " + l2_norm) \ No newline at end of file diff --git a/src/matlab_octave/jacobian2D.m b/src/matlab_octave/jacobian2D.m index 5e064335a..b2948f7ea 100644 --- a/src/matlab_octave/jacobian2D.m +++ b/src/matlab_octave/jacobian2D.m @@ -76,8 +76,6 @@ Gn = kron(IFCy,Im) * kron(Gn,Im); G = [Ge; Gn]; - X = reshape(X',[],1); - Y = reshape(Y',[],1); metrics = G * [X Y]; Xe = metrics(1:numC,1); diff --git a/src/matlab_octave/jacobian3D.m b/src/matlab_octave/jacobian3D.m index 4a85da851..623629dec 100644 --- a/src/matlab_octave/jacobian3D.m +++ b/src/matlab_octave/jacobian3D.m @@ -105,9 +105,6 @@ Gk = kron(kron(IFCz, In), Im) * kron(kron(Gk, In), Im); G = [Ge; Gn; Gk]; - X = reshape(X', [], 1); - Y = reshape(Y', [], 1); - Z = reshape(Z', [], 1); metrics = G * [X Y Z]; Xe = metrics(1:numC, 1); From f1b5f0e1a34eb05225b60ed24d7a03a01382da91 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Tue, 5 May 2026 13:40:01 -0700 Subject: [PATCH 09/10] doc strings --- examples/matlab_octave/curvConvergenceTest.m | 2 +- src/matlab_octave/div2DCurv.m | 13 +++--- src/matlab_octave/div2DCurvLegacy.m | 13 +++++- src/matlab_octave/div3DCurv.m | 27 ++++++------- src/matlab_octave/div3DCurvLegacy.m | 15 ++++++- src/matlab_octave/grad2DCurv.m | 13 +++--- src/matlab_octave/grad2DCurvLegacy.m | 15 ++++++- src/matlab_octave/grad3DCurv.m | 23 +++++------ src/matlab_octave/grad3DCurvLegacy.m | 16 +++++++- src/matlab_octave/jacobian2D.m | 22 ++++------ src/matlab_octave/jacobian2DLegacy.m | 16 ++++---- src/matlab_octave/jacobian3D.m | 42 +++++--------------- src/matlab_octave/jacobian3DLegacy.m | 26 +++++------- 13 files changed, 126 insertions(+), 117 deletions(-) diff --git a/examples/matlab_octave/curvConvergenceTest.m b/examples/matlab_octave/curvConvergenceTest.m index 564bb0d2c..385304947 100644 --- a/examples/matlab_octave/curvConvergenceTest.m +++ b/examples/matlab_octave/curvConvergenceTest.m @@ -1,4 +1,4 @@ -%% Convegence Test for Cuvilinear Operators +%% Convergence Test for Cuvilinear Operators % % In 2D space, MOLE operators can only calculate the xi and eta derivatives % On a curvilinear domain, the x and xi, and y and eta directions are diff --git a/src/matlab_octave/div2DCurv.m b/src/matlab_octave/div2DCurv.m index e3aaab8ab..3a907d1f6 100644 --- a/src/matlab_octave/div2DCurv.m +++ b/src/matlab_octave/div2DCurv.m @@ -1,13 +1,10 @@ function D = div2DCurv(k, X, Y, m, dx, n, dy, dc, nc) % PURPOSE -% Returns a 2D curvilinear mimetic divergence operator +% Returns a 2D curvilinear mimetic divergence operator. If optional +% arguments are specified, it acts on the extended faces +% (normal faces plus boundaries) % % DESCRIPTION -% Returns: -% D : 2D curvilinear mimetic divergence operator. If optional -% arguments are specified, it acts on the extended faces -% (normal faces plus boundaries) -% % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid centers if @@ -22,11 +19,11 @@ % boundaries, resp.) % (optional) nc : b0 (4x1 vector for left, right, bottom, top % boundaries, resp.) -% +% % SYNTAX % D = div2DCurv(k, X, Y) % D = div2DCurv(k, X, Y, m, dx, n, dy, dc, nc) -% +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/div2DCurvLegacy.m b/src/matlab_octave/div2DCurvLegacy.m index c90cab807..f20914847 100644 --- a/src/matlab_octave/div2DCurvLegacy.m +++ b/src/matlab_octave/div2DCurvLegacy.m @@ -1,9 +1,20 @@ function D = div2DCurvLegacy(k, X, Y) % % ---------------------------------------------------------------------------- -% !!! WARNING: DEPRICATED BY div2DCurv.m !!! +% !!! WARNING: DEPRECATED BY div2DCurv.m !!! % ---------------------------------------------------------------------------- % +% PURPOSE +% Returns a 2D curvilinear mimetic divergence operator +% +% DESCRIPTION +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid nodes +% Y : y-coordinates (physical) of meshgrid nodes +% +% SYNTAX +% D = div2DCurvLegacy(k, X, Y) % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/div3DCurv.m b/src/matlab_octave/div3DCurv.m index 58727f01b..69d5496d1 100644 --- a/src/matlab_octave/div3DCurv.m +++ b/src/matlab_octave/div3DCurv.m @@ -1,34 +1,33 @@ function D = div3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % PURPOSE -% Returns a 3D curvilinear mimetic divergence operator +% Returns a 3D curvilinear mimetic divergence operator. If optional +% arguments are specified, it acts on the extended faces +% (normal faces plus boundaries) % % DESCRIPTION -% Returns: -% D : 3D curvilinear mimetic divergence operator. If optional -% arguments are specified, it acts on the extended faces -% (normal faces plus boundaries) -% % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes % Y : y-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes +% Z : z-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes % (optional) m : Number of cells in xi direction % (optional) dx : Step size in xi direction % (optional) n : Number of cells in eta direction % (optional) dy : Step size in eta direction % (optional) o : Number of cells in kappa direction % (optional) dz : Step size in kappa direction -% (optional) dc : a0 (6x1 vector for left, right, bottom, top -% boundaries, resp.) -% (optional) nc : b0 (6x1 vector for left, right, bottom, top -% boundaries, resp.) -% +% (optional) dc : a0 (6x1 vector for left, right, bottom, top, front, and +% back boundaries, resp.) +% (optional) nc : b0 (6x1 vector for left, right, bottom, top, front, and +% back boundaries, resp.) +% % SYNTAX % D = div3DCurv(k, X, Y, Z) % D = div3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) -% +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). @@ -36,8 +35,8 @@ % ---------------------------------------------------------------------------- if nargin ~= 4 && nargin ~= 12 - error("div2DCurv:InvalidNumArgs", ... - "div2DCurv expects 4 or 12 arguments") + error("div3DCurv:InvalidNumArgs", ... + "div3DCurv expects 4 or 12 arguments") elseif nargin == 4 D = div3DCurvLegacy(k, X, Y, Z); return; diff --git a/src/matlab_octave/div3DCurvLegacy.m b/src/matlab_octave/div3DCurvLegacy.m index 6acc77be9..ce71c2f3a 100644 --- a/src/matlab_octave/div3DCurvLegacy.m +++ b/src/matlab_octave/div3DCurvLegacy.m @@ -1,9 +1,22 @@ function D = div3DCurvLegacy(k, X, Y, Z) % % ---------------------------------------------------------------------------- -% !!! WARNING: DEPRICATED BY div2DCurv.m !!! +% !!! WARNING: DEPRECATED BY div3DCurv.m !!! % ---------------------------------------------------------------------------- % +% PURPOSE +% Returns a 3D curvilinear mimetic divergence operator. +% +% DESCRIPTION +% % Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid nodes +% Y : y-coordinates (physical) of meshgrid nodes +% Z : z-coordinates (physical) of meshgrid nodes +% +% SYNTAX +% D = div3DCurvLegacy(k, X, Y, Z) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/grad2DCurv.m b/src/matlab_octave/grad2DCurv.m index d329200b9..49335821a 100644 --- a/src/matlab_octave/grad2DCurv.m +++ b/src/matlab_octave/grad2DCurv.m @@ -1,13 +1,10 @@ function G = grad2DCurv(k, X, Y, m, dx, n, dy, dc, nc) % PURPOSE -% Returns a 2D curvilinear mimetic gradient operator +% Returns a 2D curvilinear mimetic gradient operator. If optional +% arguments are specified, it outputs to the extended +% faces (normal faces plus boundaries) % % DESCRIPTION -% Returns: -% G : 2D curvilinear mimetic gradient operator. If optional -% arguments are specified, it outputs to the extended -% faces (normal faces plus boundaries) -% % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid centers if @@ -22,11 +19,11 @@ % boundaries, resp.) % (optional) nc : b0 (4x1 vector for left, right, bottom, top % boundaries, resp.) -% +% % SYNTAX % G = grad2DCurv(k, X, Y) % G = grad2DCurv(k, X, Y, m, dx, n, dy, dc, nc) -% +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/grad2DCurvLegacy.m b/src/matlab_octave/grad2DCurvLegacy.m index c30cfd744..804fbcaf3 100644 --- a/src/matlab_octave/grad2DCurvLegacy.m +++ b/src/matlab_octave/grad2DCurvLegacy.m @@ -1,10 +1,21 @@ function G = grad2DCurvLegacy(k, X, Y) % % ---------------------------------------------------------------------------- -% !!! WARNING: DEPRICATED BY grad2DCurv.m !!! +% !!! WARNING: DEPRECATED BY grad2DCurv.m !!! % ---------------------------------------------------------------------------- % -% Returns a 2D curvilinear mimetic gradient +% PURPOSE +% Returns a 2D curvilinear mimetic gradient operator +% +% DESCRIPTION +% % Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid nodes +% Y : y-coordinates (physical) of meshgrid nodes +% +% SYNTAX +% G = grad2DCurvLegacy(k, X, Y) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/grad3DCurv.m b/src/matlab_octave/grad3DCurv.m index c463d467a..a3c0322b3 100644 --- a/src/matlab_octave/grad3DCurv.m +++ b/src/matlab_octave/grad3DCurv.m @@ -1,34 +1,33 @@ function G = grad3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % PURPOSE -% Returns a 3D curvilinear mimetic gradient operator +% Returns a 3D curvilinear mimetic gradient operator. If optional +% arguments are specified, it outputs to the extended +% faces (normal faces plus boundaries) % % DESCRIPTION -% Returns: -% G : 3D curvilinear mimetic gradient operator. If optional -% arguments are specified, it outputs to the extended -% faces (normal faces plus boundaries) -% % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes % Y : y-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes +% Z : Z-coordinates (physical) of meshgrid centers if +% optional arguments are specified, else nodes % (optional) m : Number of cells in xi direction % (optional) dx : Step size in xi direction % (optional) n : Number of cells in eta direction % (optional) dy : Step size in eta direction % (optional) o : Number of cells in kappa direction % (optional) dz : Step size in kappa direction -% (optional) dc : a0 (6x1 vector for left, right, bottom, top -% boundaries, resp.) -% (optional) nc : b0 (6x1 vector for left, right, bottom, top -% boundaries, resp.) -% +% (optional) dc : a0 (6x1 vector for left, right, bottom, top, front, and +% back boundaries, resp.) +% (optional) nc : b0 (6x1 vector for left, right, bottom, top, front, and +% back boundaries, resp.) +% % SYNTAX % G = grad3DCurv(k, X, Y, Z) % G = grad3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) -% +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/grad3DCurvLegacy.m b/src/matlab_octave/grad3DCurvLegacy.m index dad96f4c1..af1f0d853 100644 --- a/src/matlab_octave/grad3DCurvLegacy.m +++ b/src/matlab_octave/grad3DCurvLegacy.m @@ -1,10 +1,22 @@ function G = grad3DCurvLegacy(k, X, Y, Z) % % ---------------------------------------------------------------------------- -% !!! WARNING: DEPRICATED BY grad2DCurv.m !!! +% !!! WARNING: DEPRECATED BY grad3DCurv.m !!! % ---------------------------------------------------------------------------- % -% Returns a 3D curvilinear mimetic gradient +% PURPOSE +% Returns a 3D curvilinear mimetic gradient operator +% +% DESCRIPTION +% Parameters: +% k : Order of accuracy +% X : x-coordinates (physical) of meshgrid nodes +% Y : y-coordinates (physical) of meshgrid nodes +% Z : Z-coordinates (physical) of meshgrid nodes +% +% SYNTAX +% G = grad3DCurvLegacy(k, X, Y, Z) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/jacobian2D.m b/src/matlab_octave/jacobian2D.m index b2948f7ea..61ad85c05 100644 --- a/src/matlab_octave/jacobian2D.m +++ b/src/matlab_octave/jacobian2D.m @@ -1,20 +1,9 @@ function [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, m, dx, n, dy, dc, nc) % PURPOSE -% 2D Jacobian Metrics for Curvilinear Operators -% -% DESCRIPTION -% Returns: -% J : Determinant of the Jacobian (XeYn - XnYe) on the -% centers if optional arguments are specified, else nodes -% Xe : dx/de metric on the centers if optional arguments are -% specified, else nodes -% Xn : dx/dn metric on the centers if optional arguments are -% specified, else nodes -% Ye : dy/de metric on the centers if optional arguments are -% specified, else nodes -% Yn : dy/dn metric on the centers if optional arguments are -% specified, else nodes +% Returns the Jacobian metrics (Xe, Xn, Ye, Yn, and J = XeYn - XnYe) of a +% mesh % +% DESCRIPTION % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid centers if @@ -29,6 +18,11 @@ % boundaries, resp.) % (optional) nc : b0 (4x1 vector for left, right, bottom, top % boundaries, resp.) +% +% SYNTAX +% [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y) +% [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, m, dx, n, dy, dc, nc) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/jacobian2DLegacy.m b/src/matlab_octave/jacobian2DLegacy.m index 25bcc7b39..ae56b622c 100644 --- a/src/matlab_octave/jacobian2DLegacy.m +++ b/src/matlab_octave/jacobian2DLegacy.m @@ -1,20 +1,22 @@ function [J, Xe, Xn, Ye, Yn] = jacobian2DLegacy(k, X, Y) % % ---------------------------------------------------------------------------- -% !!! WARNING: DEPRICATED BY jacobian2D.m !!! +% !!! WARNING: DEPRECATED BY jacobian2D.m !!! % ---------------------------------------------------------------------------- % -% Returns: -% J : Determinant of the Jacobian (XeYn - XnYe) on the nodes -% Xe : dx/de metric on the nodes -% Xn : dx/dn metric on the nodes -% Ye : dy/de metric on the nodes -% Yn : dy/dn metric on the nodes +% PURPOSE +% Returns the Jacobian metrics (Xe, Xn, Ye, Yn, and J = XeYn - XnYe) of a +% mesh % +% DESCRIPTION % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid nodes % Y : y-coordinates (physical) of meshgrid nodes +% +% SYNTAX +% [J, Xe, Xn, Ye, Yn] = jacobian2DLegacy(k, X, Y) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). diff --git a/src/matlab_octave/jacobian3D.m b/src/matlab_octave/jacobian3D.m index 623629dec..0f8544b21 100644 --- a/src/matlab_octave/jacobian3D.m +++ b/src/matlab_octave/jacobian3D.m @@ -1,30 +1,8 @@ function [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) % PURPOSE -% 3D Jacobian metrics for curvilinear operators +% Returns the 3D jacobian metrics of a mesh % % DESCRIPTION -% Returns: -% J : Determinant of the Jacobian on the centers if optional -% arguments are specified, else nodes -% Xe : dx/de metric on the centers if optional arguments are -% specified, else nodes -% Xn : dx/dn metric on the centers if optional arguments are -% specified, else nodes -% Xk : dx/dk metric on the centers if optional arguments are -% specified, else nodes -% Ye : dy/de metric on the centers if optional arguments are -% specified, else nodes -% Yn : dy/dn metric on the centers if optional arguments are -% specified, else nodes -% Yk : dy/dk metric on the centers if optional arguments are -% specified, else nodes -% Ze : dz/de metric on the centers if optional arguments are -% specified, else nodes -% Zn : dz/dn metric on the centers if optional arguments are -% specified, else nodes -% Zk : dz/dk metric on the centers if optional arguments are -% specified, else nodes -% % Parameters: % k : Order of accuracy % X : x-coordinates (physical) of meshgrid centers if optional @@ -39,15 +17,15 @@ % (optional) dy : Step size in eta direction % (optional) o : Number of cells in kappa direction % (optional) dz : Step size in kappa direction -% (optional) dc : a0 (6x1 vector for left, right, bottom, top -% boundaries, resp.) -% (optional) nc : b0 (6x1 vector for left, right, bottom, top -% boundaries, resp.) -% +% (optional) dc : a0 (6x1 vector for left, right, bottom, top, front, and +% back boundaries, resp.) +% (optional) nc : b0 (6x1 vector for left, right, bottom, top, front, and +% back boundaries, resp.) +% % SYNTAX -% [J, Xe, Xn, Xv, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z) -% [J, Xe, Xn, Xv, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) -% +% [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z) +% [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). @@ -56,7 +34,7 @@ if nargin ~= 4 && nargin ~= 12 error("jacobian3D:InvalidNumArgs", ... - "jacobian3D expects 4 or 9 arguments") + "jacobian3D expects 4 or 12 arguments") elseif nargin == 4 [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3DLegacy(k, X, Y, Z); return; diff --git a/src/matlab_octave/jacobian3DLegacy.m b/src/matlab_octave/jacobian3DLegacy.m index 562b8392f..90f7adabe 100644 --- a/src/matlab_octave/jacobian3DLegacy.m +++ b/src/matlab_octave/jacobian3DLegacy.m @@ -1,26 +1,22 @@ function [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3DLegacy(k, X, Y, Z) % % ---------------------------------------------------------------------------- -% !!! WARNING: DEPRICATED BY jacobian3D.m !!! +% !!! WARNING: DEPRECATED BY jacobian3D.m !!! % ---------------------------------------------------------------------------- % -% Returns: -% J : Determinant of the Jacobian -% Xe : dx/de metric -% Xn : dx/dn metric -% Xc : dx/dc metric -% Ye : dy/de metric -% Yn : dy/dn metric -% Yc : dy/dc metric -% Ze : dz/de metric -% Zn : dz/dn metric -% Zc : dz/dc metric +% PUPOSE +% Returns the 3D jacobian metrics of a mesh % +% DESCRIPTION % Parameters: % k : Order of accuracy -% X : x-coordinates (physical) of meshgrid -% Y : y-coordinates (physical) of meshgrid -% Z : z-coordinates (physical) of meshgrid +% X : x-coordinates (physical) of meshgrid nodes +% Y : y-coordinates (physical) of meshgrid nodes +% Z : z-coordinates (physical) of meshgrid nodes +% +% SYNTAX +% [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3DLegacy(k, X, Y, Z) +% % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later % © 2008-2024 San Diego State University Research Foundation (SDSURF). From 4bf5801a42a681212181773b77ca381718f1c2d7 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Wed, 3 Jun 2026 12:01:57 -0700 Subject: [PATCH 10/10] tweaks --- examples/matlab_octave/curvConvergenceTest.m | 20 ++++------ .../matlab_octave/elliptic2DCurvPeriodic.m | 15 +++----- examples/matlab_octave/test_div2DCurv.m | 11 +++--- examples/matlab_octave/test_div3DCurv.m | 24 ++++++++---- examples/matlab_octave/test_grad2DCurv.m | 29 +++++++------- examples/matlab_octave/test_grad3DCurv.m | 8 +++- examples/matlab_octave/test_lap3DCurv.m | 8 +++- examples/matlab_octave/test_lap3DCurv_case2.m | 8 +++- src/matlab_octave/div2DCurv.m | 27 ++++++++----- src/matlab_octave/div3DCurv.m | 33 ++++++++++------ src/matlab_octave/grad2DCurv.m | 29 +++++++++----- src/matlab_octave/grad3DCurv.m | 33 ++++++++++------ src/matlab_octave/jacobian2D.m | 31 ++++++++++----- src/matlab_octave/jacobian3D.m | 38 +++++++++++++------ 14 files changed, 197 insertions(+), 117 deletions(-) diff --git a/examples/matlab_octave/curvConvergenceTest.m b/examples/matlab_octave/curvConvergenceTest.m index 385304947..8e759b5b1 100644 --- a/examples/matlab_octave/curvConvergenceTest.m +++ b/examples/matlab_octave/curvConvergenceTest.m @@ -85,21 +85,19 @@ rsC = [1 1+dy/2:dy:2-dy/2 2]; tsC = [pi pi-dx/2:-dx:dx/2 0]; [TSC,RSC] = meshgrid(tsC,rsC); - xc = RSC .* cos(TSC); xc = reshape(xc',[],1); - yc = RSC .* sin(TSC); yc = reshape(yc',[],1); + X = RSC .* cos(TSC); + Y = RSC .* sin(TSC); % Build operators curG = grad2DCurv(k(i),xn,yn); curD = div2DCurv(k(i),xn,yn); curL = curD * curG; - newG = grad2DCurv(k(i),xc,yc,m,dx,n,dy,dc,nc); - newD = div2DCurv(k(i),xc,yc,m,dx,n,dy,dc,nc); + newG = grad2DCurv(k(i),X,Y,dc,nc); + newD = div2DCurv(k(i),X,Y,dc,nc); newL = newD * newG; % Boundary Conditions - X = reshape(xc,m+2,n+2)'; % Reshape for plotting and easier BC - Y = reshape(yc,m+2,n+2)'; u = ue(X,Y); l = u(:,1); r = u(:,end); @@ -164,21 +162,19 @@ rsC = [1 1+dy(i)/2:dy(i):2-dy(i)/2 2]; tsC = [pi pi-dx(i)/2:-dx(i):dx(i)/2 0]; [TSC,RSC] = meshgrid(tsC,rsC); - xc = RSC .* cos(TSC); xc = reshape(xc',[],1); - yc = RSC .* sin(TSC); yc = reshape(yc',[],1); + X = RSC .* cos(TSC); + Y = RSC .* sin(TSC); % Build operators curG = grad2DCurv(k,xn,yn); curD = div2DCurv(k,xn,yn); curL = curD * curG; - newG = grad2DCurv(k,xc,yc,m(i),dx(i),n(i),dy(i),dc,nc); - newD = div2DCurv(k,xc,yc,m(i),dx(i),n(i),dy(i),dc,nc); + newG = grad2DCurv(k,X,Y,dc,nc); + newD = div2DCurv(k,X,Y,dc,nc); newL = newD * newG; % Boundary Conditions - X = reshape(xc,m(i)+2,n(i)+2)'; % Reshape for plotting and easier BC - Y = reshape(yc,m(i)+2,n(i)+2)'; u = ue(X,Y); l = u(:,1); r = u(:,end); diff --git a/examples/matlab_octave/elliptic2DCurvPeriodic.m b/examples/matlab_octave/elliptic2DCurvPeriodic.m index 2060a473f..5b43ebbf5 100644 --- a/examples/matlab_octave/elliptic2DCurvPeriodic.m +++ b/examples/matlab_octave/elliptic2DCurvPeriodic.m @@ -31,20 +31,15 @@ rs = [1 (1 + dy / 2) : dy : (2 - dy / 2) 2]; ts = 0 : dx : (2*pi - dx); [T,R] = meshgrid(ts, rs); -xc = R .* cos(T); -yc = R .* sin(T); -xc = reshape(xc', [], 1); -yc = reshape(yc', [], 1); +X = R .* cos(T); +Y = R .* sin(T); % Build Operators -G = grad2DCurv(k, xc, yc, m, dx, n, dy, dc, nc); -D = div2DCurv(k, xc, yc, m, dx, n, dy, dc, nc); +G = grad2DCurv(k, X, Y, dc, nc); +D = div2DCurv(k, X, Y, dc, nc); L = D * G; % Boundary Conditions -X = reshape(xc, m, n + 2)'; % Reshape for plotting and easier BC -Y = reshape(yc, m, n + 2)'; - u = ue(X, Y); l = 0; % Left and right boundaries don't exist @@ -52,7 +47,7 @@ b = u(1, :)'; t = u(end, :)'; v = {l; r; b; t}; -B = zeros(size(xc)); +B = zeros(numel(X), 1); [L0, B0] = addScalarBC2D(L, B, k, m, dx, n, dy, dc, nc, v); ua = L0 \ B0; diff --git a/examples/matlab_octave/test_div2DCurv.m b/examples/matlab_octave/test_div2DCurv.m index cbd5f01fc..20d4dfc61 100644 --- a/examples/matlab_octave/test_div2DCurv.m +++ b/examples/matlab_octave/test_div2DCurv.m @@ -55,8 +55,11 @@ legend("Nodal Points", "u", "v", "Centers") hold off +X = reshape(Cx, m, n + 2)'; +Y = reshape(Cy, m, n + 2)'; + tic -D = div2DCurv(k, Cx, Cy, m, dx, n, dy, dc, nc); +D = div2DCurv(k, X, Y, dc, nc); toc % Vector Field @@ -72,10 +75,8 @@ Ccomp = reshape(Ccomp, m, n + 2)'; % Remove boundary points (divergence returns 0 on boundaries) -Cx = reshape(Cx, m, n + 2)'; -Cy = reshape(Cy, m, n + 2)'; -Cx = Cx(2:end-1, :); -Cy = Cy(2:end-1, :); +Cx = X(2:end-1, :); +Cy = Y(2:end-1, :); C = C(2:end-1, :); Ccomp = Ccomp(2:end-1, :); diff --git a/examples/matlab_octave/test_div3DCurv.m b/examples/matlab_octave/test_div3DCurv.m index 2186f3b3c..87e1ab7c7 100644 --- a/examples/matlab_octave/test_div3DCurv.m +++ b/examples/matlab_octave/test_div3DCurv.m @@ -16,9 +16,9 @@ dc = [1; 1; 1; 1; 1; 1]; nc = [0; 0; 0; 0; 0; 0]; -xlength = linspace(1,m,m); -ylength = linspace(1,n,n); -zlength = linspace(1,o,o); +% xlength = linspace(1,m,m); +% ylength = linspace(1,n,n); +% zlength = linspace(1,o,o); % [X, Y, Z] = meshgrid(xlength,ylength,zlength); @@ -47,11 +47,23 @@ yc = NtoC * reshape(Y, [], 1); zc = NtoC * reshape(Z, [], 1); +X = reshape(xc, m + 1, n + 1, o + 1); +Y = reshape(yc, m + 1, n + 1, o + 1); +Z = reshape(zc, m + 1, n + 1, o + 1); + +X = permute(X, [2, 1, 3]); +Y = permute(Y, [2, 1, 3]); +Z = permute(Z, [2, 1, 3]); + U = CtoU * sin(xc); V = CtoV * zeros(size(yc)); W = CtoW * zeros(size(zc)); -D = div3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); +D = div3DCurv(k, X, Y, Z, dc, nc); + +X = permute(X, [2, 1, 3]); +Y = permute(Y, [2, 1, 3]); +Z = permute(Z, [2, 1, 3]); % Apply the operator to the field Ccomp = D*[U; V; W]; @@ -60,10 +72,6 @@ Ccomp = reshape(Ccomp, m+1, n+1, o+1); Ccomp = Ccomp(2:end-1, 2:end-1, 2:end-1); -X = reshape(xc, m + 1, n + 1, o + 1); -Y = reshape(yc, m + 1, n + 1, o + 1); -Z = reshape(zc, m + 1, n + 1, o + 1); - X = X(2:end-1, 2:end-1, 2:end-1); Y = Y(2:end-1, 2:end-1, 2:end-1); Z = Z(2:end-1, 2:end-1, 2:end-1); diff --git a/examples/matlab_octave/test_grad2DCurv.m b/examples/matlab_octave/test_grad2DCurv.m index 454cc55be..313473eb6 100644 --- a/examples/matlab_octave/test_grad2DCurv.m +++ b/examples/matlab_octave/test_grad2DCurv.m @@ -4,7 +4,7 @@ % u = x^2 + y^2 % Gu = <2x, 2y> % - +clear clc close all @@ -29,15 +29,16 @@ set(gcf, 'Color', 'w') % Staggered logical grid -[Xs, Ys] = meshgrid([1 1.5 : 1 : m-0.5 m], [1 1.5 : 1 : n-0.5 n]); -Cs = Xs.^2 + Ys.^2; -% Reshape the field so it can be multiplied by the operator later on -C_ = reshape(Cs', [], 1); -Xs = reshape(Xs', [], 1); -Ys = reshape(Ys', [], 1); +NtoC = interpolNodesToCenters2D(k, m - 1, n - 1, dc, nc); +xc = NtoC * reshape(X', [], 1); +yc = NtoC * reshape(Y', [], 1); +C_ = xc.^2 + yc.^2; +Xs = reshape(xc, m + 1, n + 1)'; +Ys = reshape(yc, m + 1, n + 1)'; +Cs = reshape(C_, m + 1, n + 1)'; % Get 2D curvilinear mimetic gradient -G = grad2DCurv(k, Xs, Ys, m - 1, dx, n - 1, dy, dc, nc); +G = grad2DCurv(k, Xs, Ys, dc, nc); % Apply the operator to the field TMP = G*C_; @@ -53,10 +54,10 @@ CtoV = interpolCentersToFacesD1D(k, n - 1); CtoV = kron(CtoV, speye(m + 1)); -Ux = CtoU * Xs; -Uy = CtoU * Ys; -Vx = CtoV * Xs; -Vy = CtoV * Ys; +Ux = CtoU * xc; +Uy = CtoU * yc; +Vx = CtoV * xc; +Vy = CtoV * yc; Ux = reshape(Ux, m + 1, n)'; Uy = reshape(Uy, m + 1, n)'; @@ -67,7 +68,7 @@ figure set(gcf, 'Color', 'w') subplot(3, 1, 1) -surf(reshape(Xs, m + 1, n + 1)', reshape(Ys, m + 1, n + 1)', Cs, 'EdgeColor', 'none'); +surf(Xs, Ys, Cs, 'EdgeColor', 'none'); colorbar xlabel('x') ylabel('y') @@ -97,7 +98,7 @@ figure set(gcf, 'Color', 'w') subplot(3, 1, 1) -surf(reshape(Xs, m + 1, n + 1)', reshape(Ys, m + 1, n + 1)', Cs, 'EdgeColor', 'none'); +surf(Xs, Ys, Cs, 'EdgeColor', 'none'); colorbar xlabel('x') ylabel('y') diff --git a/examples/matlab_octave/test_grad3DCurv.m b/examples/matlab_octave/test_grad3DCurv.m index aab7afd2b..9341ee1e5 100644 --- a/examples/matlab_octave/test_grad3DCurv.m +++ b/examples/matlab_octave/test_grad3DCurv.m @@ -42,6 +42,12 @@ yc = INC * yn; zc = INC * zn; C = xc.^2 + yc.^2 + zc.^2; +X = reshape(xc, m + 1, n + 1, o + 1); +Y = reshape(yc, m + 1, n + 1, o + 1); +Z = reshape(zc, m + 1, n + 1, o + 1); +X = permute(X, [2, 1, 3]); +Y = permute(Y, [2, 1, 3]); +Z = permute(Z, [2, 1, 3]); Ux = ICFx * xc; Uy = ICFy * yc; @@ -64,7 +70,7 @@ axis equal % Get 3D curvilinear mimetic gradient -G = grad3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); +G = grad3DCurv(k, X, Y, Z, dc, nc); % Apply the operator to the field TMP = G*C; diff --git a/examples/matlab_octave/test_lap3DCurv.m b/examples/matlab_octave/test_lap3DCurv.m index e00775b08..e80a44953 100644 --- a/examples/matlab_octave/test_lap3DCurv.m +++ b/examples/matlab_octave/test_lap3DCurv.m @@ -28,13 +28,17 @@ yc = INC * reshape(permute(Y, [2, 1, 3]), [], 1); zc = INC * reshape(permute(Z, [2, 1, 3]), [], 1); +X = permute(reshape(xc, m + 1, n + 1, o + 1), [2, 1, 3]); +Y = permute(reshape(yc, m + 1, n + 1, o + 1), [2, 1, 3]); +Z = permute(reshape(zc, m + 1, n + 1, o + 1), [2, 1, 3]); + ue = xc.^2 + yc.^2 + zc.^2; RHS = ue; % Get 3D curvilinear mimetic divergence -D = div3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); +D = div3DCurv(k, X, Y, Z, dc, nc); % Get 3D curvilinear mimetic gradient -G = grad3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); +G = grad3DCurv(k, X, Y, Z, dc, nc); % Dirichlet BCs BC = robinBC3D(k, m-1, 1, n-1, 1, o-1, 1, 1, 0); % Laplacian operator with BCs diff --git a/examples/matlab_octave/test_lap3DCurv_case2.m b/examples/matlab_octave/test_lap3DCurv_case2.m index 8c62bc3bf..760230a58 100644 --- a/examples/matlab_octave/test_lap3DCurv_case2.m +++ b/examples/matlab_octave/test_lap3DCurv_case2.m @@ -32,13 +32,17 @@ yc = INC * reshape(permute(Y, [2, 1, 3]), [], 1); zc = INC * reshape(permute(Z, [2, 1, 3]), [], 1); +X = permute(reshape(xc, m + 1, n + 1, o + 1), [2, 1, 3]); +Y = permute(reshape(yc, m + 1, n + 1, o + 1), [2, 1, 3]); +Z = permute(reshape(zc, m + 1, n + 1, o + 1), [2, 1, 3]); + ue = -(xc.^2 + yc.^2 + zc.^2); % Unknown function RHS = ue; % Get 3D curvilinear mimetic divergence -D = div3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); +D = div3DCurv(k, X, Y, Z, dc, nc); % Get 3D curvilinear mimetic gradient -G = grad3DCurv(k, xc, yc, zc, m - 1, dx, n - 1, dy, o - 1, dz, dc, nc); +G = grad3DCurv(k, X, Y, Z, dc, nc); % Dirichlet BCs BC = robinBC3D(k, m-1, 1, n-1, 1, o-1, 1, 1, 0); % Laplacian operator with BCs diff --git a/src/matlab_octave/div2DCurv.m b/src/matlab_octave/div2DCurv.m index 3a907d1f6..5ec6ab203 100644 --- a/src/matlab_octave/div2DCurv.m +++ b/src/matlab_octave/div2DCurv.m @@ -1,4 +1,4 @@ -function D = div2DCurv(k, X, Y, m, dx, n, dy, dc, nc) +function D = div2DCurv(k, X, Y, dc, nc) % PURPOSE % Returns a 2D curvilinear mimetic divergence operator. If optional % arguments are specified, it acts on the extended faces @@ -11,10 +11,6 @@ % optional arguments are specified, else nodes % Y : y-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes -% (optional) m : Number of cells in xi direction -% (optional) dx : Step size in xi direction -% (optional) n : Number of cells in eta direction -% (optional) dy : Step size in eta direction % (optional) dc : a0 (4x1 vector for left, right, bottom, top % boundaries, resp.) % (optional) nc : b0 (4x1 vector for left, right, bottom, top @@ -22,7 +18,7 @@ % % SYNTAX % D = div2DCurv(k, X, Y) -% D = div2DCurv(k, X, Y, m, dx, n, dy, dc, nc) +% D = div2DCurv(k, X, Y, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -30,9 +26,9 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - if nargin ~= 3 && nargin ~= 9 + if nargin ~= 3 && nargin ~= 5 error("div2DCurv:InvalidNumArgs", ... - "div2DCurv expects 3 or 9 arguments") + "div2DCurv expects 3 or 5 arguments") elseif nargin == 3 D = div2DCurvLegacy(k,X,Y); return; @@ -41,14 +37,24 @@ assert(all(size(dc) == [4 1]), "dc is a 4x1 vector") assert(all(size(nc) == [4 1]), "nc is a 4x1 vector") + assert(size(X, 2) ~= 1, "X must be in matrix form") + assert(size(Y, 2) ~= 1, "Y must be in matrix form") + + assert(all(size(X) == size(Y)), "X and Y must be the same size") + + [n, m] = size(X); + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2,1)) + dx = 1 / m; De = divPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1DPeriodic(k,m); ICFx = interpolCentersToFacesD1DPeriodic(k,m); Im = speye(m); Bx = sparse(m, m); else + m = m - 2; + dx = 1 / m; De = divNonPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1D(k,m); ICFx = interpolCentersToFacesD1D(k,m); @@ -58,12 +64,15 @@ Bx(end, end) = 1; end if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) + dy = 1 / n; Dn = divPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1DPeriodic(k,n); ICFy = interpolCentersToFacesD1DPeriodic(k,n); In = speye(n); By = sparse(n, n); else + n = n - 2; + dy = 1 / n; Dn = divNonPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1D(k,n); ICFy = interpolCentersToFacesD1D(k,n); @@ -79,7 +88,7 @@ Dn = kron(Dn,Im) * kron(ICFy,Im); % Apply metrics - [J,Xe,Xn,Ye,Yn] = jacobian2D(k,X,Y,m,dx,n,dy,dc,nc); + [J,Xe,Xn,Ye,Yn] = jacobian2D(k,X,Y,dc,nc); Dx = (Yn ./ J) .* De - (Ye ./ J) .* Dn; Dy = (Xe ./ J) .* Dn - (Xn ./ J) .* De; diff --git a/src/matlab_octave/div3DCurv.m b/src/matlab_octave/div3DCurv.m index 69d5496d1..cc50b621c 100644 --- a/src/matlab_octave/div3DCurv.m +++ b/src/matlab_octave/div3DCurv.m @@ -1,4 +1,4 @@ -function D = div3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +function D = div3DCurv(k, X, Y, Z, dc, nc) % PURPOSE % Returns a 3D curvilinear mimetic divergence operator. If optional % arguments are specified, it acts on the extended faces @@ -13,12 +13,6 @@ % optional arguments are specified, else nodes % Z : z-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes -% (optional) m : Number of cells in xi direction -% (optional) dx : Step size in xi direction -% (optional) n : Number of cells in eta direction -% (optional) dy : Step size in eta direction -% (optional) o : Number of cells in kappa direction -% (optional) dz : Step size in kappa direction % (optional) dc : a0 (6x1 vector for left, right, bottom, top, front, and % back boundaries, resp.) % (optional) nc : b0 (6x1 vector for left, right, bottom, top, front, and @@ -26,7 +20,7 @@ % % SYNTAX % D = div3DCurv(k, X, Y, Z) -% D = div3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +% D = div3DCurv(k, X, Y, Z, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -34,9 +28,9 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - if nargin ~= 4 && nargin ~= 12 + if nargin ~= 4 && nargin ~= 6 error("div3DCurv:InvalidNumArgs", ... - "div3DCurv expects 4 or 12 arguments") + "div3DCurv expects 4 or 6 arguments") elseif nargin == 4 D = div3DCurvLegacy(k, X, Y, Z); return; @@ -45,14 +39,25 @@ assert(all(size(dc) == [6 1]), "dc is a 6x1 vector") assert(all(size(nc) == [6 1]), "nc is a 6x1 vector") + assert(size(X, 2) ~= 1, "X must be in matrix form") + assert(size(Y, 2) ~= 1, "Y must be in matrix form") + assert(size(Z, 2) ~= 1, "Z must be in matrix form") + + assert(all(size(X) == size(Y)) && all(size(X) == size(Z)), "X, Y, and Z must all be the same size") + + [n, m, o] = size(X); + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2,1)) + dx = 1 / m; De = divPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1DPeriodic(k,m); ICFx = interpolCentersToFacesD1DPeriodic(k,m); Im = speye(m); Bx = sparse(m, m); else + m = m - 2; + dx = 1 / m; De = divNonPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1D(k,m); ICFx = interpolCentersToFacesD1D(k,m); @@ -62,12 +67,15 @@ Bx(end, end) = 1; end if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) + dy = 1 / n; Dn = divPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1DPeriodic(k,n); ICFy = interpolCentersToFacesD1DPeriodic(k,n); In = speye(n); By = sparse(n, n); else + n = n - 2; + dy = 1 / n; Dn = divNonPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1D(k,n); ICFy = interpolCentersToFacesD1D(k,n); @@ -77,12 +85,15 @@ By(end, end) = 1; end if isempty(find(dc(5:6).^2 + nc(5:6).^2, 1)) + dz = 1 / o; Dk = divPeriodic(k, o, dz); IFCz = interpolFacesToCentersG1DPeriodic(k, o); ICFz = interpolCentersToFacesD1DPeriodic(k, o); Io = speye(o); Bz = sparse(o, o); else + o = o - 2; + dz = 1 / o; Dk = divNonPeriodic(k, o, dz); IFCz = interpolFacesToCentersG1D(k, o); ICFz = interpolCentersToFacesD1D(k, o); @@ -99,7 +110,7 @@ Dk = kron(kron(Dk, In), Im) * kron(kron(ICFz, In), Im); % Apply metrics - [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc); + [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, dc, nc); Dx = ((Yn .* Zk - Zn .* Yk) ./ J) .* De ... + ((Ze .* Yk - Ye .* Zk) ./ J) .* Dn ... diff --git a/src/matlab_octave/grad2DCurv.m b/src/matlab_octave/grad2DCurv.m index 49335821a..ee23ca5a5 100644 --- a/src/matlab_octave/grad2DCurv.m +++ b/src/matlab_octave/grad2DCurv.m @@ -1,4 +1,4 @@ -function G = grad2DCurv(k, X, Y, m, dx, n, dy, dc, nc) +function G = grad2DCurv(k, X, Y, dc, nc) % PURPOSE % Returns a 2D curvilinear mimetic gradient operator. If optional % arguments are specified, it outputs to the extended @@ -11,10 +11,6 @@ % optional arguments are specified, else nodes % Y : y-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes -% (optional) m : Number of cells in xi direction -% (optional) dx : Step size in xi direction -% (optional) n : Number of cells in eta direction -% (optional) dy : Step size in eta direction % (optional) dc : a0 (4x1 vector for left, right, bottom, top % boundaries, resp.) % (optional) nc : b0 (4x1 vector for left, right, bottom, top @@ -22,7 +18,7 @@ % % SYNTAX % G = grad2DCurv(k, X, Y) -% G = grad2DCurv(k, X, Y, m, dx, n, dy, dc, nc) +% G = grad2DCurv(k, X, Y, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -30,35 +26,48 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - if nargin ~=3 && nargin ~= 9 + if nargin ~= 3 && nargin ~= 5 error("grad2DCurv:InvalidNumArgs", ... - "grad2DCurv expects 3 or 9 arguments") + "grad2DCurv expects 3 or 5 arguments") elseif nargin == 3 - G = grad2DCurvLegacy(k,X,Y); + G = grad2DCurvLegacy(k, X, Y); return; end assert(all(size(dc) == [4 1]), "dc is a 4x1 vector") assert(all(size(nc) == [4 1]), "nc is a 4x1 vector") + assert(size(X, 2) ~= 1, "X must be in matrix form") + assert(size(Y, 2) ~= 1, "Y must be in matrix form") + + assert(all(size(X) == size(Y)), "X and Y must be the same size") + + [n, m] = size(X); + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) + dx = 1 / m; Ge = gradPeriodic(k,m,dx); ICFx = interpolCentersToFacesD1DPeriodic(k,m); IFCx = interpolFacesToCentersG1DPeriodic(k,m); Im = speye(m); else + m = m - 2; + dx = 1 / m; Ge = gradNonPeriodic(k,m,dx); ICFx = interpolCentersToFacesD1D(k,m); IFCx = interpolFacesToCentersG1D(k,m); Im = speye(m+2); end if isempty(find(dc(3:4).^2 + nc(3:4).^2, 1)) + dy = 1 / n; Gn = gradPeriodic(k,n,dy); ICFy = interpolCentersToFacesD1DPeriodic(k,n); IFCy = interpolFacesToCentersG1DPeriodic(k,n); In = speye(n); else + n = n - 2; + dy = 1 / n; Gn = gradNonPeriodic(k,n,dy); ICFy = interpolCentersToFacesD1D(k,n); IFCy = interpolFacesToCentersG1D(k,n); @@ -71,7 +80,7 @@ Gn = kron(IFCy,Im) * kron(Gn,Im); % Apply Metrics - [J,Xe,Xn,Ye,Yn] = jacobian2D(k,X,Y,m,dx,n,dy,dc,nc); + [J,Xe,Xn,Ye,Yn] = jacobian2D(k,X,Y,dc,nc); Gx = (Yn ./ J) .* Ge - (Ye ./ J) .* Gn; Gy = (Xe ./ J) .* Gn - (Xn ./ J) .* Ge; diff --git a/src/matlab_octave/grad3DCurv.m b/src/matlab_octave/grad3DCurv.m index a3c0322b3..82a0a0e9f 100644 --- a/src/matlab_octave/grad3DCurv.m +++ b/src/matlab_octave/grad3DCurv.m @@ -1,4 +1,4 @@ -function G = grad3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +function G = grad3DCurv(k, X, Y, Z, dc, nc) % PURPOSE % Returns a 3D curvilinear mimetic gradient operator. If optional % arguments are specified, it outputs to the extended @@ -13,12 +13,6 @@ % optional arguments are specified, else nodes % Z : Z-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes -% (optional) m : Number of cells in xi direction -% (optional) dx : Step size in xi direction -% (optional) n : Number of cells in eta direction -% (optional) dy : Step size in eta direction -% (optional) o : Number of cells in kappa direction -% (optional) dz : Step size in kappa direction % (optional) dc : a0 (6x1 vector for left, right, bottom, top, front, and % back boundaries, resp.) % (optional) nc : b0 (6x1 vector for left, right, bottom, top, front, and @@ -26,7 +20,7 @@ % % SYNTAX % G = grad3DCurv(k, X, Y, Z) -% G = grad3DCurv(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +% G = grad3DCurv(k, X, Y, Z, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -34,9 +28,9 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - if nargin ~= 4 && nargin ~= 12 + if nargin ~= 4 && nargin ~= 6 error("grad3DCurv:InvalidNumArgs", ... - "grad3DCurv expects 4 or 12 arguments") + "grad3DCurv expects 4 or 6 arguments") elseif nargin == 4 G = grad3DCurvLegacy(k, X, Y, Z); return; @@ -45,35 +39,52 @@ assert(all(size(dc) == [6 1]), "dc is a 6x1 vector") assert(all(size(nc) == [6 1]), "nc is a 6x1 vector") + assert(size(X, 2) ~= 1, "X must be in matrix form") + assert(size(Y, 2) ~= 1, "Y must be in matrix form") + assert(size(Z, 2) ~= 1, "Z must be in matrix form") + + assert(all(size(X) == size(Y)) && all(size(X) == size(Z)), "X, Y, and Z must all be the same size") + + [n, m, o] = size(X); + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) + dx = 1 / m; Ge = gradPeriodic(k,m,dx); ICFx = interpolCentersToFacesD1DPeriodic(k,m); IFCx = interpolFacesToCentersG1DPeriodic(k,m); Im = speye(m); else + m = m - 2; + dx = 1 / m; Ge = gradNonPeriodic(k,m,dx); ICFx = interpolCentersToFacesD1D(k,m); IFCx = interpolFacesToCentersG1D(k,m); Im = speye(m+2); end if isempty(find(dc(3:4).^2 + nc(3:4).^2, 1)) + dy = 1 / n; Gn = gradPeriodic(k,n,dy); ICFy = interpolCentersToFacesD1DPeriodic(k,n); IFCy = interpolFacesToCentersG1DPeriodic(k,n); In = speye(n); else + n = n - 2; + dy = 1 / n; Gn = gradNonPeriodic(k,n,dy); ICFy = interpolCentersToFacesD1D(k,n); IFCy = interpolFacesToCentersG1D(k,n); In = speye(n+2); end if isempty(find(dc(5:6).^2 + nc(5:6).^2, 1)) + dz = 1 / o; Gk = gradPeriodic(k, o, dz); ICFz = interpolCentersToFacesD1DPeriodic(k,o); IFCz = interpolFacesToCentersG1DPeriodic(k,o); Io = speye(o); else + o = o - 2; + dz = 1 / o; Gk = gradNonPeriodic(k,o,dz); ICFz = interpolCentersToFacesD1D(k,o); IFCz = interpolFacesToCentersG1D(k,o); @@ -87,7 +98,7 @@ Gk = kron(kron(IFCz, In), Im) * kron(kron(Gk, In), Im); % Apply metrics - [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc); + [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, dc, nc); Gx = ((Yn .* Zk - Zn .* Yk) ./ J) .* Ge ... + ((Ze .* Yk - Ye .* Zk) ./ J) .* Gn ... diff --git a/src/matlab_octave/jacobian2D.m b/src/matlab_octave/jacobian2D.m index 61ad85c05..c4f446656 100644 --- a/src/matlab_octave/jacobian2D.m +++ b/src/matlab_octave/jacobian2D.m @@ -1,4 +1,4 @@ -function [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, m, dx, n, dy, dc, nc) +function [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, dc, nc) % PURPOSE % Returns the Jacobian metrics (Xe, Xn, Ye, Yn, and J = XeYn - XnYe) of a % mesh @@ -10,10 +10,6 @@ % optional arguments are specified, else nodes % Y : y-coordinates (physical) of meshgrid centers if % optional arguments are specified, else nodes -% (optional) m : Number of cells in xi direction -% (optional) dx : Step size in xi direction -% (optional) n : Number of cells in eta direction -% (optional) dy : Step size in eta direction % (optional) dc : a0 (4x1 vector for left, right, bottom, top % boundaries, resp.) % (optional) nc : b0 (4x1 vector for left, right, bottom, top @@ -21,7 +17,7 @@ % % SYNTAX % [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y) -% [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, m, dx, n, dy, dc, nc) +% [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -29,9 +25,9 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - if nargin ~= 3 && nargin ~= 9 + if nargin ~= 3 && nargin ~= 5 error("jacobian2D:InvalidNumArgs", ... - "jacobian2D expects 3 or 9 arguments") + "jacobian2D expects 3 or 5 arguments") elseif nargin == 3 [J,Xe,Xn,Ye,Yn] = jacobian2DLegacy(k,X,Y); return; @@ -40,27 +36,42 @@ assert(all(size(dc) == [4 1]), "dc is a 4x1 vector") assert(all(size(nc) == [4 1]), "nc is a 4x1 vector") + assert(size(X, 2) ~= 1, "X must be in matrix form") + assert(size(Y, 2) ~= 1, "Y must be in matrix form") + + assert(all(size(X) == size(Y)), "X and Y must be the same size") + + [n, m] = size(X); + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2,1)) + dx = 1 / m; Ge = gradPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1DPeriodic(k,m); Im = speye(m); else + m = m - 2; + dx = 1 / m; Ge = gradNonPeriodic(k,m,dx); IFCx = interpolFacesToCentersG1D(k,m); Im = speye(m+2); end if isempty(find(dc(3:4).^2 + nc(3:4).^2,1)) + dy = 1 / n; Gn = gradPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1DPeriodic(k,n); In = speye(n); else + n = n - 2; + dy = 1 / n; Gn = gradNonPeriodic(k,n,dy); IFCy = interpolFacesToCentersG1D(k,n); In = speye(n+2); end - numC = numel(X); + xc = reshape(X', [], 1); + yc = reshape(Y', [], 1); + numC = numel(xc); % Get metrics on centers % We don't augment the identity matrices @@ -70,7 +81,7 @@ Gn = kron(IFCy,Im) * kron(Gn,Im); G = [Ge; Gn]; - metrics = G * [X Y]; + metrics = G * [xc yc]; Xe = metrics(1:numC,1); Xn = metrics(1+numC:end,1); diff --git a/src/matlab_octave/jacobian3D.m b/src/matlab_octave/jacobian3D.m index 0f8544b21..22c745f45 100644 --- a/src/matlab_octave/jacobian3D.m +++ b/src/matlab_octave/jacobian3D.m @@ -1,4 +1,4 @@ -function [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +function [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, dc, nc) % PURPOSE % Returns the 3D jacobian metrics of a mesh % @@ -11,12 +11,6 @@ % arguments are specified, else nodes % Z : z-coordinates (physical) of meshgrid centers if optional % arguments are specified, else nodes -% (optional) m : Number of cells in xi direction -% (optional) dx : Step size in xi direction -% (optional) n : Number of cells in eta direction -% (optional) dy : Step size in eta direction -% (optional) o : Number of cells in kappa direction -% (optional) dz : Step size in kappa direction % (optional) dc : a0 (6x1 vector for left, right, bottom, top, front, and % back boundaries, resp.) % (optional) nc : b0 (6x1 vector for left, right, bottom, top, front, and @@ -24,7 +18,7 @@ % % SYNTAX % [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z) -% [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, m, dx, n, dy, o, dz, dc, nc) +% [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3D(k, X, Y, Z, dc, nc) % % ---------------------------------------------------------------------------- % SPDX-License-Identifier: GPL-3.0-or-later @@ -32,9 +26,9 @@ % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. % ---------------------------------------------------------------------------- - if nargin ~= 4 && nargin ~= 12 + if nargin ~= 4 && nargin ~= 6 error("jacobian3D:InvalidNumArgs", ... - "jacobian3D expects 4 or 12 arguments") + "jacobian3D expects 4 or 6 arguments") elseif nargin == 4 [J, Xe, Xn, Xk, Ye, Yn, Yk, Ze, Zn, Zk] = jacobian3DLegacy(k, X, Y, Z); return; @@ -43,36 +37,56 @@ assert(all(size(dc) == [6 1]), "dc is a 6x1 vector") assert(all(size(nc) == [6 1]), "nc is a 6x1 vector") + assert(size(X, 2) ~= 1, "X must be in matrix form") + assert(size(Y, 2) ~= 1, "Y must be in matrix form") + assert(size(Z, 2) ~= 1, "Z must be in matrix form") + + assert(all(size(X) == size(Y)) && all(size(X) == size(Z)), "X, Y, and Z must all be the same size") + + [n, m, o] = size(X); + % Periodic Handling if isempty(find(dc(1:2).^2 + nc(1:2).^2, 1)) + dx = 1 / m; Ge = gradPeriodic(k, m, dx); IFCx = interpolFacesToCentersG1DPeriodic(k, m); Im = speye(m); else + m = m - 2; + dx = 1 / m; Ge = gradNonPeriodic(k, m, dx); IFCx = interpolFacesToCentersG1D(k, m); Im = speye(m + 2); end if isempty(find(dc(3:4).^2 + nc(3:4).^2, 1)) + dy = 1 / n; Gn = gradPeriodic(k, n, dy); IFCy = interpolFacesToCentersG1DPeriodic(k, n); In = speye(n); else + n = n - 2; + dy = 1 / n; Gn = gradNonPeriodic(k, n, dy); IFCy = interpolFacesToCentersG1D(k, n); In = speye(n + 2); end if isempty(find(dc(5:6).^2 + nc(5:6).^2, 1)) + dz = 1 / o; Gk = gradPeriodic(k, o, dz); IFCz = interpolFacesToCentersG1DPeriodic(k, o); Io = speye(o); else + o = o - 2; + dz = 1 / o; Gk = gradNonPeriodic(k, o, dz); IFCz = interpolFacesToCentersG1D(k, o); Io = speye(o + 2); end - numC = numel(X); + xc = reshape(permute(X, [2 1 3]), [], 1); + yc = reshape(permute(Y, [2 1 3]), [], 1); + zc = reshape(permute(Z, [2 1 3]), [], 1); + numC = numel(xc); % Get metrics on centers % We don't augment the identity matrices @@ -83,7 +97,7 @@ Gk = kron(kron(IFCz, In), Im) * kron(kron(Gk, In), Im); G = [Ge; Gn; Gk]; - metrics = G * [X Y Z]; + metrics = G * [xc yc zc]; Xe = metrics(1:numC, 1); Xn = metrics(1+numC:2*numC, 1);