-
Notifications
You must be signed in to change notification settings - Fork 4
/
TriSystemTikhonovGradientEnergy.m
124 lines (92 loc) · 2.94 KB
/
TriSystemTikhonovGradientEnergy.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Smoothed Quadratic Energies on Meshes
%% J. Martinez Esturo, C. Rössl, and H. Theisel
%%
%% ACM Transactions on Graphics 2014
%%
%% Copyright J. Martinez Esturo 2014 (MIT License)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
classdef TriSystemTikhonovGradientEnergy < handle
%TRISYSTEMTIKHONOVGRADIENTENERGY
properties (SetAccess = protected)
AAs
bbs
E = []
b = []
end
properties (Access = protected)
en = []
d = []
nv = []
nt = []
beta = []
%for later total and local energy evaluation
An = []
AT = []
T = []
end
methods
function obj = TriSystemTikhonovGradientEnergy(mesh,beta, en, d,c)
%%
% mesh - triangle mesh to apply energy smoothing to
% (mesh.p3_0 is used for integration domain / vertex coordinates)
% beta - smoothing factor in [0,1]
% en - number of quadratic error components per triangle
% d - dimension of unknown coefficients at vertices
% c - number of independent functions to setup system rhs
% (defaults to 1)
if nargin < 6, c = 1; end
obj.beta = beta;
%% number of error components
obj.en = en;
obj.d = d;
obj.nv = mesh.nv;
obj.nt = mesh.nt;
%% compute triangle areas for replicated area weighting
obj.An = sparse(1:(en*obj.nt), 1:(en*obj.nt), repmat(mesh.A0,en,1));
%% regularization by gradient operator
% (we use 3D gradient vectors)
A3 = sparse(1:(3*obj.nt), 1:(3*obj.nt), repmat(mesh.A0,3,1));
G = mesh.GP_3;
perm = reshape(1:(d*obj.nv),d,[])';
P = sparse (1:(d*obj.nv),perm,1);
for i=1:d,
Tc{i} = G;
ATc{i} = A3;
end
obj.T = spblkdiag(Tc) * P;
obj.AT = spblkdiag(ATc);
%% init rhs if b == 0
obj.bbs = zeros(d*obj.nv,c);
obj.b = zeros(en*obj.nt,c);
end
function updateLHS(obj,E)
assert(size(E,1) == obj.nt*obj.en);
assert(size(E,2) == obj.nv*obj.d);
obj.E = E;
obj.AAs = (1-obj.beta) * obj.E'*obj.An*obj.E + ...
obj.beta * obj.T'*obj.AT*obj.T;
end
function updateRHS(obj,b)
assert(size(obj.E,1) == obj.nt*obj.en);
assert(size(obj.E,1) == size(b,1));
obj.b = b;
obj.bbs = (1-obj.beta) * obj.E'*obj.An*b;
end
function [terrEnergy,errEnergy,errSmoothn,errTotal] = ...
evalTriangleErrors(obj,u)
% local problem per triangle errors
e = obj.E * u - obj.b;
terrEnergy = sum(e.*e,2);
terrEnergy = sum(reshape(terrEnergy,obj.en,[]));
if nargout < 2, return; end;
% integrated local errors
errEnergy = trace(e'*obj.An*e);
if nargout < 3, return; end;
% integrated smoothness errors
errSmoothn = obj.T*u;
errSmoothn = trace(errSmoothn'*obj.AT*errSmoothn);
errTotal = errEnergy + errSmoothn;
end
end
end