10 #if !defined(KRATOS_SOLID_MECHANICS_MATH_UTILITIES_H_INCLUDED)
11 #define KRATOS_SOLID_MECHANICS_MATH_UTILITIES_H_INCLUDED
18 #define FIND_MAX(a, b) ((a)>(b)?(a):(b))
32 template<
class TDataType>
82 std::cout<<
"This is not a cubic equation: CardanoFormula"<<std::endl;
87 double p= (3.0*
a*
c-
b*
b)/(3.0*
a*
a);
91 double discriminante=
p*
p*
p/27.0+
q*
q/4.0;
110 if(
p < 0)
return false;
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)
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)
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)
167 bool is_converged=
false;
168 unsigned int max_iters = 10;
169 unsigned int iter = 0;
171 while(is_converged ==
false && iter<max_iters )
173 double shift= HelpA((
dim-1),(
dim-1));
177 HelpA(
i,
i) = HelpA(
i,
i)- shift;
191 HelpA(
i,
j) += HelpR(
i,
k)*HelpQ(
k,
j);
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));
224 Result(
i)= HelpA(
i,
i);
226 if(fabs(Result(
i)) <
zero)
250 const double p1 =
A(0,1)*
A(0,1) +
A(0,2)*
A(0,2) +
A(1,2)*
A(1,2);
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);
264 const double inv_p = 1.0/
p;
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);
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) );
284 else if (r >= 1) {
phi = 0.0; }
285 else {
phi = acos(r) / 3.0;}
288 Result[0] =
q + 2.0 *
p * cos(
phi);
290 Result[1] = 3.0 *
q - Result[0] - Result[2];
334 std::vector<Matrix> HelpQ(
dim-1);
336 std::vector<Matrix> HelpR(
dim-1);
338 for(
int i=0;
i<
dim-1;
i++)
340 HelpQ[
i].resize(
dim,
dim,
false);
341 HelpR[
i].resize(
dim,
dim,
false);
361 double l= sqrt((normy*(normy+fabs(
y(
iteration))))/2);
395 for(
int k=0;
k<
dim-1;
k++)
404 for(
int k=1;
k<
dim-1;
k++)
409 Q(
i,
j)+= HelpQ[(
k-1)](
i,
l)*HelpQ[
k](
l,
j);
435 double zero_tolerance =1
e-9,
440 for(
int i=0;
i<3;
i++)
441 for(
int j=0;
j<3;
j++)
442 Help(
i,
j)= Help(
i,
j);
445 vectors.
resize(Help.size1(),Help.size2(),
false);
447 lambda.
resize(Help.size1(),
false);
449 Matrix HelpDummy(Help.size1(),Help.size2());
451 bool is_converged =
false;
453 Matrix unity(Help.size1(),Help.size2());
456 for(
unsigned int i=0;
i< Help.size1();
i++)
461 Matrix VDummy(Help.size1(),Help.size2());
463 Matrix Rotation(Help.size1(),Help.size2());
473 unsigned int index1= 0;
475 unsigned int index2= 1;
477 for(
unsigned int i=0;
i< Help.size1();
i++)
479 for(
unsigned int j=(
i+1);
j< Help.size2();
j++)
481 if((fabs(Help(
i,
j)) >
a ) && (fabs(Help(
i,
j)) > zero_tolerance))
500 double gamma= (Help(index2,index2)-Help(index1,index1))/(2*Help(index1,index2));
504 if(fabs(
gamma) > zero_tolerance && fabs(
gamma)< (1/zero_tolerance))
510 if (fabs(
gamma)>= (1.0/zero_tolerance))
514 double c= 1.0/(sqrt(1.0+
u*
u));
518 double teta= s/(1.0+
c);
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;
528 for(
unsigned int i=0;
i<Help.size1();
i++)
530 if((
i!= index1) && (
i!= index2))
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));
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));
545 Rotation(index2,index1)=-s;
546 Rotation(index1,index2)=s;
547 Rotation(index1,index1)=
c;
548 Rotation(index2,index2)=
c;
555 for(
unsigned int i=0;
i< Help.size1();
i++)
557 for(
unsigned int j=0;
j< Help.size1();
j++)
559 for(
unsigned int k=0;
k< Help.size1();
k++)
561 VDummy(
i,
j) +=
V(
i,
k)*Rotation(
k,
j);
595 std::cout<<
"########################################################"<<std::endl;
596 std::cout<<
"Max_Iterations exceed in Jacobi-Seidel-Iteration (eigenvectors)"<<std::endl;
597 std::cout<<
"########################################################"<<std::endl;
600 for(
unsigned int i=0;
i< Help.size1();
i++)
602 for(
unsigned int j=0;
j< Help.size1();
j++)
604 vectors(
i,
j)=
V(
j,
i);
608 for(
unsigned int i=0;
i<Help.size1();
i++)
609 lambda(
i)= Help(
i,
i);
628 if(
A.size1()!=3 ||
A.size2()!=3)
629 std::cout<<
" GIVEN MATRIX IS NOT 3x3 eigenvectors calculation "<<std::endl;
635 for (
int i = 0;
i < 3;
i++) {
636 for (
int j = 0;
j < 3;
j++) {
644 for (
int j = 0;
j < 3;
j++) {
650 for (
int i = 2;
i > 0;
i--) {
656 for (
int k = 0;
k <
i;
k++) {
657 scale = scale + fabs(
d[
k]);
661 for (
int j = 0;
j <
i;
j++) {
670 for (
int k = 0;
k <
i;
k++) {
682 for (
int j = 0;
j <
i;
j++) {
688 for (
int j = 0;
j <
i;
j++) {
692 for (
int k =
j+1;
k <=
i-1;
k++) {
699 for (
int j = 0;
j <
i;
j++) {
703 double hh =
f / (
h +
h);
704 for (
int j = 0;
j <
i;
j++) {
707 for (
int j = 0;
j <
i;
j++) {
710 for (
int k =
j;
k <=
i-1;
k++) {
722 for (
int i = 0;
i < 2;
i++) {
727 for (
int k = 0;
k <=
i;
k++) {
730 for (
int j = 0;
j <=
i;
j++) {
732 for (
int k = 0;
k <=
i;
k++) {
735 for (
int k = 0;
k <=
i;
k++) {
740 for (
int k = 0;
k <=
i;
k++) {
744 for (
int j = 0;
j < 3;
j++) {
755 for (
int i = 1;
i < 3;
i++) {
762 double eps = std::pow(2.0,-52.0);
763 for (
int l = 0;
l < 3;
l++) {
770 if (fabs(
e[
m]) <= eps*tst1) {
787 double p = (
d[
l+1] - g) / (2.0 *
e[
l]);
792 d[
l] =
e[
l] / (
p + r);
793 d[
l+1] =
e[
l] * (
p + r);
796 for (
int i =
l+2;
i < 3;
i++) {
810 for (
int i =
m-1;
i >=
l;
i--) {
820 p =
c *
d[
i] - s * g;
821 d[
i+1] =
h + s * (
c * g + s *
d[
i]);
825 for (
int k = 0;
k < 3;
k++) {
831 p = -s * s2 * c3 * el1 *
e[
l] / dl1;
837 }
while (fabs(
e[
l]) > eps*tst1);
845 for (
int i = 0;
i < 2;
i++) {
848 for (
int j =
i+1;
j < 3;
j++) {
857 for (
int j = 0;
j < 3;
j++) {
874 return std::sqrt(
x*
x+
y*
y);
892 TDataType& error_tolerance,
893 TDataType zero_tolerance)
896 TDataType
error = 1.0;
904 while(
error > error_tolerance )
906 for(
int i=0;
i<
n;
i++ )
908 for(
int j=
i+1;
j<
n;
j++ )
911 if( std::abs(
A(
i,
j) ) >= zero_tolerance )
913 if( std::abs(
A(
i,
i)-
A(
j,
j) ) > 0.0 )
915 theta = 0.5*atan(2*
A(
i,
j)/(
A(
i,
i)-
A(
j,
j)));
923 T(
i,
j) = -sin(theta);
935 for(
unsigned int i=0;
i<
A.size1();
i++ )
937 for(
unsigned int j=0;
j<
A.size2();
j++ )
939 sTot += std::abs(
A(
i,
j));
941 sDiag+= std::abs(
A(
i,
i));
943 error=(sTot-sDiag)/sDiag;
947 TDataType maxEv =
A(0,0);
948 for(
unsigned int i=0;
i<
A.size1();
i++ )
950 for(
unsigned int j=
i;
j<
A.size1();
j++ )
960 A(
i,
i) =
A(maxIndex,maxIndex);
961 A(maxIndex,maxIndex) =
dummy;
963 for(
unsigned int k=0;
k<
A.size2();
k++ )
966 V(
k,
i) =
V(
k,maxIndex);
983 for(
unsigned int i=0;
i<size ;
i++ )
998 for(
unsigned int i=0;
i<
A.size1();
i++ )
1000 for(
unsigned int j=0;
j<
A.size2();
j++ )
1017 for(
unsigned int i=0;
i<
A.size1();
i++ )
1019 for(
unsigned int j=0;
j<
B.size2();
j++ )
1021 for(
unsigned int k=0;
k<
A.size2();
k++ )
1035 for(
unsigned int i=0;
i<M.size1();
i++ )
1037 for(
unsigned int j=0;
j<M.size2();
j++ )
1049 for(
unsigned int i=0;
i<
v.size();
i++ )
1064 for(
unsigned int i=0;
i<
A.size1();
i++ )
1066 for(
unsigned int j=0;
j<
A.size2();
j++ )
1098 if (Strains.size()==3)
1100 StrainTensor.
resize(2,2,
false);
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];
1107 else if (Strains.size()==6)
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];
1122 return StrainTensor;
1132 if (Tensor.size1()==2)
1136 StrainVector(0) = Tensor(0,0);
1137 StrainVector(1) = Tensor(1,1);
1138 StrainVector(2) = 2.00*Tensor(0,1);
1140 else if (Tensor.size1()==3)
1142 StrainVector.resize(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);
1153 return StrainVector;
1165 typedef permutation_matrix<std::size_t> pmatrix;
1167 pmatrix pm(
A.size1());
1168 int singular = lu_factorize(
A,pm);
1170 lu_substitute(
A, pm, inverse);
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);
1186 result= sqrt(result);
1200 Tensor.
resize(3,3,
false);
1213 Tensor.
resize(2,2,
false);
1230 unsigned int dim = Tensor.size1();
1258 if (Tensor[0].size()== 3)
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);
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);
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);
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);
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);
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);
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);
1337 for(
unsigned int i=0;
i<3;
i++)
1339 Tensor[
i].resize(3);
1340 for(
unsigned int j=0;
j<3;
j++)
1342 Tensor[
i][
j].resize(3,3,
false);
1344 for(
unsigned int k=0;
k<3;
k++)
1345 for(
unsigned int l=0;
l<3;
l++)
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;
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;
1384 std::fill(Tensor.
begin(), Tensor.
end(), 0.0);
1385 for(
unsigned int i=0;
i<3;
i++)
1387 for(
unsigned int j=0;
j<3;
j++)
1389 for(
unsigned int k=0;
k<3;
k++)
1390 for(
unsigned int l=0;
l<3;
l++)
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;
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;
1412 Tensor[
i*27+
j*9+
k*3+
l]=
A(help1,help2)*
coeff;
1435 for(
unsigned int i=0;
i<6;
i++)
1436 for(
unsigned int j=0;
j<6;
j++)
1506 Matrix(0,3) = 2.0*Tensor[1];
1507 Matrix(0,4) = 2.0*Tensor[5];
1508 Matrix(0,5) = 2.0*Tensor[6];
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];
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];
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];
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];
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];
1556 for(
unsigned int i=0;
i<3;
i++)
1561 for(
unsigned int i=0;
i<3;
i++)
1564 for(
unsigned int j=0;
j<3;
j++)
1566 Unity[
i][
j].resize(3,3,
false);
1569 for(
unsigned int k=0;
k<3;
k++)
1571 for(
unsigned int l=0;
l<3;
l++)
1573 Unity[
i][
j](
k,
l)=kronecker(
i,
k)*kronecker(
j,
l)
1574 -1.0/3.0*kronecker(
i,
j)*kronecker(
k,
l);
1589 for(
unsigned int i=0;
i<3;
i++)
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);
1612 static bool Clipping(std::vector<Point* >& clipping_points,std::vector<Point* >& subjected_points, std::vector<Point* >& result_points,
double tolerance)
1614 result_points= subjected_points;
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++)
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;
1626 noalias(actual_edge)= *(clipping_points[clipp_edge])-*(clipping_points[index_clipp_2]);
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));
1635 actual_normal=*(clipping_points[1])-(*(clipping_points[clipp_edge])+actual_edge*
inner_prod((*(clipping_points[1])-*(clipping_points[clipp_edge])),actual_edge));
1637 noalias(actual_normal)=actual_normal/(sqrt(
inner_prod(actual_normal,actual_normal)));
1640 if(
inner_prod((*(result_points[0])-*(clipping_points[clipp_edge])), actual_normal)> tolerance)
1649 for(
unsigned int subj_edge=0; subj_edge< result_points.size(); subj_edge++)
1651 unsigned int index_subj_2=0;
1653 if(subj_edge< (result_points.size()-1))
1654 index_subj_2= subj_edge+1;
1657 if(fabs(
inner_prod((*(result_points[subj_edge])-(*(clipping_points[clipp_edge]))), actual_normal))<= tolerance)
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)
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])));
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])));
1675 A(1,1)=
inner_prod(*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge]),*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge]));
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));
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])));
1685 if(
coeff(0) > tolerance &&
coeff(0) < (1-tolerance)&& (sqrt(
inner_prod(dist_vec,dist_vec))< tolerance))
1688 temp_results.push_back(result_points[subj_edge]);
1690 temp_results.push_back(
new Point((*(clipping_points[clipp_edge])+
coeff(1)*(*(clipping_points[index_clipp_2])-*(clipping_points[clipp_edge])))));
1692 is_visible= !is_visible;
1697 temp_results.push_back(result_points[subj_edge]);
1699 result_points=temp_results;
1701 if(result_points.size()==0)
1716 for(
int i = 0;
i<size; ++
i)
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]);
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);
1752 for(
unsigned int i = 0;
i<irange.size(); ++
i)
1753 b[
i] = rVector[irange[
i]];
#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