64 std::array<double,q> discIntegral;
65 for (
int iVec = 0; iVec < q; iVec++ ) {
66 discIntegral[iVec] = 0;
67 for (
int iSum = 0; iSum < q; iSum++ ) {
68 discIntegral[iVec] += descriptors::t<double,DESCRIPTOR>(iSum) * phi[iSum][iVec];
78 std::array<double,q> discIntegral;
79 for (
int iVec = 0; iVec < q; iVec++ ) {
80 discIntegral[iVec] = 0;
81 for (
int iSum = 0; iSum < q; iSum++ ) {
82 discIntegral[iVec] += descriptors::t<double,DESCRIPTOR>(iSum) * phi[iVec][iSum];
91 const double weights[q], std::array<std::array<double,q>,q>& cosTheta)
93 std::array<double,q> discIntegral;
94 for (
int iVec = 0; iVec < q; iVec++ ) {
95 discIntegral[iVec] = 0;
96 for (
int iSum = 0; iSum < q; iSum++ ) {
97 discIntegral[iVec] += weights[iSum] * cosTheta[iVec][iSum]* phi[iVec][iSum];
133 double solution[(DESCRIPTOR::q-1)*((DESCRIPTOR::q-1)+1)/2],
134 std::array<std::array<double,DESCRIPTOR::q-1>, DESCRIPTOR::q-1>& phi,
int const breakAfter = -1)
136 using namespace descriptors;
137 if ( !DESCRIPTOR::template provides<descriptors::tag::RTLBM>() ) {
138 std::cout <<
"Warning: Compute anisotropy matrix for wrong latice stencil!" << std::endl;
139 std::cout <<
"Weight for direction 0 is required to be 0." << std::endl;
144 clout <<
"Call AnistorpiyMatrix()" << std::endl;
145#ifdef FEATURE_OPENBLAS
146 clout <<
"Compute anisotropy matrix ..." << std::endl;
147 typedef DESCRIPTOR L;
148 int const q = L::q-1;
150 int const nn = (q+1)*q/2;
153 std::array<std::array<double,q>,q> angleProb;
157 for (
int iPop=0; iPop<q; iPop++) {
158 for (
int jPop=0; jPop<q; jPop++) {
161 dotProduct = c<L>(iPop+1)*c<L>(jPop+1);
162 normI =
util::sqrt( c<L>(iPop+1)*c<L>(iPop+1) );
163 normJ =
util::sqrt( c<L>(jPop+1)*c<L>(jPop+1) );
164 angleProb[iPop][jPop] = dotProduct / (normI*normJ);
168 for (
int i=0; i<q; i++) {
169 for (
int j=0; j<q; j++) {
174 for (
int i = 0; i < nn; i++ ) {
180 double U_array[nn*mm];
181 std::vector<double> anisoIterVector =
linespace( stepsize, 0, anisotropyFactor );
185 for ( ; index < anisoIterVector.size() && index != (std::size_t)(breakAfter); index++) {
187 for (
int m = 0; m < mm; m++) {
188 for (
int n = 0; n < nn; n++) {
195 for (
int m=0; m<q; m++) {
196 for (
int n=m; n<q; n++) {
197 U[m][iter] = t<double,L>(m+1)*phi[n][m];
198 U[n][iter] = t<double,L>(n+1)*phi[m][n];
199 U[m+q][iter] = t<double,L>(m+1)*phi[n][m]*angleProb[n][m];
200 U[n+q][iter] = t<double,L>(n+1)*phi[m][n]*angleProb[m][n];
208 for (
int n=0; n<q; n++) {
212 for (
int k=0; k<q; k++) {
213 sum1 += t<double,L>(k+1)*phi[k][n];
214 sum2 += t<double,L>(k+1)*phi[k][n]*angleProb[k][n];
216 solution[n] = 1 - sum1;
217 solution[q+n] = anisoIterVector[index] - sum2;
221 for (
int n = 0; n < nn; ++n) {
222 for (
int m = 0; m < mm; ++m ) {
223 U_array[n*mm +m] = U[m][n];
231 LAPACKE_dgels( LAPACK_COL_MAJOR,
'N', mm,
237 for (
int m=0; m<q; m++) {
238 for (
int n=m; n<q; n++) {
239 phi[n][m] = phi[n][m]*(1+solution[iter]);
240 phi[m][n] = phi[n][m];
248 clout <<
"Terminate after " << index <<
" iterations" << std::endl;
255 std::array<std::array<double,DESCRIPTOR::q>,DESCRIPTOR::q>& phi )
258 clout <<
"Compute anisotropy matrix ..." << std::endl;
259 typedef DESCRIPTOR L;
263 std::array<std::array<double,q>, q> cosTheta;
267 for (
int iPop=0; iPop<q; iPop++) {
268 for (
int jPop=0; jPop<q; jPop++) {
269 dotProduct = descriptors::c<L>(iPop) * descriptors::c<L>(jPop);
270 normI =
util::sqrt( util::normSqr<int,3>(descriptors::c<L>(iPop)) );
271 normJ =
util::sqrt( util::normSqr<int,3>(descriptors::c<L>(jPop)) );
272 cosTheta[iPop][jPop] = dotProduct /(normI*normJ);
273 if ( util::normSqr<int,3>(descriptors::c<L>(iPop)) == 0 || util::normSqr<int,3>(descriptors::c<L>(jPop)) == 0) {
274 cosTheta[iPop][jPop] = 0.0;
279 std::array<std::array<double,q>, q> phaseFunction;
280 for (
int m=0; m<q; m++) {
281 for (
int n=0; n<q; n++) {
286 for (
int m=0; m<q; m++) {
287 for (
int n=0; n<q; n++) {
289 for (
int i = 0; i < q; ++i) {
290 dotProduct += phaseFunction[m][i]*descriptors::t<double,L>(i);
292 phi[m][n] = phaseFunction[m][n] / dotProduct;
std::vector< double > linespace(double const stepsize, double const a, double const b)
Computes vector [a, a+stepsize, a+2*stepsize, ..., b].
std::array< double, q > testAnisotropyConservationColumn(const std::array< std::array< double, q >, q > &phi, const double weights[q], std::array< std::array< double, q >, q > &cosTheta)
void computeAnisotropyMatrix(double const stepsize, double const anisotropyFactor, double solution[(DESCRIPTOR::q-1) *((DESCRIPTOR::q-1)+1)/2], std::array< std::array< double, DESCRIPTOR::q-1 >, DESCRIPTOR::q-1 > &phi, int const breakAfter=-1)