207 typename PointType::Pointer pPoint2,
208 typename PointType::Pointer pPoint3,
209 typename PointType::Pointer pPoint4,
210 typename PointType::Pointer pPoint5,
211 typename PointType::Pointer pPoint6,
212 typename PointType::Pointer pPoint7,
213 typename PointType::Pointer pPoint8,
214 typename PointType::Pointer pPoint9 )
229 :
BaseType( ThisPoints, &msGeometryData )
238 ) :
BaseType(GeometryId, rThisPoints, &msGeometryData)
245 const std::string& rGeometryName,
247 ) :
BaseType(rGeometryName, rThisPoints, &msGeometryData)
332 template<
class TOtherPo
intType>
355 return typename BaseType::Pointer(
new Quadrilateral3D9(NewGeometryId, rThisPoints ) );
369 auto p_geometry =
typename BaseType::Pointer(
new Quadrilateral3D9( NewGeometryId, rGeometry.
Points() ) );
370 p_geometry->SetData(rGeometry.
GetData());
381 if ((LocalDirectionIndex == 0) || (LocalDirectionIndex == 1)) {
384 KRATOS_ERROR <<
"Possible direction index reaches from 0-1. Given direction index: "
385 << LocalDirectionIndex << std::endl;
407 return( std::sqrt(
d[0]*
d[0] +
d[1]*
d[1] +
d[2]*
d[2] ) );
437 KRATOS_WARNING(
"Quadrilateral3D9") <<
"Method not well defined. Replace with DomainSize() instead. This method preserves current behaviour but will be changed in June 2023 (returning error instead)" << std::endl;
471 const double Tolerance = std::numeric_limits<double>::epsilon()
476 if ( std::abs( rResult[0] ) <= (1.0 + Tolerance) )
478 if ( std::abs( rResult[1] ) <= (1.0 + Tolerance) )
498 const double tol = 1.0e-8;
499 const int maxiter = 1000;
501 std::vector< unsigned int> orientation( 3 );
518 switch ( orientation[0] )
545 if ( rResult.
size() != 2 )
546 rResult.
resize( 2,
false );
556 for (
int k = 0;
k < maxiter;
k++ )
560 noalias( CurrentGlobalCoords ) = rPoint - CurrentGlobalCoords;
566 Matrix shape_functions_gradients;
568 shape_functions_gradients, rResult );
575 J( 0, 0 ) += ( dummyPoint[orientation[0]] ) * ( shape_functions_gradients(
i, 0 ) );
576 J( 0, 1 ) += ( dummyPoint[orientation[0]] ) * ( shape_functions_gradients(
i, 1 ) );
577 J( 1, 0 ) += ( dummyPoint[orientation[1]] ) * ( shape_functions_gradients(
i, 0 ) );
578 J( 1, 1 ) += ( dummyPoint[orientation[1]] ) * ( shape_functions_gradients(
i, 1 ) );
582 double det_j =
J( 0, 0 ) *
J( 1, 1 ) -
J( 0, 1 ) *
J( 1, 0 );
585 invJ( 0, 0 ) = (
J( 1, 1 ) ) / ( det_j );
587 invJ( 1, 0 ) = -(
J( 1, 0 ) ) / ( det_j );
589 invJ( 0, 1 ) = -(
J( 0, 1 ) ) / ( det_j );
591 invJ( 1, 1 ) = (
J( 0, 0 ) ) / ( det_j );
593 DeltaXi( 0 ) = invJ( 0, 0 ) * CurrentGlobalCoords( orientation[0] ) + invJ( 0, 1 ) * CurrentGlobalCoords( orientation[1] );
595 DeltaXi( 1 ) = invJ( 1, 0 ) * CurrentGlobalCoords( orientation[0] ) + invJ( 1, 1 ) * CurrentGlobalCoords( orientation[1] );
606 if ( !( std::abs( CurrentGlobalCoords( orientation[2] ) ) <=
tol ) )
643 CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
645 Matrix shape_functions_values =
646 CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
649 if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
666 jacobian( 0, 0 ) += ( this->
GetPoint(
i ).X() ) * ( shape_functions_gradients[pnt](
i, 0 ) );
667 jacobian( 0, 1 ) += ( this->
GetPoint(
i ).X() ) * ( shape_functions_gradients[pnt](
i, 1 ) );
668 jacobian( 1, 0 ) += ( this->
GetPoint(
i ).Y() ) * ( shape_functions_gradients[pnt](
i, 0 ) );
669 jacobian( 1, 1 ) += ( this->
GetPoint(
i ).Y() ) * ( shape_functions_gradients[pnt](
i, 1 ) );
670 jacobian( 2, 0 ) += ( this->
GetPoint(
i ).Z() ) * ( shape_functions_gradients[pnt](
i, 0 ) );
671 jacobian( 2, 1 ) += ( this->
GetPoint(
i ).Z() ) * ( shape_functions_gradients[pnt](
i, 1 ) );
674 rResult[pnt] = jacobian;
699 Matrix & DeltaPosition )
const override
703 CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
705 Matrix shape_functions_values =
706 CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
709 if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
726 jacobian( 0, 0 ) += ( this->
GetPoint(
i ).X() - DeltaPosition(
i,0) ) * ( shape_functions_gradients[pnt](
i, 0 ) );
727 jacobian( 0, 1 ) += ( this->
GetPoint(
i ).X() - DeltaPosition(
i,0) ) * ( shape_functions_gradients[pnt](
i, 1 ) );
728 jacobian( 1, 0 ) += ( this->
GetPoint(
i ).Y() - DeltaPosition(
i,1) ) * ( shape_functions_gradients[pnt](
i, 0 ) );
729 jacobian( 1, 1 ) += ( this->
GetPoint(
i ).Y() - DeltaPosition(
i,1) ) * ( shape_functions_gradients[pnt](
i, 1 ) );
730 jacobian( 2, 0 ) += ( this->
GetPoint(
i ).Z() - DeltaPosition(
i,2) ) * ( shape_functions_gradients[pnt](
i, 0 ) );
731 jacobian( 2, 1 ) += ( this->
GetPoint(
i ).Z() - DeltaPosition(
i,2) ) * ( shape_functions_gradients[pnt](
i, 1 ) );
734 rResult[pnt] = jacobian;
766 if (rResult.size1() != 3 || rResult.size2() != 2)
767 rResult.
resize( 3, 2 ,
false);
771 CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
772 Matrix ShapeFunctionsGradientInIntegrationPoint = shape_functions_gradients( IntegrationPointIndex );
776 ShapeFunctionsValuesInIntegrationPoint =
row(
777 CalculateShapeFunctionsIntegrationPointsValues( ThisMethod ), IntegrationPointIndex );
784 rResult( 0, 0 ) += ( this->
GetPoint(
i ).X() ) * ( ShapeFunctionsGradientInIntegrationPoint(
i, 0 ) );
785 rResult( 0, 1 ) += ( this->
GetPoint(
i ).X() ) * ( ShapeFunctionsGradientInIntegrationPoint(
i, 1 ) );
786 rResult( 1, 0 ) += ( this->
GetPoint(
i ).Y() ) * ( ShapeFunctionsGradientInIntegrationPoint(
i, 0 ) );
787 rResult( 1, 1 ) += ( this->
GetPoint(
i ).Y() ) * ( ShapeFunctionsGradientInIntegrationPoint(
i, 1 ) );
788 rResult( 2, 0 ) += ( this->
GetPoint(
i ).Z() ) * ( ShapeFunctionsGradientInIntegrationPoint(
i, 0 ) );
789 rResult( 2, 1 ) += ( this->
GetPoint(
i ).Z() ) * ( ShapeFunctionsGradientInIntegrationPoint(
i, 1 ) );
812 if (rResult.size1() != 3 || rResult.size2() != 2)
813 rResult.
resize( 3, 2 ,
false);
816 Matrix shape_functions_gradients;
818 shape_functions_gradients, rPoint );
824 const auto& coordinates = this->
GetPoint(
i).Coordinates();
825 rResult( 0, 0 ) += ( coordinates[0] ) * ( shape_functions_gradients(
i, 0 ) );
826 rResult( 0, 1 ) += ( coordinates[0] ) * ( shape_functions_gradients(
i, 1 ) );
827 rResult( 1, 0 ) += ( coordinates[1] ) * ( shape_functions_gradients(
i, 0 ) );
828 rResult( 1, 1 ) += ( coordinates[1] ) * ( shape_functions_gradients(
i, 1 ) );
829 rResult( 2, 0 ) += ( coordinates[2] ) * ( shape_functions_gradients(
i, 0 ) );
830 rResult( 2, 1 ) += ( coordinates[2] ) * ( shape_functions_gradients(
i, 1 ) );
893 double fx1 = 0.5 * ( rPoint[0] - 1 ) * rPoint[0];
894 double fx2 = 0.5 * ( rPoint[0] + 1 ) * rPoint[0];
895 double fx3 = 1 - rPoint[0] * rPoint[0];
896 double fy1 = 0.5 * ( rPoint[1] - 1 ) * rPoint[1];
897 double fy2 = 0.5 * ( rPoint[1] + 1 ) * rPoint[1];
898 double fy3 = 1 - rPoint[1] * rPoint[1];
900 switch ( ShapeFunctionIndex )
921 KRATOS_ERROR <<
"Wrong index of shape function!" << *
this << std::endl;
940 if(rResult.size() != 9) rResult.
resize(9,
false);
942 double fx1 = 0.5 * ( rCoordinates[0] - 1 ) * rCoordinates[0];
943 double fx2 = 0.5 * ( rCoordinates[0] + 1 ) * rCoordinates[0];
944 double fx3 = 1 - rCoordinates[0] * rCoordinates[0];
945 double fy1 = 0.5 * ( rCoordinates[1] - 1 ) * rCoordinates[1];
946 double fy2 = 0.5 * ( rCoordinates[1] + 1 ) * rCoordinates[1];
947 double fy3 = 1 - rCoordinates[1] * rCoordinates[1];
949 rResult[0] = fx1*fy1 ;
950 rResult[1] = fx2*fy1 ;
951 rResult[2] = fx2*fy2 ;
952 rResult[3] = fx1*fy2 ;
953 rResult[4] = fx3*fy1 ;
954 rResult[5] = fx2*fy3 ;
955 rResult[6] = fx3*fy2 ;
956 rResult[7] = fx1*fy3 ;
957 rResult[8] = fx3*fy3 ;
975 std::string
Info()
const override
977 return "2 dimensional quadrilateral with nine nodes in 3D space";
988 rOStream <<
"2 dimensional quadrilateral with nine nodes in 3D space";
1004 std::cout << std::endl;
1010 rOStream <<
" Jacobian in the origin\t : " << jacobian;
1022 = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1023 const int integration_points_number
1027 for (
int pnt = 0; pnt < integration_points_number; pnt++ )
1029 Result[pnt] = localGradients[pnt];
1043 = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1044 const int integration_points_number
1048 for (
int pnt = 0; pnt < integration_points_number; pnt++ )
1050 Result[pnt] = localGradients[pnt];
1067 double fx1 = 0.5 * ( rPoint[0] - 1 ) * rPoint[0];
1068 double fx2 = 0.5 * ( rPoint[0] + 1 ) * rPoint[0];
1069 double fx3 = 1 - rPoint[0] * rPoint[0];
1070 double fy1 = 0.5 * ( rPoint[1] - 1 ) * rPoint[1];
1071 double fy2 = 0.5 * ( rPoint[1] + 1 ) * rPoint[1];
1072 double fy3 = 1 - rPoint[1] * rPoint[1];
1074 double gx1 = 0.5 * ( 2 * rPoint[0] - 1 );
1075 double gx2 = 0.5 * ( 2 * rPoint[0] + 1 );
1076 double gx3 = -2.0 * rPoint[0];
1077 double gy1 = 0.5 * ( 2 * rPoint[1] - 1 );
1078 double gy2 = 0.5 * ( 2 * rPoint[1] + 1 );
1079 double gy3 = -2.0 * rPoint[1];
1081 rResult.
resize( 9, 2,
false );
1083 rResult( 0, 0 ) = gx1 * fy1;
1084 rResult( 0, 1 ) = fx1 * gy1;
1085 rResult( 1, 0 ) = gx2 * fy1;
1086 rResult( 1, 1 ) = fx2 * gy1;
1087 rResult( 2, 0 ) = gx2 * fy2;
1088 rResult( 2, 1 ) = fx2 * gy2;
1089 rResult( 3, 0 ) = gx1 * fy2;
1090 rResult( 3, 1 ) = fx1 * gy2;
1091 rResult( 4, 0 ) = gx3 * fy1;
1092 rResult( 4, 1 ) = fx3 * gy1;
1093 rResult( 5, 0 ) = gx2 * fy3;
1094 rResult( 5, 1 ) = fx2 * gy3;
1095 rResult( 6, 0 ) = gx3 * fy2;
1096 rResult( 6, 1 ) = fx3 * gy2;
1097 rResult( 7, 0 ) = gx1 * fy3;
1098 rResult( 7, 1 ) = fx1 * gy3;
1099 rResult( 8, 0 ) = gx3 * fy3;
1100 rResult( 8, 1 ) = fx3 * gy3;
1112 rResult.
resize( 9, 2,
false );
1114 rResult( 0, 0 ) = -1.0;
1115 rResult( 0, 1 ) = -1.0;
1116 rResult( 1, 0 ) = 1.0;
1117 rResult( 1, 1 ) = -1.0;
1118 rResult( 2, 0 ) = 1.0;
1119 rResult( 2, 1 ) = 1.0;
1120 rResult( 3, 0 ) = -1.0;
1121 rResult( 3, 1 ) = 1.0;
1122 rResult( 4, 0 ) = 0.0;
1123 rResult( 4, 1 ) = -1.0;
1124 rResult( 5, 0 ) = 1.0;
1125 rResult( 5, 1 ) = 0.0;
1126 rResult( 6, 0 ) = 0.0;
1127 rResult( 6, 1 ) = 1.0;
1128 rResult( 7, 0 ) = -1.0;
1129 rResult( 7, 1 ) = 0.0;
1130 rResult( 8, 0 ) = 0.0;
1131 rResult( 8, 1 ) = 0.0;
1145 double fx1 = 0.5 * ( rPoint.
X() - 1 ) * rPoint.
X();
1146 double fx2 = 0.5 * ( rPoint.
X() + 1 ) * rPoint.
X();
1147 double fx3 = 1 - rPoint.
X() * rPoint.
X();
1148 double fy1 = 0.5 * ( rPoint.
Y() - 1 ) * rPoint.
Y();
1149 double fy2 = 0.5 * ( rPoint.
Y() + 1 ) * rPoint.
Y();
1150 double fy3 = 1 - rPoint.
Y() * rPoint.
Y();
1152 double gx1 = 0.5 * ( 2 * rPoint.
X() - 1 );
1153 double gx2 = 0.5 * ( 2 * rPoint.
X() + 1 );
1154 double gx3 = -2.0 * rPoint.
X();
1155 double gy1 = 0.5 * ( 2 * rPoint.
Y() - 1 );
1156 double gy2 = 0.5 * ( 2 * rPoint.
Y() + 1 );
1157 double gy3 = -2.0 * rPoint.
Y();
1159 rResult.
resize( 9, 2,
false );
1161 rResult( 0, 0 ) = gx1 * fy1;
1162 rResult( 0, 1 ) = fx1 * gy1;
1163 rResult( 1, 0 ) = gx2 * fy1;
1164 rResult( 1, 1 ) = fx2 * gy1;
1165 rResult( 2, 0 ) = gx2 * fy2;
1166 rResult( 2, 1 ) = fx2 * gy2;
1167 rResult( 3, 0 ) = gx1 * fy2;
1168 rResult( 3, 1 ) = fx1 * gy2;
1169 rResult( 4, 0 ) = gx3 * fy1;
1170 rResult( 4, 1 ) = fx3 * gy1;
1171 rResult( 5, 0 ) = gx2 * fy3;
1172 rResult( 5, 1 ) = fx2 * gy3;
1173 rResult( 6, 0 ) = gx3 * fy2;
1174 rResult( 6, 1 ) = fx3 * gy2;
1175 rResult( 7, 0 ) = gx1 * fy3;
1176 rResult( 7, 1 ) = fx1 * gy3;
1177 rResult( 8, 0 ) = gx3 * fy3;
1178 rResult( 8, 1 ) = fx3 * gy3;
1191 if ( rResult.size() != this->PointsNumber() )
1200 rResult[
i].
resize( 2, 2,
false );
1204 double fx1 = 0.5 * ( rPoint[0] - 1 ) * rPoint[0];
1205 double fx2 = 0.5 * ( rPoint[0] + 1 ) * rPoint[0];
1206 double fx3 = 1 - rPoint[0] * rPoint[0];
1207 double fy1 = 0.5 * ( rPoint[1] - 1 ) * rPoint[1];
1208 double fy2 = 0.5 * ( rPoint[1] + 1 ) * rPoint[1];
1209 double fy3 = 1 - rPoint[1] * rPoint[1];
1211 double gx1 = 0.5 * ( 2 * rPoint[0] - 1 );
1212 double gx2 = 0.5 * ( 2 * rPoint[0] + 1 );
1213 double gx3 = -2.0 * rPoint[0];
1214 double gy1 = 0.5 * ( 2 * rPoint[1] - 1 );
1215 double gy2 = 0.5 * ( 2 * rPoint[1] + 1 );
1216 double gy3 = -2.0 * rPoint[1];
1225 rResult[0]( 0, 0 ) = hx1 * fy1;
1226 rResult[0]( 0, 1 ) = gx1 * gy1;
1227 rResult[0]( 1, 0 ) = gx1 * gy1;
1228 rResult[0]( 1, 1 ) = fx1 * hy1;
1230 rResult[1]( 0, 0 ) = hx2 * fy1;
1231 rResult[1]( 0, 1 ) = gx2 * gy1;
1232 rResult[1]( 1, 0 ) = gx2 * gy1;
1233 rResult[1]( 1, 1 ) = fx2 * hy1;
1235 rResult[2]( 0, 0 ) = hx2 * fy2;
1236 rResult[2]( 0, 1 ) = gx2 * gy2;
1237 rResult[2]( 1, 0 ) = gx2 * gy2;
1238 rResult[2]( 1, 1 ) = fx2 * hy2;
1240 rResult[3]( 0, 0 ) = hx1 * fy2;
1241 rResult[3]( 0, 1 ) = gx1 * gy2;
1242 rResult[3]( 1, 0 ) = gx1 * gy2;
1243 rResult[3]( 1, 1 ) = fx1 * hy2;
1245 rResult[4]( 0, 0 ) = hx3 * fy1;
1246 rResult[4]( 0, 1 ) = gx3 * gy1;
1247 rResult[4]( 1, 0 ) = gx3 * gy1;
1248 rResult[4]( 1, 1 ) = fx3 * hy1;
1250 rResult[5]( 0, 0 ) = hx2 * fy3;
1251 rResult[5]( 0, 1 ) = gx2 * gy3;
1252 rResult[5]( 1, 0 ) = gx2 * gy3;
1253 rResult[5]( 1, 1 ) = fx2 * hy3;
1255 rResult[6]( 0, 0 ) = hx3 * fy2;
1256 rResult[6]( 0, 1 ) = gx3 * gy2;
1257 rResult[6]( 1, 0 ) = gx3 * gy2;
1258 rResult[6]( 1, 1 ) = fx3 * hy2;
1260 rResult[7]( 0, 0 ) = hx1 * fy3;
1261 rResult[7]( 0, 1 ) = gx1 * gy3;
1262 rResult[7]( 1, 0 ) = gx1 * gy3;
1263 rResult[7]( 1, 1 ) = fx1 * hy3;
1265 rResult[8]( 0, 0 ) = hx3 * fy3;
1266 rResult[8]( 0, 1 ) = gx3 * gy3;
1267 rResult[8]( 1, 0 ) = gx3 * gy3;
1268 rResult[8]( 1, 1 ) = fx3 * hy3;
1283 if ( rResult.size() != this->PointsNumber() )
1300 for (
unsigned int j = 0;
j < 2;
j++ )
1302 rResult[
i][
j].
resize( 2, 2,
false );
1314 double gx1 = 0.5 * ( 2 * rPoint[0] - 1 );
1316 double gx2 = 0.5 * ( 2 * rPoint[0] + 1 );
1318 double gx3 = -2.0 * rPoint[0];
1320 double gy1 = 0.5 * ( 2 * rPoint[1] - 1 );
1322 double gy2 = 0.5 * ( 2 * rPoint[1] + 1 );
1324 double gy3 = -2.0 * rPoint[1];
1338 rResult[0][0]( 0, 0 ) = 0.0;
1340 rResult[0][0]( 0, 1 ) = hx1 * gy1;
1342 rResult[0][0]( 1, 0 ) = hx1 * gy1;
1344 rResult[0][0]( 1, 1 ) = gx1 * hy1;
1346 rResult[0][1]( 0, 0 ) = hx1 * gy1;
1348 rResult[0][1]( 0, 1 ) = gx1 * hy1;
1350 rResult[0][1]( 1, 0 ) = gx1 * hy1;
1352 rResult[0][1]( 1, 1 ) = 0.0;
1354 rResult[1][0]( 0, 0 ) = 0.0;
1356 rResult[1][0]( 0, 1 ) = hx2 * gy1;
1358 rResult[1][0]( 1, 0 ) = hx2 * gy1;
1360 rResult[1][0]( 1, 1 ) = gx2 * hy1;
1362 rResult[1][1]( 0, 0 ) = hx2 * gy1;
1364 rResult[1][1]( 0, 1 ) = gx2 * hy1;
1366 rResult[1][1]( 1, 0 ) = gx2 * hy1;
1368 rResult[1][1]( 1, 1 ) = 0.0;
1370 rResult[2][0]( 0, 0 ) = 0.0;
1372 rResult[2][0]( 0, 1 ) = hx2 * gy2;
1374 rResult[2][0]( 1, 0 ) = hx2 * gy2;
1376 rResult[2][0]( 1, 1 ) = gx2 * hy2;
1378 rResult[2][1]( 0, 0 ) = hx2 * gy2;
1380 rResult[2][1]( 0, 1 ) = gx2 * hy2;
1382 rResult[2][1]( 1, 0 ) = gx2 * hy2;
1384 rResult[2][1]( 1, 1 ) = 0.0;
1386 rResult[3][0]( 0, 0 ) = 0.0;
1388 rResult[3][0]( 0, 1 ) = hx1 * gy2;
1390 rResult[3][0]( 1, 0 ) = hx1 * gy2;
1392 rResult[3][0]( 1, 1 ) = gx1 * hy2;
1394 rResult[3][1]( 0, 0 ) = hx1 * gy2;
1396 rResult[3][1]( 0, 1 ) = gx1 * hy2;
1398 rResult[3][1]( 1, 0 ) = gx1 * hy2;
1400 rResult[3][1]( 1, 1 ) = 0.0;
1402 rResult[4][0]( 0, 0 ) = 0.0;
1404 rResult[4][0]( 0, 1 ) = hx3 * gy1;
1406 rResult[4][0]( 1, 0 ) = hx3 * gy1;
1408 rResult[4][0]( 1, 1 ) = gx3 * hy1;
1410 rResult[4][1]( 0, 0 ) = hx3 * gy1;
1412 rResult[4][1]( 0, 1 ) = gx3 * hy1;
1414 rResult[4][1]( 1, 0 ) = gx3 * hy1;
1416 rResult[4][1]( 1, 1 ) = 0.0;
1418 rResult[5][0]( 0, 0 ) = 0.0;
1420 rResult[5][0]( 0, 1 ) = hx2 * gy3;
1422 rResult[5][0]( 1, 0 ) = hx2 * gy3;
1424 rResult[5][0]( 1, 1 ) = gx2 * hy3;
1426 rResult[5][1]( 0, 0 ) = hx2 * gy3;
1428 rResult[5][1]( 0, 1 ) = gx2 * hy3;
1430 rResult[5][1]( 1, 0 ) = gx2 * hy3;
1432 rResult[5][1]( 1, 1 ) = 0.0;
1434 rResult[6][0]( 0, 0 ) = 0.0;
1436 rResult[6][0]( 0, 1 ) = hx3 * gy2;
1438 rResult[6][0]( 1, 0 ) = hx3 * gy2;
1440 rResult[6][0]( 1, 1 ) = gx3 * hy2;
1442 rResult[6][1]( 0, 0 ) = hx3 * gy2;
1444 rResult[6][1]( 0, 1 ) = gx3 * hy2;
1446 rResult[6][1]( 1, 0 ) = gx3 * hy2;
1448 rResult[6][1]( 1, 1 ) = 0.0;
1450 rResult[7][0]( 0, 0 ) = 0.0;
1452 rResult[7][0]( 0, 1 ) = hx1 * gy3;
1454 rResult[7][0]( 1, 0 ) = hx1 * gy3;
1456 rResult[7][0]( 1, 1 ) = gx1 * hy3;
1458 rResult[7][1]( 0, 0 ) = hx1 * gy3;
1460 rResult[7][1]( 0, 1 ) = gx1 * hy3;
1462 rResult[7][1]( 1, 0 ) = gx1 * hy3;
1464 rResult[7][1]( 1, 1 ) = 0.0;
1466 rResult[8][0]( 0, 0 ) = 0.0;
1468 rResult[8][0]( 0, 1 ) = hx3 * gy3;
1470 rResult[8][0]( 1, 0 ) = hx3 * gy3;
1472 rResult[8][0]( 1, 1 ) = gx3 * hy3;
1474 rResult[8][1]( 0, 0 ) = hx3 * gy3;
1476 rResult[8][1]( 0, 1 ) = gx3 * hy3;
1478 rResult[8][1]( 1, 0 ) = gx3 * hy3;
1480 rResult[8][1]( 1, 1 ) = 0.0;
1510 void save(
Serializer& rSerializer )
const override
1540 static Matrix CalculateShapeFunctionsIntegrationPointsValues(
1546 const int integration_points_number = integration_points.size();
1548 const int points_number = 9;
1550 Matrix shape_function_values( integration_points_number, points_number );
1553 for (
int pnt = 0; pnt < integration_points_number; pnt++ )
1555 double fx1 = 0.5 * ( integration_points[pnt].X() - 1 ) * integration_points[pnt].
X();
1556 double fx2 = 0.5 * ( integration_points[pnt].X() + 1 ) * integration_points[pnt].
X();
1557 double fx3 = 1 - integration_points[pnt].X() * integration_points[pnt].X();
1558 double fy1 = 0.5 * ( integration_points[pnt].Y() - 1 ) * integration_points[pnt].
Y();
1559 double fy2 = 0.5 * ( integration_points[pnt].Y() + 1 ) * integration_points[pnt].
Y();
1560 double fy3 = 1 - integration_points[pnt].Y() * integration_points[pnt].Y();
1562 shape_function_values( pnt, 0 ) = ( fx1 * fy1 );
1563 shape_function_values( pnt, 1 ) = ( fx2 * fy1 );
1564 shape_function_values( pnt, 2 ) = ( fx2 * fy2 );
1565 shape_function_values( pnt, 3 ) = ( fx1 * fy2 );
1566 shape_function_values( pnt, 4 ) = ( fx3 * fy1 );
1567 shape_function_values( pnt, 5 ) = ( fx2 * fy3 );
1568 shape_function_values( pnt, 6 ) = ( fx3 * fy2 );
1569 shape_function_values( pnt, 7 ) = ( fx1 * fy3 );
1570 shape_function_values( pnt, 8 ) = ( fx3 * fy3 );
1573 return shape_function_values;
1594 const int integration_points_number = integration_points.size();
1600 for (
int pnt = 0; pnt < integration_points_number; pnt++ )
1602 double fx1 = 0.5 * ( integration_points[pnt].X() - 1 ) * integration_points[pnt].
X();
1603 double fx2 = 0.5 * ( integration_points[pnt].X() + 1 ) * integration_points[pnt].
X();
1604 double fx3 = 1 - integration_points[pnt].X() * integration_points[pnt].X();
1605 double fy1 = 0.5 * ( integration_points[pnt].Y() - 1 ) * integration_points[pnt].
Y();
1606 double fy2 = 0.5 * ( integration_points[pnt].Y() + 1 ) * integration_points[pnt].
Y();
1607 double fy3 = 1 - integration_points[pnt].Y() * integration_points[pnt].Y();
1609 double gx1 = 0.5 * ( 2 * integration_points[pnt].X() - 1 );
1610 double gx2 = 0.5 * ( 2 * integration_points[pnt].X() + 1 );
1611 double gx3 = -2.0 * integration_points[pnt].X();
1612 double gy1 = 0.5 * ( 2 * integration_points[pnt].Y() - 1 );
1613 double gy2 = 0.5 * ( 2 * integration_points[pnt].Y() + 1 );
1614 double gy3 = -2.0 * integration_points[pnt].Y();
1617 result( 0, 0 ) = gx1 * fy1;
1618 result( 0, 1 ) = fx1 * gy1;
1619 result( 1, 0 ) = gx2 * fy1;
1620 result( 1, 1 ) = fx2 * gy1;
1621 result( 2, 0 ) = gx2 * fy2;
1622 result( 2, 1 ) = fx2 * gy2;
1623 result( 3, 0 ) = gx1 * fy2;
1624 result( 3, 1 ) = fx1 * gy2;
1625 result( 4, 0 ) = gx3 * fy1;
1626 result( 4, 1 ) = fx3 * gy1;
1627 result( 5, 0 ) = gx2 * fy3;
1628 result( 5, 1 ) = fx2 * gy3;
1629 result( 6, 0 ) = gx3 * fy2;
1630 result( 6, 1 ) = fx3 * gy2;
1631 result( 7, 0 ) = gx1 * fy3;
1632 result( 7, 1 ) = fx1 * gy3;
1633 result( 8, 0 ) = gx3 * fy3;
1634 result( 8, 1 ) = fx3 * gy3;
1636 d_shape_f_values[pnt] = result;
1639 return d_shape_f_values;
1650 Quadrature < QuadrilateralGaussLegendreIntegrationPoints1,
1651 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1652 Quadrature < QuadrilateralGaussLegendreIntegrationPoints2,
1653 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1654 Quadrature < QuadrilateralGaussLegendreIntegrationPoints3,
1655 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1656 Quadrature < QuadrilateralGaussLegendreIntegrationPoints4,
1657 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1658 Quadrature < QuadrilateralGaussLegendreIntegrationPoints5,
1659 2, IntegrationPoint<3> >::GenerateIntegrationPoints()
1662 return integration_points;
1673 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1675 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1677 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1679 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1681 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1685 return shape_functions_values;
1696 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1698 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1700 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1702 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1704 Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1708 return shape_functions_local_gradients;
1730 std::istream& rIStream,
1737 std::ostream& rOStream,
1741 rOStream << std::endl;
1746 template<
class TPo
intType>
1747 const GeometryData Quadrilateral3D9<TPointType>::msGeometryData(
1750 Quadrilateral3D9<TPointType>::AllIntegrationPoints(),
1751 Quadrilateral3D9<TPointType>::AllShapeFunctionsValues(),
1752 AllShapeFunctionsLocalGradients()
1755 template<
class TPo
intType>
Definition: geometry_data.h:60
KratosGeometryType
Definition: geometry_data.h:110
@ Kratos_Quadrilateral3D9
SizeType IntegrationPointsNumber() const
Definition: geometry_data.h:430
IntegrationMethod
Definition: geometry_data.h:76
KratosGeometryFamily
Definition: geometry_data.h:91
IntegrationMethod DefaultIntegrationMethod() const
Definition: geometry_data.h:425
Definition: geometry_dimension.h:42
Geometry base class.
Definition: geometry.h:71
SizeType PointsNumber() const
Definition: geometry.h:528
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: geometry.h:3834
Geometry & operator=(const Geometry &rOther)
Definition: geometry.h:400
virtual CoordinatesArrayType & GlobalCoordinates(CoordinatesArrayType &rResult, CoordinatesArrayType const &LocalCoordinates) const
Definition: geometry.h:2377
std::vector< IntegrationPointType > IntegrationPointsArrayType
Definition: geometry.h:161
DataValueContainer & GetData()
Definition: geometry.h:591
std::size_t SizeType
Definition: geometry.h:144
const TPointType::Pointer pGetPoint(const int Index) const
Definition: geometry.h:1790
std::size_t IndexType
Definition: geometry.h:137
std::array< Matrix, static_cast< int >GeometryData::IntegrationMethod::NumberOfIntegrationMethods)> ShapeFunctionsValuesContainerType
Definition: geometry.h:172
const PointsArrayType & Points() const
Definition: geometry.h:1768
bool AllPointsAreValid() const
Checks if the geometry points are valid Checks if the geometry points are valid from the pointer valu...
Definition: geometry.h:4093
std::array< IntegrationPointsArrayType, static_cast< int >GeometryData::IntegrationMethod::NumberOfIntegrationMethods)> IntegrationPointsContainerType
Definition: geometry.h:167
TPointType const & GetPoint(const int Index) const
Definition: geometry.h:1816
GeometryData::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: geometry.h:177
SizeType IntegrationPointsNumber() const
Definition: geometry.h:2257
GeometryData::IntegrationMethod IntegrationMethod
Definition: geometry.h:122
Short class definition.
Definition: integration_point.h:52
static double ComputeDomainSize(const TGeometryType &rGeometry)
This method calculates and returns the domain size of the geometry from any geometry in a generic man...
Definition: integration_utilities.h:63
Definition: amatrix_interface.h:41
void swap(Matrix &Other)
Definition: amatrix_interface.h:289
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
void clear()
Definition: amatrix_interface.h:284
An three node 3D line geometry with quadratic shape functions.
Definition: line_3d_3.h:66
Various mathematical utilitiy functions.
Definition: math_utils.h:62
This class defines the node.
Definition: node.h:65
Point class.
Definition: point.h:59
double Y() const
Definition: point.h:187
double X() const
Definition: point.h:181
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
void push_back(const TPointerType &x)
Definition: pointer_vector.h:270
A nine node 3D quadrilateral geometry with quadratic shape functions.
Definition: quadrilateral_3d_9.h:48
BaseType::Pointer Create(const IndexType NewGeometryId, const BaseType &rGeometry) const override
Creates a new geometry pointer.
Definition: quadrilateral_3d_9.h:364
JacobiansType & Jacobian(JacobiansType &rResult, IntegrationMethod ThisMethod, Matrix &DeltaPosition) const override
Definition: quadrilateral_3d_9.h:697
void PrintInfo(std::ostream &rOStream) const override
Definition: quadrilateral_3d_9.h:986
GeometryData::KratosGeometryType GetGeometryType() const override
Definition: quadrilateral_3d_9.h:294
Matrix & PointsLocalCoordinates(Matrix &rResult) const override
Definition: quadrilateral_3d_9.h:1110
GeometryData::KratosGeometryFamily GetGeometryFamily() const override
Definition: quadrilateral_3d_9.h:289
BaseType::NormalType NormalType
Definition: quadrilateral_3d_9.h:182
Matrix & ShapeFunctionsLocalGradients(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_3d_9.h:1064
KRATOS_CLASS_POINTER_DEFINITION(Quadrilateral3D9)
~Quadrilateral3D9() override
Definition: quadrilateral_3d_9.h:287
BaseType::IntegrationPointsContainerType IntegrationPointsContainerType
Definition: quadrilateral_3d_9.h:132
Quadrilateral3D9(typename PointType::Pointer pPoint1, typename PointType::Pointer pPoint2, typename PointType::Pointer pPoint3, typename PointType::Pointer pPoint4, typename PointType::Pointer pPoint5, typename PointType::Pointer pPoint6, typename PointType::Pointer pPoint7, typename PointType::Pointer pPoint8, typename PointType::Pointer pPoint9)
Definition: quadrilateral_3d_9.h:206
ShapeFunctionsSecondDerivativesType & ShapeFunctionsSecondDerivatives(ShapeFunctionsSecondDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_3d_9.h:1189
BaseType::JacobiansType JacobiansType
Definition: quadrilateral_3d_9.h:153
BaseType::GeometriesArrayType GeometriesArrayType
Definition: quadrilateral_3d_9.h:78
GeometriesArrayType GenerateEdges() const override
This method gives you all edges of this geometry.
Definition: quadrilateral_3d_9.h:864
Quadrilateral3D9 & operator=(Quadrilateral3D9< TOtherPointType > const &rOther)
Definition: quadrilateral_3d_9.h:333
GeometryData::IntegrationMethod IntegrationMethod
Definition: quadrilateral_3d_9.h:72
Line3D3< TPointType > EdgeType
Definition: quadrilateral_3d_9.h:62
BaseType::CoordinatesArrayType CoordinatesArrayType
Definition: quadrilateral_3d_9.h:88
TPointType PointType
Definition: quadrilateral_3d_9.h:83
SizeType PointsNumberInDirection(IndexType LocalDirectionIndex) const override
Returns number of points per direction.
Definition: quadrilateral_3d_9.h:379
virtual ShapeFunctionsGradientsType ShapeFunctionsLocalGradients()
Definition: quadrilateral_3d_9.h:1039
Quadrilateral3D9(Quadrilateral3D9< TOtherPointType > const &rOther)
Definition: quadrilateral_3d_9.h:278
Geometry< TPointType > BaseType
Definition: quadrilateral_3d_9.h:57
double ShapeFunctionValue(IndexType ShapeFunctionIndex, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_3d_9.h:890
double Length() const override
Definition: quadrilateral_3d_9.h:404
BaseType::IndexType IndexType
Definition: quadrilateral_3d_9.h:97
BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType
Definition: quadrilateral_3d_9.h:139
Quadrilateral3D9(const IndexType GeometryId, const PointsArrayType &rThisPoints)
Constructor with Geometry Id.
Definition: quadrilateral_3d_9.h:235
BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: quadrilateral_3d_9.h:160
BaseType::Pointer Create(const IndexType NewGeometryId, PointsArrayType const &rThisPoints) const override
Creates a new geometry pointer.
Definition: quadrilateral_3d_9.h:350
JacobiansType & Jacobian(JacobiansType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_3d_9.h:638
BaseType::IntegrationPointType IntegrationPointType
Definition: quadrilateral_3d_9.h:117
BaseType::IntegrationPointsArrayType IntegrationPointsArrayType
Definition: quadrilateral_3d_9.h:125
Matrix & Jacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_3d_9.h:761
virtual ShapeFunctionsGradientsType ShapeFunctionsLocalGradients(IntegrationMethod ThisMethod)
Definition: quadrilateral_3d_9.h:1018
friend class Quadrilateral3D9
Definition: quadrilateral_3d_9.h:1715
ShapeFunctionsThirdDerivativesType & ShapeFunctionsThirdDerivatives(ShapeFunctionsThirdDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_3d_9.h:1280
Matrix & Jacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_3d_9.h:809
BaseType::SizeType SizeType
Definition: quadrilateral_3d_9.h:104
BaseType::PointsArrayType PointsArrayType
Definition: quadrilateral_3d_9.h:110
Quadrilateral3D9(Quadrilateral3D9 const &rOther)
Definition: quadrilateral_3d_9.h:261
void PrintData(std::ostream &rOStream) const override
Definition: quadrilateral_3d_9.h:1000
SizeType EdgesNumber() const override
This method gives you number of all edges of this geometry.
Definition: quadrilateral_3d_9.h:851
bool IsInside(const CoordinatesArrayType &rPoint, CoordinatesArrayType &rResult, const double Tolerance=std::numeric_limits< double >::epsilon()) const override
Returns whether given arbitrary point is inside the Geometry and the respective local point for the g...
Definition: quadrilateral_3d_9.h:468
BaseType::ShapeFunctionsSecondDerivativesType ShapeFunctionsSecondDerivativesType
Definition: quadrilateral_3d_9.h:168
Quadrilateral3D9(const PointsArrayType &ThisPoints)
Definition: quadrilateral_3d_9.h:228
virtual Matrix & ShapeFunctionsGradients(Matrix &rResult, PointType &rPoint)
Definition: quadrilateral_3d_9.h:1143
Vector & ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType &rCoordinates) const override
Definition: quadrilateral_3d_9.h:938
Quadrilateral3D9(const std::string &rGeometryName, const PointsArrayType &rThisPoints)
Constructor with Geometry Name.
Definition: quadrilateral_3d_9.h:244
BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: quadrilateral_3d_9.h:146
CoordinatesArrayType & PointLocalCoordinates(CoordinatesArrayType &rResult, const CoordinatesArrayType &rPoint) const override
Returns the local coordinates of a given arbitrary point.
Definition: quadrilateral_3d_9.h:493
BaseType::ShapeFunctionsThirdDerivativesType ShapeFunctionsThirdDerivativesType
Definition: quadrilateral_3d_9.h:176
std::string Info() const override
Definition: quadrilateral_3d_9.h:975
Quadrilateral3D9 & operator=(const Quadrilateral3D9 &rOther)
Definition: quadrilateral_3d_9.h:314
double Area() const override
Definition: quadrilateral_3d_9.h:422
double DomainSize() const override
Definition: quadrilateral_3d_9.h:455
double Volume() const override
This method calculates and returns the volume of this geometry.
Definition: quadrilateral_3d_9.h:435
The serialization consists in storing the state of an object into a storage format like data file or ...
Definition: serializer.h:123
Short class definition.
Definition: array_1d.h:61
BOOST_UBLAS_INLINE void resize(size_type array_size, bool preserve=true)
Definition: array_1d.h:242
BOOST_UBLAS_INLINE size_type size() const
Definition: array_1d.h:370
#define KRATOS_SERIALIZE_SAVE_BASE_CLASS(Serializer, BaseType)
Definition: define.h:812
#define KRATOS_SERIALIZE_LOAD_BASE_CLASS(Serializer, BaseType)
Definition: define.h:815
#define KRATOS_ERROR
Definition: exception.h:161
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
#define KRATOS_WARNING(label)
Definition: logger.h:265
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
const GeometryData Quadrilateral3D9< TPointType >::msGeometryData & msGeometryDimension(), Quadrilateral3D9< TPointType >::AllShapeFunctionsValues(), AllShapeFunctionsLocalGradients()
Definition: brep_curve.h:483
const GeometryData BrepCurve< TContainerPointType, TContainerPointEmbeddedType >::msGeometryData & msGeometryDimension
Definition: brep_curve.h:483
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
AMatrix::MatrixRow< const TExpressionType > row(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression, std::size_t RowIndex)
Definition: amatrix_interface.h:649
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
int tol
Definition: hinsberg_optimization.py:138
def load(f)
Definition: ode_solve.py:307
int d
Definition: ode_solve.py:397
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
float temp
Definition: rotating_cone.py:85
dummy
Definition: script.py:194
J
Definition: sensitivityMatrix.py:58
integer i
Definition: TensorModule.f:17