KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
solid_mechanics_math_utilities.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosPoromechanicsApplication $
3 // Created by: $Author: JMCarbonell $
4 // Last modified by: $Co-Author: $
5 // Date: $Date: July 2013 $
6 // Revision: $Revision: 0.0 $
7 //
8 //
9 
10 #if !defined(KRATOS_SOLID_MECHANICS_MATH_UTILITIES_H_INCLUDED)
11 #define KRATOS_SOLID_MECHANICS_MATH_UTILITIES_H_INCLUDED
12 
13 
14 #ifdef FIND_MAX
15 #undef FIND_MAX
16 #endif
17 
18 #define FIND_MAX(a, b) ((a)>(b)?(a):(b))
19 
20 // System includes
21 #include <cmath>
22 
23 // External includes
24 
25 // Project includes
26 #include "utilities/math_utils.h"
27 #include "geometries/point.h"
28 
29 
30 namespace Kratos
31 {
32 template<class TDataType>
34 {
35 public:
40  typedef Matrix MatrixType;
41 
42  typedef Vector VectorType;
43 
44  typedef unsigned int IndexType;
45 
46  typedef unsigned int SizeType;
47 
49 
51 
53 
55 
57 
58 
75  static inline bool CardanoFormula(double a, double b, double c, double d, Vector& solution)
76  {
77  solution.resize(3,false);
79 
80  if(a==0)
81  {
82  std::cout<<"This is not a cubic equation: CardanoFormula"<<std::endl;
83 
84  return false;
85  }
86 
87  double p= (3.0*a*c-b*b)/(3.0*a*a);
88 
89  double q= 2.0*b*b*b/(27.0*a*a*a)-b*c/(3.0*a*a)+d/a;
90 
91  double discriminante= p*p*p/27.0+q*q/4.0;
92 
93  if(discriminante>0)
94  {
95  return false;
96  }
97 
98  if(discriminante==0)
99  {
100  if( a == 0 )
101  return false;
102 
103  solution(0)= pow(q/2.0, 1.0/3.0)-b/(3*a);
104  solution(1)= pow(q/2.0, 1.0/3.0)-b/(3*a);
105  solution(2)= pow(-4.0*q, 1.0/3.0)-b/(3*a);
106 
107  return true;
108  }
109 
110  if(p < 0) return false; //in this case the square roots below will be negative. This substitutes with better efficiency lines 107-110
111 
112  solution(0)=
113  -sqrt(-4.0/3.0*p)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p)))+ Globals::Pi / 3.0)
114  -b/(3*a);
115  solution(1)=
116  sqrt(-4.0/3.0*p)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p))))-b/(3*a)
117  ;
118  solution(2)=
119  -sqrt(-4.0/3.0*p)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p)))- Globals::Pi / 3.0)
120  -b/(3*a);
121 
122 // if(std::isnan<double>(solution(0)) || std::isnan<double>(solution(1))|| std::isnan<double>(solution(2)))
123 // {
124 // return false;
125 // }
126 
127  return true;
128  }
142  static inline Vector EigenValues(const Matrix& A, double tolerance, double zero)
143  {
144  int dim= A.size1();
145 
146  Matrix Convergence(2,dim);
147  noalias(Convergence) = ZeroMatrix(2,dim);
148 
149  double delta = 0.0;
150 
151  double abs = 0.0;
152 
153  Vector Result(dim);
154  noalias(Result) = ZeroVector(dim);
155 
156  Matrix HelpA(dim,dim);
157  noalias(HelpA) = ZeroMatrix(dim, dim);
158 
159  Matrix HelpQ(dim,dim);
160  noalias(HelpQ) = ZeroMatrix(dim, dim);
161 
162  Matrix HelpR(dim,dim);
163  noalias(HelpR) = ZeroMatrix(dim, dim);
164 
165  HelpA=A;
166 
167  bool is_converged=false;
168  unsigned int max_iters = 10;
169  unsigned int iter = 0;
170 
171  while(is_converged == false && iter<max_iters )
172  {
173  double shift= HelpA((dim-1),(dim-1));
174  //
175  for(int i=0; i<dim; i++)
176  {
177  HelpA(i,i) = HelpA(i,i)- shift;
178  }
179 
181 
182  HelpA= ZeroMatrix(dim, dim);
183 
184  for(int i=0; i<dim; i++)
185  {
186  HelpA(i,i) += shift;
187  for(int j=0; j< dim; j++)
188  {
189  for(int k=0; k< dim; k++)
190  {
191  HelpA(i,j) += HelpR(i,k)*HelpQ(k,j);
192  }
193  }
194  }
195 
196  delta= 0.0;
197 
198  abs = 0.0;
199 
200  for(int i=0; i<dim; i++)
201  {
202  Convergence(0,i)=Convergence(1,i);
203  Convergence(1,i)=HelpA(i,i);
204  delta+= (Convergence(1,i)-Convergence(0,i))*(Convergence(1,i)-Convergence(0,i));
205  abs+=(Convergence(1,i))*(Convergence(1,i));
206  }
207 
208  delta= sqrt(delta);
209 
210  abs=sqrt(abs);
211 
212  if(abs < zero)
213  abs=1.0;
214 
215  if(delta < zero || (delta/abs) < tolerance)
216  is_converged=true;
217 
218  iter++;
219  }
220 
221 
222  for(int i=0; i<dim; i++)
223  {
224  Result(i)= HelpA(i,i);
225 
226  if(fabs(Result(i)) <zero)
227  Result(i)=0.0;
228  }
229 
230  return Result;
231  }
232 
243  static inline Vector EigenValuesDirectMethod(const Matrix& A)
244  {
245  // Given a real symmetric 3x3 matrix A, compute the eigenvalues
246  int dim= A.size1();
247  Vector Result(dim);
248  noalias(Result) = ZeroVector(dim);
249 
250  const double p1 = A(0,1)*A(0,1) + A(0,2)*A(0,2) + A(1,2)*A(1,2);
251  if (p1 == 0)
252  {//A is diagonal.
253  Result[0] = A(0,0);
254  Result[1] = A(1,1);
255  Result[2] = A(2,2);
256  return Result;
257  }
258 
259  const double q = (A(0,0) + A(1,1) + A(2,2)) / 3.0;
260  const double p2 = (A(0,0) - q) * (A(0,0) - q) + (A(1,1) - q) * (A(1,1) - q) + (A(2,2) - q) * (A(2,2) - q) + 2.0 * p1;
261  const double p = sqrt(p2 / 6.0);
262 
263  Matrix B(3,3);
264  const double inv_p = 1.0/p;
265 
266  // B = (1 / p) * (A - q * I) where I is the identity matrix
267  B(0,0) = inv_p * (A(0,0) - q);
268  B(1,1) = inv_p * (A(1,1) - q);
269  B(2,2) = inv_p * (A(2,2) - q);
270  B(0,1) = inv_p * A(0,1);
271  B(1,0) = inv_p * A(1,0);
272  B(0,2) = inv_p * A(0,2);
273  B(2,0) = inv_p * A(2,0);
274  B(1,2) = inv_p * A(1,2);
275  B(2,1) = inv_p * A(2,1);
276 
277  //r = det(B) / 2
278  double r = 0.5 * ( B(0,0)*B(1,1)*B(2,2) + B(0,1)*B(1,2)*B(2,0) + B(1,0)*B(2,1)*B(0,2) - B(2,0)*B(1,1)*B(0,2) - B(1,0)*B(0,1)*B(2,2) - B(0,0)*B(2,1)*B(1,2) );
279 
280  // In exact arithmetic for a symmetric matrix -1 <= r <= 1
281  // but computation error can leave it slightly outside this range.
282  double phi = 0.0;
283  if (r <= -1) { phi = Globals::Pi / 3.0; }
284  else if (r >= 1) { phi = 0.0; }
285  else { phi = acos(r) / 3.0;}
286 
287  // the eigenvalues satisfy eig3 <= eig2 <= eig1
288  Result[0] = q + 2.0 * p * cos(phi);
289  Result[2] = q + 2.0 * p * cos(phi + (2.0 * Globals::Pi / 3.0));
290  Result[1] = 3.0 * q - Result[0] - Result[2]; //% since trace(A) = eig1 + eig2 + eig3
291 
292  return Result;
293  }
294 
295 
307  static inline void QRFactorization(const MatrixType& A, MatrixType& Q, MatrixType& R)
308  {
309 
310  //QR Factorization with Householder-Algo
311  int dim= A.size1();
312 
313  Vector y(dim);
314 
315  Vector w(dim);
316 
317  R.resize(dim,dim,false);
318 
320 
321  Q.resize(dim,dim,false);
322 
324 
325  Matrix Help(A.size1(),A.size2());
326  noalias(Help) = A;
327 
328  Matrix unity(dim,dim);
329  noalias(unity) = ZeroMatrix(dim,dim);
330 
331  for(int j=0; j<dim; j++)
332  unity(j,j)=1.0;
333 
334  std::vector<Matrix> HelpQ(dim-1);
335 
336  std::vector<Matrix> HelpR(dim-1);
337 
338  for(int i=0; i< dim-1; i++)
339  {
340  HelpQ[i].resize(dim,dim,false);
341  HelpR[i].resize(dim,dim,false);
342  noalias(HelpQ[i])= unity;
343  noalias(HelpR[i])= ZeroMatrix(dim,dim);
344  }
345 
346  for(int iteration=0; iteration< dim-1; iteration++)
347  {
348  //Vector y
349  for(int i=iteration; i<dim; i++)
350  y(i)= Help(i,iteration);
351 
352 
353  //Helpvalue l
354  double normy=0.0;
355 
356  for(int i=iteration; i<dim; i++)
357  normy += y(i)*y(i);
358 
359  normy= sqrt(normy);
360 
361  double l= sqrt((normy*(normy+fabs(y(iteration))))/2);
362 
363  double k=0.0;
364 
365  if(y[iteration] !=0)
366  k= - y(iteration)/fabs(y(iteration))*normy;
367  else
368  k= -normy;
369 
370  for(int i=iteration; i<dim; i++)
371  {
372  double e=0;
373 
374  if(i==iteration)
375  e=1;
376 
377  w(i)= 1/(2*l)*(y(i)-k*e);
378  }
379 
380  for(int i=iteration; i<dim; i++)
381  for(int j=iteration; j<dim; j++)
382  HelpQ[iteration](i,j)= unity(i,j)- 2*w(i)*w(j);
383 
384 
385  for(int i=iteration; i<dim; i++)
386  for(int j=iteration; j<dim; j++)
387  for(int k=iteration; k<dim; k++)
388  HelpR[iteration](i,j)+= HelpQ[iteration](i,k)*Help(k,j);
389 
390  Help= HelpR[iteration];
391 
392  }
393 
394  //Assembling R
395  for(int k=0; k<dim-1; k++)
396  {
397  for(int i=k; i<dim; i++)
398  for(int j=k; j<dim; j++)
399  R(i,j) =HelpR[k](i,j);
400 
401  }
402 
403 
404  for(int k=1; k<dim-1; k++)
405  {
406  for(int i=0; i<dim; i++)
407  for(int j=0; j<dim; j++)
408  for(int l=0; l<dim; l++)
409  Q(i,j)+= HelpQ[(k-1)](i,l)*HelpQ[k](l,j);
410  noalias(HelpQ[k])=Q;
411  }
412  if(dim-1==1)
413  noalias(Q)=HelpQ[0];
414 
415  }
416 
432  static inline void EigenVectors(const MatrixType& A,
433  MatrixType& vectors,
434  VectorType& lambda,
435  double zero_tolerance =1e-9,
436  int max_iterations = 10)
437  {
438  Matrix Help= A;
439 
440  for(int i=0; i<3; i++)
441  for(int j=0; j<3; j++)
442  Help(i,j)= Help(i,j);
443 
444 
445  vectors.resize(Help.size1(),Help.size2(),false);
446 
447  lambda.resize(Help.size1(),false);
448 
449  Matrix HelpDummy(Help.size1(),Help.size2());
450 
451  bool is_converged = false;
452 
453  Matrix unity(Help.size1(),Help.size2());
454  noalias(unity) = ZeroMatrix(Help.size1(),Help.size2());
455 
456  for(unsigned int i=0; i< Help.size1(); i++)
457  unity(i,i)= 1.0;
458 
459  Matrix V= unity;
460 
461  Matrix VDummy(Help.size1(),Help.size2());
462 
463  Matrix Rotation(Help.size1(),Help.size2());
464 
465 
466  for(int iterations=0; iterations<max_iterations; iterations++)
467  {
468 
469  is_converged= true;
470 
471  double a= 0.0;
472 
473  unsigned int index1= 0;
474 
475  unsigned int index2= 1;
476 
477  for(unsigned int i=0; i< Help.size1(); i++)
478  {
479  for(unsigned int j=(i+1); j< Help.size2(); j++)
480  {
481  if((fabs(Help(i,j)) > a ) && (fabs(Help(i,j)) > zero_tolerance))
482  {
483  a= fabs(Help(i,j));
484 
485  index1= i;
486  index2= j;
487 
488  is_converged= false;
489  }
490  }
491  }
492 
493 // KRATOS_WATCH( Help )
494 
495  if(is_converged)
496  break;
497 
498  //Calculation of Rotationangle
499 
500  double gamma= (Help(index2,index2)-Help(index1,index1))/(2*Help(index1,index2));
501 
502  double u=1.0;
503 
504  if(fabs(gamma) > zero_tolerance && fabs(gamma)< (1/zero_tolerance))
505  {
506  u= gamma/fabs(gamma)*1.0/(fabs(gamma)+sqrt(1.0+gamma*gamma));
507  }
508  else
509  {
510  if (fabs(gamma)>= (1.0/zero_tolerance))
511  u= 0.5/gamma;
512  }
513 
514  double c= 1.0/(sqrt(1.0+u*u));
515 
516  double s= c*u;
517 
518  double teta= s/(1.0+c);
519 
520  //Ratotion of the Matrix
521  HelpDummy= Help;
522 
523  HelpDummy(index2,index2)= Help(index2,index2)+u*Help(index1,index2);
524  HelpDummy(index1,index1)= Help(index1,index1)-u*Help(index1,index2);
525  HelpDummy(index1,index2)= 0.0;
526  HelpDummy(index2,index1)= 0.0;
527 
528  for(unsigned int i=0; i<Help.size1(); i++)
529  {
530  if((i!= index1) && (i!= index2))
531  {
532  HelpDummy(index2,i)=Help(index2,i)+s*(Help(index1,i)- teta*Help(index2,i));
533  HelpDummy(i,index2)=Help(index2,i)+s*(Help(index1,i)- teta*Help(index2,i));
534 
535  HelpDummy(index1,i)=Help(index1,i)-s*(Help(index2,i)+ teta*Help(index1,i));
536  HelpDummy(i,index1)=Help(index1,i)-s*(Help(index2,i)+ teta*Help(index1,i));
537  }
538  }
539 
540 
541  Help= HelpDummy;
542 
543  //Calculation of the eigenvectors V
544  Rotation =unity;
545  Rotation(index2,index1)=-s;
546  Rotation(index1,index2)=s;
547  Rotation(index1,index1)=c;
548  Rotation(index2,index2)=c;
549 
550 // Help.resize(A.size1(),A.size1(),false);
551 // noalias(Help)=ZeroMatrix(A.size1(),A.size1());
552 
553  noalias(VDummy) = ZeroMatrix(Help.size1(), Help.size2());
554 
555  for(unsigned int i=0; i< Help.size1(); i++)
556  {
557  for(unsigned int j=0; j< Help.size1(); j++)
558  {
559  for(unsigned int k=0; k< Help.size1(); k++)
560  {
561  VDummy(i,j) += V(i,k)*Rotation(k,j);
562  }
563  }
564  }
565  V= VDummy;
566 
567 // Matrix VTA(3,3);
568 // noalias(VTA) = ZeroMatrix(3,3);
569 // for(int i=0; i< Help.size1(); i++)
570 // {
571 // for(int j=0; j< Help.size1(); j++)
572 // {
573 // for(int k=0; k< Help.size1(); k++)
574 // {
575 // VTA(i,j) += V(k,i)*A(k,j);
576 // }
577 // }
578 // }
579 //
580 // for(int i=0; i< Help.size1(); i++)
581 // {
582 // for(int j=0; j< Help.size1(); j++)
583 // {
584 // for(int k=0; k< Help.size1(); k++)
585 // {
586 // Help(i,j) += VTA(i,k)*V(k,j);
587 // }
588 // }
589 // }
590 
591  }
592 
593  if(!(is_converged))
594  {
595  std::cout<<"########################################################"<<std::endl;
596  std::cout<<"Max_Iterations exceed in Jacobi-Seidel-Iteration (eigenvectors)"<<std::endl;
597  std::cout<<"########################################################"<<std::endl;
598  }
599 
600  for(unsigned int i=0; i< Help.size1(); i++)
601  {
602  for(unsigned int j=0; j< Help.size1(); j++)
603  {
604  vectors(i,j)= V(j,i);
605  }
606  }
607 
608  for(unsigned int i=0; i<Help.size1(); i++)
609  lambda(i)= Help(i,i);
610 
611  }
612 
613 
625  static inline void EigenVectors3x3(const Matrix& A, Matrix&V, Vector& d)
626  {
627 
628  if(A.size1()!=3 || A.size2()!=3)
629  std::cout<<" GIVEN MATRIX IS NOT 3x3 eigenvectors calculation "<<std::endl;
630 
631  Vector e(3);
632  noalias(e) = ZeroVector(3);
633  noalias(V) = ZeroMatrix(3,3);
634 
635  for (int i = 0; i < 3; i++) {
636  for (int j = 0; j < 3; j++) {
637  V(i,j) = A(i,j);
638  }
639  }
640 
641  // *******************//
642  //Symmetric Housholder reduction to tridiagonal form
643 
644  for (int j = 0; j < 3; j++) {
645  d[j] = V(2,j);
646  }
647 
648  // Householder reduction to tridiagonal form.
649 
650  for (int i = 2; i > 0; i--) {
651 
652  // Scale to avoid under/overflow.
653 
654  double scale = 0.0;
655  double h = 0.0;
656  for (int k = 0; k < i; k++) {
657  scale = scale + fabs(d[k]);
658  }
659  if (scale == 0.0) {
660  e[i] = d[i-1];
661  for (int j = 0; j < i; j++) {
662  d[j] = V(i-1,j);
663  V(i,j) = 0.0;
664  V(j,i) = 0.0;
665  }
666  } else {
667 
668  // Generate Householder vector.
669 
670  for (int k = 0; k < i; k++) {
671  d[k] /= scale;
672  h += d[k] * d[k];
673  }
674  double f = d[i-1];
675  double g = sqrt(h);
676  if (f > 0) {
677  g = -g;
678  }
679  e[i] = scale * g;
680  h = h - f * g;
681  d[i-1] = f - g;
682  for (int j = 0; j < i; j++) {
683  e[j] = 0.0;
684  }
685 
686  // Apply similarity transformation to remaining columns.
687 
688  for (int j = 0; j < i; j++) {
689  f = d[j];
690  V(j,i) = f;
691  g = e[j] + V(j,j) * f;
692  for (int k = j+1; k <= i-1; k++) {
693  g += V(k,j) * d[k];
694  e[k] += V(k,j) * f;
695  }
696  e[j] = g;
697  }
698  f = 0.0;
699  for (int j = 0; j < i; j++) {
700  e[j] /= h;
701  f += e[j] * d[j];
702  }
703  double hh = f / (h + h);
704  for (int j = 0; j < i; j++) {
705  e[j] -= hh * d[j];
706  }
707  for (int j = 0; j < i; j++) {
708  f = d[j];
709  g = e[j];
710  for (int k = j; k <= i-1; k++) {
711  V(k,j) -= (f * e[k] + g * d[k]);
712  }
713  d[j] = V(i-1,j);
714  V(i,j) = 0.0;
715  }
716  }
717  d[i] = h;
718  }
719 
720  // Accumulate transformations.
721 
722  for (int i = 0; i < 2; i++) {
723  V(2,i) = V(i,i);
724  V(i,i) = 1.0;
725  double h = d[i+1];
726  if (h != 0.0) {
727  for (int k = 0; k <= i; k++) {
728  d[k] = V(k,i+1) / h;
729  }
730  for (int j = 0; j <= i; j++) {
731  double g = 0.0;
732  for (int k = 0; k <= i; k++) {
733  g += V(k,i+1) * V(k,j);
734  }
735  for (int k = 0; k <= i; k++) {
736  V(k,j) -= g * d[k];
737  }
738  }
739  }
740  for (int k = 0; k <= i; k++) {
741  V(k,i+1) = 0.0;
742  }
743  }
744  for (int j = 0; j < 3; j++) {
745  d[j] = V(2,j);
746  V(2,j) = 0.0;
747  }
748  V(2,2) = 1.0;
749  e[0] = 0.0;
750 
751  // *******************//
752 
753  // Symmetric tridiagonal QL algorithm.
754 
755  for (int i = 1; i < 3; i++) {
756  e[i-1] = e[i];
757  }
758  e[2] = 0.0;
759 
760  double f = 0.0;
761  double tst1 = 0.0;
762  double eps = std::pow(2.0,-52.0);
763  for (int l = 0; l < 3; l++) {
764 
765  // Find small subdiagonal element
766 
767  tst1 = FIND_MAX(tst1,fabs(d[l]) + fabs(e[l]));
768  int m = l;
769  while (m < 3) {
770  if (fabs(e[m]) <= eps*tst1) {
771  break;
772  }
773  m++;
774  }
775 
776  // If m == l, d[l] is an eigenvalue,
777  // otherwise, iterate.
778 
779  if (m > l) {
780  int iter = 0;
781  do {
782  iter = iter + 1; // (Could check iteration count here.)
783 
784  // Compute implicit shift
785 
786  double g = d[l];
787  double p = (d[l+1] - g) / (2.0 * e[l]);
788  double r = hypot2(p,1.0);
789  if (p < 0) {
790  r = -r;
791  }
792  d[l] = e[l] / (p + r);
793  d[l+1] = e[l] * (p + r);
794  double dl1 = d[l+1];
795  double h = g - d[l];
796  for (int i = l+2; i < 3; i++) {
797  d[i] -= h;
798  }
799  f = f + h;
800 
801  // Implicit QL transformation.
802 
803  p = d[m];
804  double c = 1.0;
805  double c2 = c;
806  double c3 = c;
807  double el1 = e[l+1];
808  double s = 0.0;
809  double s2 = 0.0;
810  for (int i = m-1; i >= l; i--) {
811  c3 = c2;
812  c2 = c;
813  s2 = s;
814  g = c * e[i];
815  h = c * p;
816  r = hypot2(p,e[i]);
817  e[i+1] = s * r;
818  s = e[i] / r;
819  c = p / r;
820  p = c * d[i] - s * g;
821  d[i+1] = h + s * (c * g + s * d[i]);
822 
823  // Accumulate transformation.
824 
825  for (int k = 0; k < 3; k++) {
826  h = V(k,i+1);
827  V(k,i+1) = s * V(k,i) + c * h;
828  V(k,i) = c * V(k,i) - s * h;
829  }
830  }
831  p = -s * s2 * c3 * el1 * e[l] / dl1;
832  e[l] = s * p;
833  d[l] = c * p;
834 
835  // Check for convergence.
836 
837  } while (fabs(e[l]) > eps*tst1);
838  }
839  d[l] = d[l] + f;
840  e[l] = 0.0;
841  }
842 
843  // Sort eigenvalues and corresponding vectors.
844 
845  for (int i = 0; i < 2; i++) {
846  int k = i;
847  double p = d[i];
848  for (int j = i+1; j < 3; j++) {
849  if (d[j] < p) {
850  k = j;
851  p = d[j];
852  }
853  }
854  if (k != i) {
855  d[k] = d[i];
856  d[i] = p;
857  for (int j = 0; j < 3; j++) {
858  p = V(j,i);
859  V(j,i) = V(j,k);
860  V(j,k) = p;
861  }
862  }
863  }
864 
865  V = trans(V);
866 
867  // *******************//
868 
869 
870  }
871 
872 
873  static double hypot2(double x, double y) {
874  return std::sqrt(x*x+y*y);
875  }
876 
890  static inline void EigenVectors( MatrixType& A,
891  MatrixType& V,
892  TDataType& error_tolerance,
893  TDataType zero_tolerance)
894  {
895  //initial error
896  TDataType error = 1.0;
897  int n = A.size2();
898  //setting V to identity matrix
899  if( V.size1() != n )
900  V.resize(n,n,false);
902  //calculation loop (as long as there is no convergence)
903  //WARNING: iteration never exceeds
904  while( error > error_tolerance )
905  {
906  for( int i=0; i<n; i++ )
907  {
908  for( int j=i+1; j<n; j++ )
909  {
910  double theta = 0.0;
911  if( std::abs( A(i,j) ) >= zero_tolerance )
912  {
913  if( std::abs( A(i,i)-A(j,j) ) > 0.0 )
914  {
915  theta = 0.5*atan(2*A(i,j)/(A(i,i)-A(j,j)));
916  }
917  else theta = 0.25 * Globals::Pi;
918  }
919  MatrixType T(n,n);
920  noalias(T) = IdentityMatrix( n );
921 
922  T(i,i) = cos(theta);
923  T(i,j) = -sin(theta);
924  T(j,i) = -T(i,j);
925  T(j,j) = T(i,i);
926 
927  A = Mult( A, T );
928  MatrixType TT = Transpose(T);
929  A = Mult( TT, A );
930  V = Mult( V, T );
931  }
932  }
933  double sTot = 0.0;
934  double sDiag = 0.0;
935  for( unsigned int i=0; i<A.size1(); i++ )
936  {
937  for( unsigned int j=0; j<A.size2(); j++ )
938  {
939  sTot += std::abs(A(i,j));
940  }
941  sDiag+= std::abs(A(i,i));
942  }
943  error=(sTot-sDiag)/sDiag;
944  }
945  //sorting eigenvalues
946  int maxIndex = 0;
947  TDataType maxEv = A(0,0);
948  for( unsigned int i=0; i<A.size1(); i++ )
949  {
950  for( unsigned int j=i; j<A.size1(); j++ )
951  {
952  //searching current maximum
953  if( A(j,j) > maxEv )
954  {
955  maxIndex = j;
956  maxEv = A(j,j);
957  }
958  //swapping eigenvalue matrix
959  TDataType dummy = A(i,i);
960  A(i,i) = A(maxIndex,maxIndex);
961  A(maxIndex,maxIndex) = dummy;
962  //swapping eigenvector matrix
963  for( unsigned int k=0; k<A.size2(); k++ )
964  {
965  dummy = V(k,i);
966  V(k,i) = V(k,maxIndex);
967  V(k,maxIndex) = dummy;
968  }
969  }
970 
971  }
972  }//EigenVectors
973 
979  static inline MatrixType IdentityMatrix( SizeType size )
980  {
981  MatrixType A(size,size);
982  noalias(A) = ZeroMatrix( size, size );
983  for( unsigned int i=0; i<size ; i++ )
984  {
985  A(i,i) = 1.0;
986  }
987  return A;
988  }//IdentityMatrix
989 
996  static inline void Add( MatrixType& A, MatrixType& B )
997  {
998  for( unsigned int i=0; i<A.size1(); i++ )
999  {
1000  for( unsigned int j=0; j<A.size2(); j++ )
1001  {
1002  A(i,j) += B(i,j);
1003  }
1004  }
1005  }
1006 
1013  static inline MatrixType Mult( MatrixType& A,
1014  MatrixType& B)
1015  {
1016  MatrixType C(A.size1(),B.size2());
1017  for( unsigned int i=0; i<A.size1(); i++ )
1018  {
1019  for( unsigned int j=0; j<B.size2(); j++ )
1020  {
1021  for( unsigned int k=0; k<A.size2(); k++ )
1022  {
1023  C(i,j) += B(k,j)*A(i,k);
1024  }
1025  }
1026  }
1027  return C;
1028  }//Mult
1029 
1033  static inline void Mult( MatrixType& M, TDataType a )
1034  {
1035  for( unsigned int i=0; i<M.size1(); i++ )
1036  {
1037  for( unsigned int j=0; j<M.size2(); j++ )
1038  {
1039  M(i,j) = M(i,j)*a;
1040  }
1041  }
1042  }
1043 
1047  static inline void Mult( VectorType& v, TDataType a )
1048  {
1049  for( unsigned int i=0; i<v.size(); i++ )
1050  {
1051  v(i) = v(i)*a;
1052  }
1053  }//Mult
1054 
1060  static inline MatrixType Transpose( MatrixType& A )
1061  {
1062  MatrixType AT(A.size2(),A.size1());
1063  noalias(AT) = ZeroMatrix(A.size2(),A.size1());
1064  for( unsigned int i=0; i<A.size1(); i++ )
1065  {
1066  for( unsigned int j=0; j<A.size2(); j++ )
1067  {
1068  AT(j,i) = A(i,j);
1069  }
1070  }
1071  return AT;
1072  }//Transpose
1073 
1077  static inline void Normalize( VectorType& v )
1078  {
1079  Mult( v, 1.0/(MathUtilsType::Norm( v )) );
1080  }
1081 
1092  static inline MatrixType StrainVectorToTensor( const VectorType& Strains )
1093  {
1094  KRATOS_TRY
1095 
1096  Matrix StrainTensor;
1097  //KRATOS_WATCH( Strains )
1098  if (Strains.size()==3)
1099  {
1100  StrainTensor.resize(2,2, false);
1101  //KRATOS_WATCH( StrainTensor )
1102  StrainTensor(0,0) = Strains[0];
1103  StrainTensor(0,1) = 0.5*Strains[2];
1104  StrainTensor(1,0) = 0.5*Strains[2];
1105  StrainTensor(1,1) = Strains[1];
1106  }
1107  else if (Strains.size()==6)
1108  {
1109  StrainTensor.resize(3,3, false);
1110  StrainTensor(0,0) = Strains[0];
1111  StrainTensor(0,1) = 0.5*Strains[3];
1112  StrainTensor(0,2) = 0.5*Strains[5];
1113  StrainTensor(1,0) = 0.5*Strains[3];
1114  StrainTensor(1,1) = Strains[1];
1115  StrainTensor(1,2) = 0.5*Strains[4];
1116  StrainTensor(2,0) = 0.5*Strains[5];
1117  StrainTensor(2,1) = 0.5*Strains[4];
1118  StrainTensor(2,2) = Strains[2];
1119  }
1120 
1121  //KRATOS_WATCH( StrainTensor )
1122  return StrainTensor;
1123  KRATOS_CATCH( "" )
1124  }
1125 
1126  static inline Vector TensorToStrainVector( const Matrix& Tensor )
1127  {
1128  KRATOS_TRY
1129  Vector StrainVector;
1130 
1131 
1132  if (Tensor.size1()==2)
1133  {
1134  StrainVector.resize(3);
1135  noalias(StrainVector) = ZeroVector(3);
1136  StrainVector(0) = Tensor(0,0);
1137  StrainVector(1) = Tensor(1,1);
1138  StrainVector(2) = 2.00*Tensor(0,1);
1139  }
1140  else if (Tensor.size1()==3)
1141  {
1142  StrainVector.resize(6);
1143  noalias(StrainVector) = ZeroVector(6);
1144  StrainVector(0) = Tensor(0,0);
1145  StrainVector(1) = Tensor(1,1);
1146  StrainVector(2) = Tensor(2,2);
1147  StrainVector(3) = 2.00*Tensor(0,1);
1148  StrainVector(4) = 2.00*Tensor(1,2);
1149  StrainVector(5) = 2.00*Tensor(0,2);
1150  }
1151 
1152 
1153  return StrainVector;
1154  KRATOS_CATCH( "" )
1155  }
1156 
1157 
1163  static int InvertMatrix( const MatrixType& input, MatrixType& inverse )
1164  {
1165  typedef permutation_matrix<std::size_t> pmatrix;
1166  Matrix A(input);
1167  pmatrix pm(A.size1());
1168  int singular = lu_factorize(A,pm);
1169  inverse.assign( IdentityMatrix(A.size1()));
1170  lu_substitute(A, pm, inverse);
1171  return singular;
1172  }
1173 
1179  static double NormTensor(Matrix& Tensor)
1180  {
1181  double result=0.0;
1182  for(unsigned int i=0; i< Tensor.size1(); i++)
1183  for(unsigned int j=0; j< Tensor.size2(); j++)
1184  result+= Tensor(i,j)*Tensor(i,j);
1185 
1186  result= sqrt(result);
1187 
1188  return result;
1189  }
1190 
1196  static inline void VectorToTensor(const Vector& Stress,Matrix& Tensor)
1197  {
1198  if(Stress.size()==6)
1199  {
1200  Tensor.resize(3,3,false);
1201  Tensor(0,0)= Stress(0);
1202  Tensor(0,1)= Stress(3);
1203  Tensor(0,2)= Stress(5);
1204  Tensor(1,0)= Stress(3);
1205  Tensor(1,1)= Stress(1);
1206  Tensor(1,2)= Stress(4);
1207  Tensor(2,0)= Stress(5);
1208  Tensor(2,1)= Stress(4);
1209  Tensor(2,2)= Stress(2);
1210  }
1211  if(Stress.size()==3)
1212  {
1213  Tensor.resize(2,2,false);
1214  Tensor(0,0)= Stress(0);
1215  Tensor(0,1)= Stress(2);
1216  Tensor(1,0)= Stress(2);
1217  Tensor(1,1)= Stress(1);
1218  }
1219 
1220  }
1221 
1227  static void TensorToVector( const Matrix& Tensor, Vector& Vector)
1228  {
1229  //if(Vector.size()!= 6)
1230  unsigned int dim = Tensor.size1();
1231  if (dim==3)
1232  {
1233  Vector.resize(6,false);
1234  Vector(0)= Tensor(0,0);
1235  Vector(1)= Tensor(1,1);
1236  Vector(2)= Tensor(2,2);
1237  Vector(3)= Tensor(0,1);
1238  Vector(4)= Tensor(1,2);
1239  Vector(5)= Tensor(2,0);
1240  }
1241  else if(dim==2)
1242  {
1243  Vector.resize(3,false);
1244  Vector(0)= Tensor(0,0);
1245  Vector(1)= Tensor(1,1);
1246  Vector(2)= Tensor(0,1);
1247  }
1248 
1249  }
1250 
1251  static inline void TensorToMatrix(Fourth_Order_Tensor& Tensor,Matrix& Matrix)
1252  {
1253 
1254 
1255  // Simetrias seguras
1256  // Cijkl = Cjilk;
1257  // Cijkl = Cklji;
1258  if (Tensor[0].size()== 3)
1259  {
1260  // Tensor de cuarto orden cuyos componentes correspondes a una matriz de 3x3
1261  if(Matrix.size1()!=6 || Matrix.size2()!=6)
1262  Matrix.resize(6,6,false);
1263  Matrix(0,0) = Tensor[0][0](0,0);
1264  Matrix(0,1) = Tensor[0][0](1,1);
1265  Matrix(0,2) = Tensor[0][0](2,2);
1266  Matrix(0,3) = Tensor[0][0](0,1);
1267  Matrix(0,4) = Tensor[0][0](0,2);
1268  Matrix(0,5) = Tensor[0][0](1,2);
1269 
1270  Matrix(1,0) = Tensor[1][1](0,0);
1271  Matrix(1,1) = Tensor[1][1](1,1);
1272  Matrix(1,2) = Tensor[1][1](2,2);
1273  Matrix(1,3) = Tensor[1][1](0,1);
1274  Matrix(1,4) = Tensor[1][1](0,2);
1275  Matrix(1,5) = Tensor[1][1](1,2);
1276 
1277  Matrix(2,0) = Tensor[2][2](0,0);
1278  Matrix(2,1) = Tensor[2][2](1,1);
1279  Matrix(2,2) = Tensor[2][2](2,2);
1280  Matrix(2,3) = Tensor[2][2](0,1);
1281  Matrix(2,4) = Tensor[2][2](0,2);
1282  Matrix(2,5) = Tensor[2][2](1,2);
1283 
1284  Matrix(3,0) = Tensor[0][1](0,0);
1285  Matrix(3,1) = Tensor[0][1](1,1);
1286  Matrix(3,2) = Tensor[0][1](2,2);
1287  Matrix(3,3) = Tensor[0][1](0,1);
1288  Matrix(3,4) = Tensor[0][1](0,2);
1289  Matrix(3,5) = Tensor[0][1](1,2);
1290 
1291  Matrix(4,0) = Tensor[0][2](0,0);
1292  Matrix(4,1) = Tensor[0][2](1,1);
1293  Matrix(4,2) = Tensor[0][2](2,2);
1294  Matrix(4,3) = Tensor[0][2](0,1);
1295  Matrix(4,4) = Tensor[0][2](0,2);
1296  Matrix(4,5) = Tensor[0][2](1,2);
1297 
1298  Matrix(5,0) = Tensor[1][2](0,0);
1299  Matrix(5,1) = Tensor[1][2](1,1);
1300  Matrix(5,2) = Tensor[1][2](2,2);
1301  Matrix(5,3) = Tensor[1][2](0,1);
1302  Matrix(5,4) = Tensor[1][2](0,2);
1303  Matrix(5,5) = Tensor[1][2](1,2);
1304  }
1305  else
1306  {
1307  // Tensor de cuarto orden cuyos componentes correspondes a una matriz de 2x2
1308  if(Matrix.size1()!=3 || Matrix.size2()!=3)
1309  Matrix.resize(3,3,false);
1310  Matrix(0,0) = Tensor[0][0](0,0);
1311  Matrix(0,1) = Tensor[0][0](1,1);
1312  Matrix(0,2) = Tensor[0][0](0,1);
1313  Matrix(1,0) = Tensor[1][1](0,0);
1314  Matrix(1,1) = Tensor[1][1](1,1);
1315  Matrix(1,2) = Tensor[1][1](0,1);
1316  Matrix(2,0) = Tensor[0][1](0,0);
1317  Matrix(2,1) = Tensor[0][1](1,1);
1318  Matrix(2,2) = Tensor[0][1](0,1);
1319 
1320  }
1321 
1322  }
1323 
1329  static void MatrixToTensor(MatrixType& A,std::vector<std::vector<Matrix> >& Tensor)
1330  {
1331  int help1 = 0;
1332  int help2 = 0;
1333  double coeff = 1.0;
1334 
1335  Tensor.resize(3);
1336 
1337  for(unsigned int i=0; i<3; i++)
1338  {
1339  Tensor[i].resize(3);
1340  for(unsigned int j=0; j<3; j++)
1341  {
1342  Tensor[i][j].resize(3,3,false);
1343  noalias(Tensor[i][j])= ZeroMatrix(3,3);
1344  for(unsigned int k=0; k<3; k++)
1345  for(unsigned int l=0; l<3; l++)
1346  {
1347  if(i==j) help1= i;
1348  else
1349  {
1350  if((i==0 && j==1) || (i==1 && j==0)) help1= 3;
1351  if((i==1 && j==2) || (i==2 && j==1)) help1= 4;
1352  if((i==2 && j==0) || (i==0 && j==2)) help1= 5;
1353  }
1354  if(k==l)
1355  {
1356  help2= k;
1357  coeff=1.0;
1358  }
1359  else
1360  {
1361  coeff=0.5;
1362  if((k==0 && l==1) || (k==1 && l==0)) help2= 3;
1363  if((k==1 && l==2) || (k==2 && l==1)) help2= 4;
1364  if((k==2 && l==0) || (k==0 && l==2)) help2= 5;
1365  }
1366 
1367  Tensor[i][j](k,l)= A(help1,help2)*coeff;
1368  }
1369  }
1370  }
1371 
1372 
1373  }
1380  {
1381  int help1 = 0;
1382  int help2 = 0;
1383  double coeff = 1.0;
1384  std::fill(Tensor.begin(), Tensor.end(), 0.0);
1385  for(unsigned int i=0; i<3; i++)
1386  {
1387  for(unsigned int j=0; j<3; j++)
1388  {
1389  for(unsigned int k=0; k<3; k++)
1390  for(unsigned int l=0; l<3; l++)
1391  {
1392  if(i==j) help1= i;
1393  else
1394  {
1395  if((i==0 && j==1) || (i==1 && j==0)) help1= 3;
1396  if((i==1 && j==2) || (i==2 && j==1)) help1= 4;
1397  if((i==2 && j==0) || (i==0 && j==2)) help1= 5;
1398  }
1399  if(k==l)
1400  {
1401  help2= k;
1402  coeff=1.0;
1403  }
1404  else
1405  {
1406  coeff=0.5;
1407  if((k==0 && l==1) || (k==1 && l==0)) help2= 3;
1408  if((k==1 && l==2) || (k==2 && l==1)) help2= 4;
1409  if((k==2 && l==0) || (k==0 && l==2)) help2= 5;
1410  }
1411 
1412  Tensor[i*27+j*9+k*3+l]= A(help1,help2)*coeff;
1413  }
1414  }
1415  }
1416 
1417 
1418  }
1424  static void TensorToMatrix(std::vector<std::vector<Matrix> >& Tensor,Matrix& Matrix)
1425  {
1426  int help1 = 0;
1427  int help2 = 0;
1428  int help3 = 0;
1429  int help4 = 0;
1430  double coeff = 1.0;
1431 
1432  if(Matrix.size1()!=6 || Matrix.size2()!=6)
1433  Matrix.resize(6,6,false);
1434 
1435  for(unsigned int i=0; i<6; i++)
1436  for(unsigned int j=0; j<6; j++)
1437  {
1438  if(i<3)
1439  {
1440  help1= i;
1441  help2= i;
1442  }
1443  else
1444  {
1445  if(i==3)
1446  {
1447  help1= 0;
1448  help2= 1;
1449  }
1450  if(i==4)
1451  {
1452  help1= 1;
1453  help2= 2;
1454  }
1455  if(i==5)
1456  {
1457  help1= 2;
1458  help2= 0;
1459  }
1460  }
1461 
1462  if(j<3)
1463  {
1464  help3= j;
1465  help4= j;
1466  coeff= 1.0;
1467  }
1468  else
1469  {
1470  if(j==3)
1471  {
1472  help3= 0;
1473  help4= 1;
1474  }
1475  if(j==4)
1476  {
1477  help3= 1;
1478  help4= 2;
1479  }
1480  if(j==5)
1481  {
1482  help3= 2;
1483  help4= 0;
1484  }
1485  coeff= 2.0;
1486  }
1487 
1488  Matrix(i,j)= Tensor[help1][help2](help3,help4)*coeff;
1489  }
1490 
1491  }
1492 
1498  static void TensorToMatrix( const array_1d<double, 81>& Tensor, Matrix& Matrix )
1499  {
1500  if(Matrix.size1()!=6 || Matrix.size2()!=6)
1501  Matrix.resize(6,6,false);
1502 
1503  Matrix(0,0) = Tensor[0];
1504  Matrix(0,1) = Tensor[4];
1505  Matrix(0,2) = Tensor[8];
1506  Matrix(0,3) = 2.0*Tensor[1];
1507  Matrix(0,4) = 2.0*Tensor[5];
1508  Matrix(0,5) = 2.0*Tensor[6];
1509 
1510  Matrix(1,0) = Tensor[36];
1511  Matrix(1,1) = Tensor[40];
1512  Matrix(1,2) = Tensor[44];
1513  Matrix(1,3) = 2.0*Tensor[37];
1514  Matrix(1,4) = 0.0*Tensor[41];
1515  Matrix(1,5) = 0.0*Tensor[42];
1516 
1517  Matrix(2,0) = Tensor[72];
1518  Matrix(2,1) = Tensor[76];
1519  Matrix(2,2) = Tensor[80];
1520  Matrix(2,3) = 2.0*Tensor[73];
1521  Matrix(2,4) = 2.0*Tensor[77];
1522  Matrix(2,5) = 2.0*Tensor[78];
1523 
1524  Matrix(3,0) = Tensor[9];
1525  Matrix(3,1) = Tensor[13];
1526  Matrix(3,2) = Tensor[18];
1527  Matrix(3,3) = 2.0*Tensor[10];
1528  Matrix(3,4) = 2.0*Tensor[14];
1529  Matrix(3,5) = 2.0*Tensor[15];
1530 
1531  Matrix(4,0) = Tensor[45];
1532  Matrix(4,1) = Tensor[49];
1533  Matrix(4,2) = Tensor[53];
1534  Matrix(4,3) = 2.0*Tensor[46];
1535  Matrix(4,4) = 0.0*Tensor[50];
1536  Matrix(4,5) = 0.0*Tensor[51];
1537 
1538  Matrix(5,0) = Tensor[54];
1539  Matrix(5,1) = Tensor[58];
1540  Matrix(5,2) = Tensor[62];
1541  Matrix(5,3) = 2.0*Tensor[55];
1542  Matrix(5,4) = 2.0*Tensor[59];
1543  Matrix(5,5) = 2.0*Tensor[60];
1544 
1545  }
1550  static void DeviatoricUnity(std::vector<std::vector<Matrix> >& Unity)
1551  {
1552  Unity.resize(3);
1553 
1554  Matrix kronecker(3,3);
1555  noalias(kronecker)=ZeroMatrix(3,3);
1556  for(unsigned int i=0; i<3; i++)
1557  {
1558  kronecker(i,i)=1;
1559  }
1560 
1561  for(unsigned int i=0; i<3; i++)
1562  {
1563  Unity[i].resize(3);
1564  for(unsigned int j=0; j<3; j++)
1565  {
1566  Unity[i][j].resize(3,3,false);
1567  noalias(Unity[i][j])= ZeroMatrix(3,3);
1568 
1569  for(unsigned int k=0; k<3; k++)
1570  {
1571  for(unsigned int l=0; l<3; l++)
1572  {
1573  Unity[i][j](k,l)=kronecker(i,k)*kronecker(j,l)
1574  -1.0/3.0*kronecker(i,j)*kronecker(k,l);
1575  }
1576  }
1577  }
1578  }
1579  }
1580 
1586  {
1587  Matrix kronecker(3,3);
1588  noalias(kronecker)=ZeroMatrix(3,3);
1589  for(unsigned int i=0; i<3; i++)
1590  {
1591  kronecker(i,i)=1;
1592  }
1593 
1594  for(unsigned int i=0; i<3; i++)
1595  for(unsigned int j=0; j<3; j++)
1596  for(unsigned int k=0; k<3; k++)
1597  for(unsigned int l=0; l<3; l++)
1598  Unity[27*i+9*j+3*k+l]=kronecker(i,k)*kronecker(j,l)
1599  -1.0/3.0*kronecker(i,j)*kronecker(k,l);
1600  }
1601 
1612  static bool Clipping(std::vector<Point* >& clipping_points,std::vector<Point* >& subjected_points, std::vector<Point* >& result_points, double tolerance)
1613  {
1614  result_points= subjected_points;
1615  Vector actual_edge(3);
1616  Vector actual_normal(3);
1617  std::vector<Point* > temp_results;
1618  bool is_visible= false;
1619  for(unsigned int clipp_edge=0; clipp_edge<clipping_points.size(); clipp_edge++)
1620  {
1621  temp_results.clear();
1622  unsigned int index_clipp_2=0;
1623  if(clipp_edge< (clipping_points.size()-1))
1624  index_clipp_2= clipp_edge+1;
1625  //define clipping edge vector
1626  noalias(actual_edge)= *(clipping_points[clipp_edge])-*(clipping_points[index_clipp_2]);
1627  noalias(actual_edge)= actual_edge/sqrt(inner_prod(actual_edge,actual_edge));
1628 
1629  //define normal on clipping-edge vector towards visible side
1630  if(clipp_edge< (clipping_points.size()-2))
1631  actual_normal=*(clipping_points[clipp_edge+2])-(*(clipping_points[clipp_edge])+actual_edge*inner_prod((*(clipping_points[clipp_edge+2])-*(clipping_points[clipp_edge])),actual_edge));
1632  else if(clipp_edge< (clipping_points.size()-1))
1633  actual_normal=*(clipping_points[0])-(*(clipping_points[clipp_edge])+actual_edge*inner_prod((*(clipping_points[0])-*(clipping_points[clipp_edge])),actual_edge));
1634  else
1635  actual_normal=*(clipping_points[1])-(*(clipping_points[clipp_edge])+actual_edge*inner_prod((*(clipping_points[1])-*(clipping_points[clipp_edge])),actual_edge));
1636 
1637  noalias(actual_normal)=actual_normal/(sqrt(inner_prod(actual_normal,actual_normal)));
1638 
1639  //test if the first point is visible or unvisible
1640  if( inner_prod((*(result_points[0])-*(clipping_points[clipp_edge])), actual_normal)> tolerance)
1641  {
1642  is_visible= true;
1643  }
1644  else
1645  {
1646  is_visible= false;
1647  }
1648 
1649  for(unsigned int subj_edge=0; subj_edge< result_points.size(); subj_edge++)
1650  {
1651  unsigned int index_subj_2=0;
1652 
1653  if(subj_edge< (result_points.size()-1))
1654  index_subj_2= subj_edge+1;
1655 
1656  //Test whether the points of the actual subj_edge lay on clipp_edge
1657  if(fabs(inner_prod((*(result_points[subj_edge])-(*(clipping_points[clipp_edge]))), actual_normal))<= tolerance)
1658  {
1659  temp_results.push_back(result_points[subj_edge]);
1660  if( inner_prod((*(result_points[index_subj_2])-*(clipping_points[clipp_edge])), actual_normal)> tolerance)
1661  is_visible= true;
1662  else
1663  is_visible= false;
1664 
1665  continue;
1666  }
1667  //Calculate minimal distance between the two points
1668  Vector b(2);
1669  b(0)= -inner_prod((*(result_points[index_subj_2])-*(result_points[subj_edge])),(*(result_points[subj_edge])-*(clipping_points[clipp_edge])));
1670  b(1)= inner_prod((*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge])),(*(result_points[subj_edge])-*(clipping_points[clipp_edge])));
1671  Matrix A(2,2);
1672  A(0,0)=inner_prod((*(result_points[index_subj_2])-*(result_points[subj_edge])),(*(result_points[index_subj_2])-*(result_points[subj_edge])));
1673  A(0,1)=-inner_prod((*(result_points[index_subj_2])-*(result_points[subj_edge])),(*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge])));
1674  A(1,0)= A(0,1);
1675  A(1,1)=inner_prod(*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge]),*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge]));
1676  Vector coeff(2);
1677  coeff(0)=1.0/A(0,0)*(b(0)-A(0,1)/(A(1,1)-A(0,1)*A(1,0)/A(0,0))*(b(1)-b(0)*A(1,0)/A(0,0)));
1678  coeff(1)=1.0/(A(1,1)-A(0,1)*A(1,0)/A(0,0))*(b(1)-b(0)*A(1,0)/A(0,0));
1679 
1680 
1681  //TEST on distance to endpoints of the line
1682  Vector dist_vec(3);
1683  noalias(dist_vec)= *(result_points[subj_edge])+coeff(0)*(*(result_points[index_subj_2])-*(result_points[subj_edge]))-(*(clipping_points[clipp_edge])+coeff(1)*(*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge])));
1684 
1685  if( coeff(0) > tolerance && coeff(0) < (1-tolerance)&& (sqrt(inner_prod(dist_vec,dist_vec))< tolerance))
1686  {
1687  if(is_visible)
1688  temp_results.push_back(result_points[subj_edge]);
1689 
1690  temp_results.push_back(new Point((*(clipping_points[clipp_edge])+coeff(1)*(*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge])))));
1691 
1692  is_visible= !is_visible;
1693 
1694  continue;
1695  }
1696  if(is_visible)
1697  temp_results.push_back(result_points[subj_edge]);
1698  }
1699  result_points=temp_results;
1700  }
1701  if(result_points.size()==0)
1702  return false;
1703  else
1704  return true;
1705  }
1706 
1707 
1708  inline static VectorType range(const unsigned int start, const unsigned int end)
1709  {
1710  int size = end-start;
1711  if(size<0)
1712  KRATOS_ERROR<<" range out of bounds start:"<<start<<" end:"<<end<<std::endl;
1713 
1714  VectorType b(size);
1715 
1716  for(int i = 0; i<size; ++i)
1717  b[i] = start + i;
1718 
1719  return b;
1720  }
1721 
1722 
1723  inline static MatrixType project(const MatrixType& rMatrix, VectorType irange, VectorType jrange)
1724  {
1725 
1726  MatrixType A(irange.size(), jrange.size());
1727 
1728  for(unsigned int i = 0; i<irange.size(); ++i)
1729  for(unsigned int j = 0; j<jrange.size(); ++j)
1730  A(i,j) = rMatrix(irange[i],jrange[j]);
1731 
1732  return A;
1733  }
1734 
1735 
1736  inline static MatrixType& project(MatrixType& rMatrixA, VectorType irange, VectorType jrange, const MatrixType& rMatrixB)
1737  {
1738 
1739  for(unsigned int i = 0; i<irange.size(); ++i)
1740  for(unsigned int j = 0; j<jrange.size(); ++j)
1741  rMatrixA(irange[i],jrange[j]) = rMatrixB(i,j);
1742 
1743  return rMatrixA;
1744  }
1745 
1746 
1747  inline static VectorType project(const VectorType& rVector, VectorType irange)
1748  {
1749 
1750  VectorType b(irange.size());
1751 
1752  for(unsigned int i = 0; i<irange.size(); ++i)
1753  b[i] = rVector[irange[i]];
1754 
1755  return b;
1756  }
1757 
1758 
1759 
1760 private:
1761 };// class SolidMechanicsMathUtilities
1762 }
1763 #endif /* KRATOS_SOLID_MECHANICS_MATH_UTILITIESS_H_INCLUDED defined */
#define FIND_MAX(a, b)
Definition: solid_mechanics_math_utilities.hpp:18
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
Various mathematical utilitiy functions.
Definition: math_utils.h:62
static double Norm(const Vector &a)
Calculates the norm of vector "a".
Definition: math_utils.h:703
Point class.
Definition: point.h:59
Definition: solid_mechanics_math_utilities.hpp:34
static void TensorToVector(const Matrix &Tensor, Vector &Vector)
Definition: solid_mechanics_math_utilities.hpp:1227
static VectorType project(const VectorType &rVector, VectorType irange)
Definition: solid_mechanics_math_utilities.hpp:1747
static void TensorToMatrix(std::vector< std::vector< Matrix > > &Tensor, Matrix &Matrix)
Definition: solid_mechanics_math_utilities.hpp:1424
static MatrixType IdentityMatrix(SizeType size)
Definition: solid_mechanics_math_utilities.hpp:979
static void Normalize(VectorType &v)
Definition: solid_mechanics_math_utilities.hpp:1077
static MatrixType StrainVectorToTensor(const VectorType &Strains)
Definition: solid_mechanics_math_utilities.hpp:1092
static void QRFactorization(const MatrixType &A, MatrixType &Q, MatrixType &R)
Definition: solid_mechanics_math_utilities.hpp:307
static void EigenVectors(const MatrixType &A, MatrixType &vectors, VectorType &lambda, double zero_tolerance=1e-9, int max_iterations=10)
Definition: solid_mechanics_math_utilities.hpp:432
static void MatrixToTensor(MatrixType &A, array_1d< double, 81 > &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1379
static MatrixType project(const MatrixType &rMatrix, VectorType irange, VectorType jrange)
Definition: solid_mechanics_math_utilities.hpp:1723
static Vector EigenValues(const Matrix &A, double tolerance, double zero)
Definition: solid_mechanics_math_utilities.hpp:142
static void DeviatoricUnity(array_1d< double, 81 > &Unity)
Definition: solid_mechanics_math_utilities.hpp:1585
static Vector EigenValuesDirectMethod(const Matrix &A)
Definition: solid_mechanics_math_utilities.hpp:243
DenseVector< Vector > Second_Order_Tensor
Definition: solid_mechanics_math_utilities.hpp:50
Matrix MatrixType
Definition: solid_mechanics_math_utilities.hpp:40
unsigned int SizeType
Definition: solid_mechanics_math_utilities.hpp:46
static void DeviatoricUnity(std::vector< std::vector< Matrix > > &Unity)
Definition: solid_mechanics_math_utilities.hpp:1550
static VectorType range(const unsigned int start, const unsigned int end)
Definition: solid_mechanics_math_utilities.hpp:1708
static MatrixType Transpose(MatrixType &A)
Definition: solid_mechanics_math_utilities.hpp:1060
static void TensorToMatrix(const array_1d< double, 81 > &Tensor, Matrix &Matrix)
Definition: solid_mechanics_math_utilities.hpp:1498
static bool CardanoFormula(double a, double b, double c, double d, Vector &solution)
Definition: solid_mechanics_math_utilities.hpp:75
DenseMatrix< Second_Order_Tensor > Matrix_Second_Tensor
Definition: solid_mechanics_math_utilities.hpp:56
static void Mult(MatrixType &M, TDataType a)
Definition: solid_mechanics_math_utilities.hpp:1033
static MatrixType Mult(MatrixType &A, MatrixType &B)
Definition: solid_mechanics_math_utilities.hpp:1013
static void EigenVectors3x3(const Matrix &A, Matrix &V, Vector &d)
Definition: solid_mechanics_math_utilities.hpp:625
static Vector TensorToStrainVector(const Matrix &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1126
MathUtils< TDataType > MathUtilsType
Definition: solid_mechanics_math_utilities.hpp:48
static int InvertMatrix(const MatrixType &input, MatrixType &inverse)
Definition: solid_mechanics_math_utilities.hpp:1163
static void MatrixToTensor(MatrixType &A, std::vector< std::vector< Matrix > > &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1329
static bool Clipping(std::vector< Point * > &clipping_points, std::vector< Point * > &subjected_points, std::vector< Point * > &result_points, double tolerance)
Definition: solid_mechanics_math_utilities.hpp:1612
static void Add(MatrixType &A, MatrixType &B)
Definition: solid_mechanics_math_utilities.hpp:996
static void Mult(VectorType &v, TDataType a)
Definition: solid_mechanics_math_utilities.hpp:1047
static void EigenVectors(MatrixType &A, MatrixType &V, TDataType &error_tolerance, TDataType zero_tolerance)
Definition: solid_mechanics_math_utilities.hpp:890
unsigned int IndexType
Definition: solid_mechanics_math_utilities.hpp:44
Vector VectorType
Definition: solid_mechanics_math_utilities.hpp:42
static double NormTensor(Matrix &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1179
static void VectorToTensor(const Vector &Stress, Matrix &Tensor)
Definition: solid_mechanics_math_utilities.hpp:1196
static double hypot2(double x, double y)
Definition: solid_mechanics_math_utilities.hpp:873
DenseVector< Second_Order_Tensor > Third_Order_Tensor
Definition: solid_mechanics_math_utilities.hpp:52
static MatrixType & project(MatrixType &rMatrixA, VectorType irange, VectorType jrange, const MatrixType &rMatrixB)
Definition: solid_mechanics_math_utilities.hpp:1736
static void TensorToMatrix(Fourth_Order_Tensor &Tensor, Matrix &Matrix)
Definition: solid_mechanics_math_utilities.hpp:1251
DenseVector< DenseVector< Matrix > > Fourth_Order_Tensor
Definition: solid_mechanics_math_utilities.hpp:54
Short class definition.
Definition: array_1d.h:61
BOOST_UBLAS_INLINE const_iterator end() const
Definition: array_1d.h:611
BOOST_UBLAS_INLINE const_iterator begin() const
Definition: array_1d.h:606
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_TRY
Definition: define.h:109
#define KRATOS_ERROR
Definition: exception.h:161
start
Definition: DEM_benchmarks.py:17
end
Definition: DEM_benchmarks.py:180
iteration
Definition: DEM_benchmarks.py:172
constexpr double Pi
Definition: global_variables.h:25
solution
Definition: KratosDEM.py:9
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Internals::Matrix< double, AMatrix::dynamic, 1 > Vector
Definition: amatrix_interface.h:472
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
AMatrix::TransposeMatrix< const T > trans(const T &TheMatrix)
Definition: amatrix_interface.h:486
float error
Definition: PecletTest.py:102
int max_iterations
Definition: ProjectParameters.py:53
list coeff
Definition: bombardelli_test.py:41
y
Other simbols definition.
Definition: generate_axisymmetric_navier_stokes_element.py:54
v
Definition: generate_convection_diffusion_explicit_element.py:114
w
Definition: generate_convection_diffusion_explicit_element.py:108
q
Definition: generate_convection_diffusion_explicit_element.py:109
f
Definition: generate_convection_diffusion_explicit_element.py:112
V
Definition: generate_droplet_dynamics.py:256
h
Definition: generate_droplet_dynamics.py:91
input
Definition: generate_frictional_mortar_condition.py:435
int C
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:27
a
Definition: generate_stokes_twofluid_element.py:77
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
u
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:30
float gamma
Definition: generate_two_fluid_navier_stokes.py:131
c
Definition: generate_weakly_compressible_navier_stokes_element.py:108
phi
Definition: isotropic_damage_automatic_differentiation.py:168
R
Definition: isotropic_damage_automatic_differentiation.py:172
Stress
Definition: isotropic_damage_automatic_differentiation.py:135
tuple Q
Definition: isotropic_damage_automatic_differentiation.py:235
int d
Definition: ode_solve.py:397
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
int m
Definition: run_marine_rain_substepping.py:8
dummy
Definition: script.py:194
A
Definition: sensitivityMatrix.py:70
p
Definition: sensitivityMatrix.py:52
x
Definition: sensitivityMatrix.py:49
int dim
Definition: sensitivityMatrix.py:25
B
Definition: sensitivityMatrix.py:76
integer i
Definition: TensorModule.f:17
double precision, dimension(3, 3), public delta
Definition: TensorModule.f:16
integer l
Definition: TensorModule.f:17
zero
Definition: test_pureconvectionsolver_benchmarking.py:94
e
Definition: run_cpp_mpi_tests.py:31