convert the tensor of mixed volumes for application to hC-vectors THIS FUNCTION IS NO USER FUNCTION.
0001 function tensorHc = convertMeasureTensorSe(tensorFull, indexVector, ... 0002 mapping, dim) 0003 % convert the tensor of mixed volumes for application to hC-vectors 0004 % 0005 % THIS FUNCTION IS NO USER FUNCTION. 0006 0007 % The elk-library: convex geometry applied to crystallization modeling. 0008 % Copyright (C) 2012 Alexander Reinhold 0009 % 0010 % This program is free software: you can redistribute it and/or modify it 0011 % under the terms of the GNU General Public License as published by the 0012 % Free Software Foundation, either version 3 of the License, or (at your 0013 % option) any later version. 0014 % 0015 % This program is distributed in the hope that it will be useful, but 0016 % WITHOUT ANY WARRANTY; without even the implied warranty of 0017 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0018 % General Public License for more details. 0019 % 0020 % You should have received a copy of the GNU General Public License along 0021 % with this program. If not, see <http://www.gnu.org/licenses/>. 0022 0023 % substitute nan with zero 0024 tensorFull(isnan(tensorFull)) = 0; 0025 0026 switch dim 0027 case 1 0028 % mapping: nC x nC 0029 % tensorFull: 1 x nS 0030 % tensor: 1 x nC 0031 % tensorHc: 1 x nC 0032 tensor = tensorFull(indexVector); 0033 0034 % meas = sum(i): mixV(i) * lam(i) 0035 % meas = sum(i): mixV(i) * (sum(q): mapping(i, q) * hc(q)) 0036 % meas = sum(q): hc(j) * (sum(i): mapping(i, q) * mixV(i)) 0037 % meas = sum(q): hc(q) * tensor(q) 0038 tensorHc = tensor * mapping; 0039 0040 case 2 0041 % mapping: nC x nC 0042 % tensorFull: nS x nS 0043 % tensor: nC x nC 0044 % tensorHc: nC x nC 0045 tensor = tensorFull(indexVector, indexVector); 0046 nCol = size(mapping, 2); 0047 0048 % meas = sum(i,j): mixV(i,j) * lam(i) * lam(j) 0049 % meas = sum(i,j): mixV(i,j) * (sum(q): mapping(i, q) * hc(q)) 0050 % * (sum(r): mapping(j, r) * hc(r)) 0051 % meas = sum(q,r): hc(q) * hc(r) 0052 % * sum(i,j): mapping(i, q) * mixV(i,j) * mapping(j, r) 0053 % meas = sum(q,r): hc(q) * hc(r) * tensor(q,r) 0054 tensorHc = mapping' * tensor * mapping; 0055 0056 % to make the hC-tensors comparable, we must see that the 0057 % application is for i~=j the effective value is: 0058 % tensor(i,j) + tensor(j,i). 0059 % so that tensors with j<i are forced to 0 with the following 0060 filterLower = logical(tril(ones(nCol)) - eye(nCol)); 0061 tensorHc = tensorHc + (tensorHc.*filterLower)'; 0062 tensorHc(filterLower) = 0; 0063 0064 case 3 0065 % mapping: nC x nC 0066 % tensorFull: nS x nS x nS 0067 % tensor: nC x nC x nC 0068 % tensorHc: nC x nC x nC 0069 tensor = tensorFull(indexVector, indexVector, indexVector); 0070 nCol = size(mapping, 2); 0071 nRow = size(mapping, 1); 0072 0073 % meas = sum(i,j,k): mixV(i,j,k) * lam(i) * lam(j) * lam(k) 0074 % meas = sum(i,j,k): mixV(i,j,k) * (sum(q): mapping(i, q) * hc(q)) 0075 % * (sum(r): mapping(j, r) * hc(r)) 0076 % * (sum(s): mapping(k, s) * hc(s)) 0077 % meas = sum(q,r,s): hc(q) * hc(r) * hc(s) 0078 % * sum(k) mapping(k,s) * ( 0079 % sum(i,j): mapping(i, q) * mixV(i,j,k) * mapping(j, s) 0080 % meas = sum(q,r,s): hc(q) * hc(r) * hc(s) * tensor(q,r,s) 0081 % 0082 % where the inner most sum can be calculated as above, and is 0083 % summed up in a for loop 0084 0085 tensorHc = zeros([nCol, nCol, nCol]); 0086 for iResult = 1:nCol 0087 0088 tensorHcTmp = zeros([nCol, nCol, nCol]); 0089 0090 for iSum = 1:nRow 0091 tensorHcTmp(:,:,iSum) = mapping(iSum, iResult) * ... 0092 (mapping' * tensor(:,:,iSum) * mapping); 0093 end 0094 0095 tensorHc(:,:,iResult) = sum(tensorHcTmp, 3); 0096 end 0097 0098 % to make the hC-tensors comparable, we must see that the 0099 % effective coefficient for (i=j=k) is: 0100 % V(i,j,k) 0101 % and for (i = j, j < k) is: 0102 % V(i,j,k) + V(i,k,j) 0103 % and for (i < j, j = k) is: 0104 % V(i,j,k) + V(j,i,k) 0105 % and for (i < j < k) is: 0106 % V(i,j,k) + V(i,k,j) + V(j,i,k) + V(j,k,i) + V(k,i,j) + 0107 % V(k,j,i) 0108 0109 iTensor = ones([nCol, nCol, nCol]); 0110 jTensor = iTensor; 0111 kTensor = iTensor; 0112 for i = 1:nCol 0113 iTensor(i,:,:) = iTensor(i,:,:) * i; 0114 jTensor(:,i,:) = jTensor(:,i,:) * i; 0115 kTensor(:,:,i) = kTensor(:,:,i) * i; 0116 end 0117 0118 % We assign the values that should not be zero! 0119 % i = j = k (nothing to do) 0120 % i < j < k 0121 filter = (iTensor < jTensor) & (jTensor < kTensor); 0122 tensorHc = tensorHc ... 0123 + permute(tensorHc.*permute(filter, [1 3 2]), [1 3 2]) ... 0124 + permute(tensorHc.*permute(filter, [2 1 3]), [2 1 3]) ... 0125 + permute(tensorHc.*permute(filter, [3 1 2]), [2 3 1]) ... 0126 + permute(tensorHc.*permute(filter, [2 3 1]), [3 1 2]) ... 0127 + permute(tensorHc.*permute(filter, [3 2 1]), [3 2 1]); 0128 0129 % i = j, j < k 0130 filter = (iTensor == jTensor) & (jTensor < kTensor); 0131 tensorHc = tensorHc + ... 0132 + permute(tensorHc.*permute(filter, [1 3 2]), [1 3 2]) ... 0133 + permute(tensorHc.*permute(filter, [3 1 2]), [2 3 1]); 0134 0135 % j = k, i < j 0136 filter = (iTensor < jTensor) & (jTensor == kTensor); 0137 tensorHc = tensorHc + ... 0138 + permute(tensorHc.*permute(filter, [2 3 1]), [3 1 2]) ... 0139 + permute(tensorHc.*permute(filter, [2 1 3]), [2 1 3]); 0140 0141 % And we delete the values that are not permitted 0142 % i > j or j > k 0143 tensorHc((iTensor > jTensor) | (jTensor > kTensor)) = 0; 0144 0145 end 0146 0147 end