1+ function [gradE , E , resid ] = fasthelmholtz(grid ,V )
2+
3+ % Recreate the original spacing vectors defining the grid
4+ gridSpacing = extractGridSpacing(grid );
5+
6+ % % Convert from ndgrid to meshgrid
7+ % dimcount = numel(grid); % Count how many dimensions field is defined over
8+ % dimlist = [2 1 3:dimcount]; % Permutation mapping from ndgrid to meshgrid
9+ % grid = permute(grid,dimlist); % Perform permutation on grid component ordering
10+ % V = permute(V,dimlist); % Perform permutation on vector component ordering
11+ % gridSpacing = permute(gridSpacing,dimlist); % Perform permutation on grid spacing component ordering
12+ % %grid = cellfun(@(g) permute(g,dimlist),grid,'UniformOutput',false); % Perform permutation on grid value ordering
13+ % %V = cellfun(@(g) permute(g,dimlist),V,'UniformOutput',false); % Perform permutation on vector value ordering
14+
15+
16+ grid = {grid{2 }; grid{1 }; grid{3 : end }};
17+ V = {V{2 }; V{1 }; V{3 : end }};
18+
19+
20+ % %%%
21+ % Take the gradient of each vector component
22+
23+ % Make a cell matrix the same size as the vector cell matrix
24+ gradV = cell(size(V ));
25+
26+ % Loop over each element of this cell (one cell per component of the
27+ % vector)
28+ for idx = 1 : numel(grid )
29+
30+ % Each vector component has a derivative along each vector component
31+ gradV{idx } = cell(size(grid ));
32+
33+ % Calculate the gradient
34+ [gradV{idx }{: }] = gradient(V{idx },gridSpacing{: });
35+ end
36+
37+
38+ % %%%%
39+ % Convert the gradients into gradients calculated at each location
40+
41+ % One matrix at each grid point, with as many rows/columns as there are
42+ % dimensions
43+ gradV_local = repmat({zeros(numel(grid ),numel(grid ))},size(grid{1 }));
44+
45+ for idx1 = 1 : numel(grid ) % Loop over all vector components
46+ for idx2 = 1 : numel(grid ) % Loop over all derivative directions
47+ for idx3 = 1 : numel(grid{1 }) % Loop over all grid points
48+
49+ gradV_local{idx3 }(idx1 ,idx2 ) = gradV{idx1 }{idx2 }(idx3 );
50+
51+ end
52+ end
53+ end
54+
55+
56+ % %%%%%%%%%
57+ % Remove the skew-symmetry in the gradient, remainder is derivative of
58+ % conservative function
59+
60+ gradV_conservative = gradV ;
61+
62+ for idx1 = 1 : numel(grid ) % Loop over all vector components
63+ for idx2 = 1 : numel(grid ) % Loop over all derivative directions
64+ for idx3 = 1 : numel(grid{1 }) % Loop over all grid points
65+
66+ gradV_conservative{idx1 }{idx2 }(idx3 ) = ...
67+ (gradV{idx1 }{idx2 }(idx3 ) + ...
68+ gradV{idx2 }{idx1 }(idx3 ))/2 ;
69+
70+ end
71+ end
72+ end
73+
74+
75+ % %%%%%%
76+ % Integrate conservative component
77+ V_conservative = V ;
78+ for idx = 1 : numel(grid ) % Loop over all vector components
79+ V_conservative{idx } = intgrad2(gradV_conservative{idx }{: },gridSpacing{: });
80+ end
81+
82+ % %%%%%%
83+ % Subtract conservative component from original vector field to get the
84+ % residual non-conservative component
85+ V_nonconservative = V ;
86+ for idx = 1 : numel(grid ) % Loop over all vector components
87+ V_nonconservative{idx } = V{idx } - V_conservative{idx };
88+ end
89+
90+ % %%%%%%
91+ % Calculate the mean of the nonconservative field, so that we can extract
92+ % this harmonic component
93+ % residual non-conservative component
94+ V_nonconservative_mean = V ;
95+ for idx = 1 : numel(grid ) % Loop over all vector components
96+ V_nonconservative_mean{idx } = mean(V_nonconservative{idx }(: ));
97+ end
98+
99+ % %%%%%%%
100+ % Add the mean of the non-conservative field back into the conservative
101+ % field
102+ V_conservative_fixed = V ;
103+ V_nonconservative_fixed = V ;
104+ for idx = 1 : numel(grid ) % Loop over all vector components
105+ V_conservative_fixed{idx } = V_conservative{idx }+V_nonconservative_mean{idx };
106+ V_nonconservative_fixed{idx } = V_nonconservative{idx }-V_nonconservative_mean{idx };
107+ end
108+
109+
110+ %
111+ % %%%%
112+ % % Split out the conservative component of the vector gradient
113+ %
114+ % % Predefine a cell array to hold conservative components
115+ % gradV_local_conservative = gradV_local;
116+ %
117+ % for idx = 1:numel(gradV_local) % Loop over all grid points
118+ %
119+ % % Conservative component is average of gradient and its transpose
120+ % gradV_local_conservative{idx} = (gradV_local{idx} + transpose(gradV_local{idx}))/2;
121+ %
122+ % end
123+
124+
125+ % %%%%
126+ % Make the
127+
128+
129+ % % Convert from meshgrid to ndgrid
130+ V_conservative_fixed = {V_conservative_fixed{2 }; V_conservative_fixed{1 }; V_conservative_fixed{3 : end }};
131+ % V_conservative_fixed = cellfun(@(g) permute(g,dimlist),V_conservative_fixed,'UniformOutput',false); % Perform permutation on vector value ordering
132+
133+ V_nonconservative_fixed = {V_nonconservative_fixed{2 }; V_nonconservative_fixed{1 }; V_nonconservative_fixed{3 : end }};
134+ % V_nonconservative_fixed = cellfun(@(g) permute(g,dimlist),V_nonconservative_fixed,'UniformOutput',false); % Perform permutation on vector value ordering
135+
136+ E = [];
137+
138+ gradE = V_conservative_fixed ;
139+ resid = V_nonconservative_fixed ;
140+
141+
142+ end
0 commit comments