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.
quadrilateral_2d_8.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Riccardo Rossi
11 // Janosch Stascheit
12 // Felix Nagel
13 // contributors: Hoang Giang Bui
14 // Josep Maria Carbonell
15 //
16 
17 #pragma once
18 
19 // System includes
20 
21 // External includes
22 
23 // Project includes
24 #include "geometries/line_2d_3.h"
27 
28 namespace Kratos
29 {
47 template<class TPointType> class Quadrilateral2D8
48  : public Geometry<TPointType>
49 {
50 public:
58 
63 
68 
73 
79 
83  typedef TPointType PointType;
84 
88  //typedef typename BaseType::CoordinatesArrayType CoordinatesArrayType;
89 
90 
97  typedef typename BaseType::IndexType IndexType;
98 
104  typedef typename BaseType::SizeType SizeType;
105 
111 
118 
126 
133 
140 
147 
154 
161 
169 
170 
178 
183 
188 
189 
193  Quadrilateral2D8( const PointType& Point1, const PointType& Point2,
194  const PointType& Point3, const PointType& Point4,
195  const PointType& Point5, const PointType& Point6,
196  const PointType& Point7, const PointType& Point8
197  )
198  : BaseType( PointsArrayType(), &msGeometryData )
199  {
200  this->Points().push_back( typename PointType::Pointer( new PointType( Point1 ) ) );
201  this->Points().push_back( typename PointType::Pointer( new PointType( Point2 ) ) );
202  this->Points().push_back( typename PointType::Pointer( new PointType( Point3 ) ) );
203  this->Points().push_back( typename PointType::Pointer( new PointType( Point4 ) ) );
204  this->Points().push_back( typename PointType::Pointer( new PointType( Point5 ) ) );
205  this->Points().push_back( typename PointType::Pointer( new PointType( Point6 ) ) );
206  this->Points().push_back( typename PointType::Pointer( new PointType( Point7 ) ) );
207  this->Points().push_back( typename PointType::Pointer( new PointType( Point8 ) ) );
208  }
209 
210  Quadrilateral2D8( typename PointType::Pointer pPoint1, typename PointType::Pointer pPoint2,
211  typename PointType::Pointer pPoint3, typename PointType::Pointer pPoint4,
212  typename PointType::Pointer pPoint5, typename PointType::Pointer pPoint6,
213  typename PointType::Pointer pPoint7, typename PointType::Pointer pPoint8 )
214  : BaseType( PointsArrayType(), &msGeometryData )
215  {
216  this->Points().push_back( pPoint1 );
217  this->Points().push_back( pPoint2 );
218  this->Points().push_back( pPoint3 );
219  this->Points().push_back( pPoint4 );
220  this->Points().push_back( pPoint5 );
221  this->Points().push_back( pPoint6 );
222  this->Points().push_back( pPoint7 );
223  this->Points().push_back( pPoint8 );
224  }
225 
226  Quadrilateral2D8( const PointsArrayType& ThisPoints )
227  : BaseType( ThisPoints, &msGeometryData )
228  {
229  if ( this->PointsNumber() != 8 )
230  KRATOS_ERROR << "Invalid points number. Expected 8, given " << this->PointsNumber() << std::endl;
231  }
232 
235  const IndexType GeometryId,
236  const PointsArrayType& rThisPoints
237  ) : BaseType(GeometryId, rThisPoints, &msGeometryData)
238  {
239  KRATOS_ERROR_IF( this->PointsNumber() != 8 ) << "Invalid points number. Expected 8, given " << this->PointsNumber() << std::endl;
240  }
241 
244  const std::string& rGeometryName,
245  const PointsArrayType& rThisPoints
246  ) : BaseType(rGeometryName, rThisPoints, &msGeometryData)
247  {
248  KRATOS_ERROR_IF(this->PointsNumber() != 8) << "Invalid points number. Expected 8, given " << this->PointsNumber() << std::endl;
249  }
250 
261  : BaseType( rOther )
262  {
263  }
264 
277  template<class TOtherPointType> Quadrilateral2D8( Quadrilateral2D8<TOtherPointType> const& rOther )
278  : BaseType( rOther )
279  {
280  }
281 
285  ~Quadrilateral2D8() override {}
286 
288  {
290  }
291 
293  {
295  }
296 
313  {
314  BaseType::operator=( rOther );
315 
316  return *this;
317  }
318 
330  template<class TOtherPointType>
332  {
333  BaseType::operator=( rOther );
334 
335  return *this;
336  }
337 
348  typename BaseType::Pointer Create(
349  const IndexType NewGeometryId,
350  PointsArrayType const& rThisPoints
351  ) const override
352  {
353  return typename BaseType::Pointer( new Quadrilateral2D8( NewGeometryId, rThisPoints ) );
354  }
355 
362  typename BaseType::Pointer Create(
363  const IndexType NewGeometryId,
364  const BaseType& rGeometry
365  ) const override
366  {
367  auto p_geometry = typename BaseType::Pointer( new Quadrilateral2D8( NewGeometryId, rGeometry.Points() ) );
368  p_geometry->SetData(rGeometry.GetData());
369  return p_geometry;
370  }
371 
377  SizeType PointsNumberInDirection(IndexType LocalDirectionIndex) const override
378  {
379  if ((LocalDirectionIndex == 0) || (LocalDirectionIndex == 1)) {
380  return 3;
381  }
382  KRATOS_ERROR << "Possible direction index reaches from 0-1. Given direction index: "
383  << LocalDirectionIndex << std::endl;
384  }
385 
403  double Length() const override
404  {
405  array_1d<double,3> point;
406  point[0] = 1.0/3.0; point[1] = 1.0/3.0; point[2] = 1.0/3.0;
407  return sqrt( fabs( DeterminantOfJacobian( point ) ) );
408  }
409 
422  double Area() const override
423  {
425  }
426 
427  // TODO: Code activated in June 2023
428  // /**
429  // * @brief This method calculates and returns the volume of this geometry.
430  // * @return Error, the volume of a 2D geometry is not defined
431  // * @see Length()
432  // * @see Area()
433  // * @see Volume()
434  // */
435  // double Volume() const override
436  // {
437  // KRATOS_ERROR << "Quadrilateral2D8:: Method not well defined. Replace with DomainSize() instead." << std::endl;
438  // return 0.0;
439  // }
440 
452  double DomainSize() const override
453  {
454  return Area();
455  }
456 
465  bool IsInside(
466  const CoordinatesArrayType& rPoint,
467  CoordinatesArrayType& rResult,
468  const double Tolerance = std::numeric_limits<double>::epsilon()
469  ) const override
470  {
471  this->PointLocalCoordinates( rResult, rPoint );
472 
473  if ( std::abs( rResult[0] ) <= (1.0 + Tolerance) )
474  {
475  if ( std::abs( rResult[1] ) <= (1.0 + Tolerance) )
476  {
477  return true;
478  }
479  }
480 
481  return false;
482  }
483 
506  IntegrationMethod ThisMethod ) const override
507  {
508  //getting derivatives of shape functions
509  ShapeFunctionsGradientsType shape_functions_gradients =
510  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
511  //getting values of shape functions
512  Matrix shape_functions_values =
513  CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
514  //workaround by riccardo...
515 
516  if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
517  {
518  // KLUDGE: While there is a bug in ublas
519  // vector resize, I have to put this beside resizing!!
520  JacobiansType temp( this->IntegrationPointsNumber( ThisMethod ) );
521  rResult.swap( temp );
522  }
523 
524  //loop over all integration points
525  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
526  {
527  //defining single jacobian matrix
528  Matrix jacobian = ZeroMatrix( 2, 2 );
529  //loop over all nodes
530 
531  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
532  {
533  jacobian( 0, 0 ) += ( this->GetPoint( i ).X() ) * ( shape_functions_gradients[pnt]( i, 0 ) );
534  jacobian( 0, 1 ) += ( this->GetPoint( i ).X() ) * ( shape_functions_gradients[pnt]( i, 1 ) );
535  jacobian( 1, 0 ) += ( this->GetPoint( i ).Y() ) * ( shape_functions_gradients[pnt]( i, 0 ) );
536  jacobian( 1, 1 ) += ( this->GetPoint( i ).Y() ) * ( shape_functions_gradients[pnt]( i, 1 ) );
537  }
538 
539  rResult[pnt] = jacobian;
540  }//end of loop over all integration points
541 
542  return rResult;
543  }
544 
564  IntegrationMethod ThisMethod,
565  Matrix & DeltaPosition ) const override
566  {
567  //getting derivatives of shape functions
568  ShapeFunctionsGradientsType shape_functions_gradients =
569  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
570  //getting values of shape functions
571  Matrix shape_functions_values =
572  CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
573  //workaround by riccardo...
574 
575  if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
576  {
577  // KLUDGE: While there is a bug in ublas
578  // vector resize, I have to put this beside resizing!!
579  JacobiansType temp( this->IntegrationPointsNumber( ThisMethod ) );
580  rResult.swap( temp );
581  }
582 
583  //loop over all integration points
584  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
585  {
586  //defining single jacobian matrix
587  Matrix jacobian = ZeroMatrix( 2, 2 );
588  //loop over all nodes
589 
590  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
591  {
592  jacobian( 0, 0 ) += ( this->GetPoint( i ).X() - DeltaPosition(i,0) ) * ( shape_functions_gradients[pnt]( i, 0 ) );
593  jacobian( 0, 1 ) += ( this->GetPoint( i ).X() - DeltaPosition(i,0) ) * ( shape_functions_gradients[pnt]( i, 1 ) );
594  jacobian( 1, 0 ) += ( this->GetPoint( i ).Y() - DeltaPosition(i,1) ) * ( shape_functions_gradients[pnt]( i, 0 ) );
595  jacobian( 1, 1 ) += ( this->GetPoint( i ).Y() - DeltaPosition(i,1) ) * ( shape_functions_gradients[pnt]( i, 1 ) );
596  }
597 
598  rResult[pnt] = jacobian;
599  }//end of loop over all integration points
600 
601  return rResult;
602  }
603 
624  Matrix& Jacobian( Matrix& rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod ) const override
625  {
626  // Setting up size of jacobian matrix
627  if (rResult.size1() != 2 || rResult.size2() != 2)
628  rResult.resize( 2, 2 , false);
629  rResult.clear();
630  //derivatives of shape functions
631  ShapeFunctionsGradientsType shape_functions_gradients = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
632  Matrix ShapeFunctionsGradientInIntegrationPoint = shape_functions_gradients( IntegrationPointIndex );
633  //values of shape functions in integration points
634  DenseVector<double> ShapeFunctionsValuesInIntegrationPoint = ZeroVector( 8 );
635  ShapeFunctionsValuesInIntegrationPoint =row( CalculateShapeFunctionsIntegrationPointsValues( ThisMethod ), IntegrationPointIndex );
636 
637  //Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
638  // Loop over all nodes
639  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
640  {
641  rResult( 0, 0 ) +=
642  ( this->GetPoint( i ).X() ) * ( ShapeFunctionsGradientInIntegrationPoint( i, 0 ) );
643  rResult( 0, 1 ) +=
644  ( this->GetPoint( i ).X() ) * ( ShapeFunctionsGradientInIntegrationPoint( i, 1 ) );
645  rResult( 1, 0 ) +=
646  ( this->GetPoint( i ).Y() ) * ( ShapeFunctionsGradientInIntegrationPoint( i, 0 ) );
647  rResult( 1, 1 ) +=
648  ( this->GetPoint( i ).Y() ) * ( ShapeFunctionsGradientInIntegrationPoint( i, 1 ) );
649  }
650 
651  return rResult;
652  }
653 
670  Matrix& Jacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
671  {
672  // Setting up size of jacobian matrix
673  if (rResult.size1() != 2 || rResult.size2() != 2)
674  rResult.resize( 2, 2 , false);
675  rResult.clear();
676  //derivatives of shape functions
677  Matrix shape_functions_gradients;
678  shape_functions_gradients = ShapeFunctionsLocalGradients( shape_functions_gradients, rPoint );
679  //Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
680  //loop over all nodes
681 
682  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
683  {
684  rResult( 0, 0 ) += ( this->GetPoint( i ).X() ) * ( shape_functions_gradients( i, 0 ) );
685  rResult( 0, 1 ) += ( this->GetPoint( i ).X() ) * ( shape_functions_gradients( i, 1 ) );
686  rResult( 1, 0 ) += ( this->GetPoint( i ).Y() ) * ( shape_functions_gradients( i, 0 ) );
687  rResult( 1, 1 ) += ( this->GetPoint( i ).Y() ) * ( shape_functions_gradients( i, 1 ) );
688  }
689 
690  return rResult;
691  }
692 
706  Vector& rResult,
707  IntegrationMethod ThisMethod
708  ) const override
709  {
710  if( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
711  rResult.resize( this->IntegrationPointsNumber( ThisMethod ), false );
712 
714  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ ) {
715  this->Jacobian( J, pnt, ThisMethod);
716  rResult[pnt] = (( J( 0, 0 )*J( 1, 1 ) ) - ( J( 0, 1 )*J( 1, 0 ) ) );
717  }
718  return rResult;
719  }
720 
743  double DeterminantOfJacobian( IndexType IntegrationPointIndex, IntegrationMethod ThisMethod ) const override
744  {
745  Matrix jacobian( 2, 2 );
746  jacobian = Jacobian( jacobian, IntegrationPointIndex, ThisMethod );
747  return(( jacobian( 0, 0 )*jacobian( 1, 1 ) ) - ( jacobian( 0, 1 )*jacobian( 1, 0 ) ) );
748  }
749 
767  double DeterminantOfJacobian( const CoordinatesArrayType& rPoint ) const override
768  {
769  Matrix jacobian( 2, 2 );
770  jacobian = Jacobian( jacobian, rPoint );
771  return(( jacobian( 0, 0 )*jacobian( 1, 1 ) ) - ( jacobian( 0, 1 )*jacobian( 1, 0 ) ) );
772  }
773 
794  JacobiansType& InverseOfJacobian( JacobiansType& rResult, IntegrationMethod ThisMethod ) const override
795  {
796  //workaround by riccardo
797  if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
798  {
799  // KLUDGE: While there is a bug in ublas
800  // vector resize, I have to put this beside resizing!!
801  JacobiansType temp( this->IntegrationPointsNumber( ThisMethod ) );
802  rResult.swap( temp );
803  }
804 
805  //loop over all integration points
806  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
807  {
808  Matrix tempMatrix( 2, 2 );
809  rResult[pnt] = InverseOfJacobian( tempMatrix, pnt, ThisMethod );
810  }
811 
812  return rResult;
813  }
814 
839  Matrix& InverseOfJacobian( Matrix& rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod ) const override
840  {
841  // Current jacobian
842  Matrix tempMatrix = ZeroMatrix( 2, 2 );
843  tempMatrix = Jacobian( tempMatrix, IntegrationPointIndex, ThisMethod );
844 
845  // Determinant of jacobian
846  const double det_j = DeterminantOfJacobian( IntegrationPointIndex, ThisMethod );
847 
848  // Checking for singularity
849  if ( det_j == 0.00 )
850  {
851  KRATOS_ERROR << "Zero determinant of jacobian." << *this << std::endl;
852  }
853 
854  // Setting up result matrix
855  rResult.resize( 2, 2 , false);
856 
857  // Filling matrix
858  rResult( 0, 0 ) = ( tempMatrix( 1, 1 ) ) / ( det_j );
859 
860  rResult( 1, 0 ) = -( tempMatrix( 1, 0 ) ) / ( det_j );
861 
862  rResult( 0, 1 ) = -( tempMatrix( 0, 1 ) ) / ( det_j );
863 
864  return rResult;
865  }
866 
884  Matrix& InverseOfJacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
885  {
886  //current jacobian
887  Matrix tempMatrix;
888  tempMatrix.resize( 2, 2, false );
889  tempMatrix = this->Jacobian( tempMatrix, rPoint );
890  //deteminant of Jacobian
891  double det_j = DeterminantOfJacobian( rPoint );
892  //checking for singularity
893 
894  if ( det_j == 0.0 )
895  {
896  KRATOS_ERROR << "Zero determinant of jacobian." << *this << std::endl;
897  }
898 
899  //setting up result matrix
900  rResult.resize( 2, 2 , false);
901 
902  //filling matrix
903  rResult( 0, 0 ) = ( tempMatrix( 1, 1 ) ) / ( det_j );
904 
905  rResult( 1, 0 ) = -( tempMatrix( 1, 0 ) ) / ( det_j );
906 
907  rResult( 0, 1 ) = -( tempMatrix( 0, 1 ) ) / ( det_j );
908 
909  rResult( 1, 1 ) = ( tempMatrix( 0, 0 ) ) / ( det_j );
910 
911  return rResult;
912  }
913 
917 
929  SizeType EdgesNumber() const override
930  {
931  return 4;
932  }
933 
943  {
945  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 0 ), this->pGetPoint( 1 ), this->pGetPoint( 4 ) ) );
946  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 1 ), this->pGetPoint( 2 ), this->pGetPoint( 5 ) ) );
947  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 2 ), this->pGetPoint( 3 ), this->pGetPoint( 6 ) ) );
948  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 3 ), this->pGetPoint( 0 ), this->pGetPoint( 7 ) ) );
949  return edges;
950  }
951 
968  double ShapeFunctionValue( IndexType ShapeFunctionIndex,
969  const CoordinatesArrayType& rPoint ) const override
970  {
971  switch ( ShapeFunctionIndex )
972  {
973  // Primary nodes
974  case 0:
975  return -(( 1.0 - rPoint[0] )*( 1.0 - rPoint[1] )
976  *( 1.0 + rPoint[0] + rPoint[1] ) ) / 4.0;
977  case 1:
978  return -(( 1.0 + rPoint[0] )*( 1.0 - rPoint[1] )
979  *( 1.0 - rPoint[0] + rPoint[1] ) ) / 4.0;
980  case 2:
981  return -(( 1.0 + rPoint[0] )*( 1.0 + rPoint[1] )
982  *( 1.0 - rPoint[0] - rPoint[1] ) ) / 4.0;
983  case 3:
984  return -(( 1.0 - rPoint[0] )*( 1.0 + rPoint[1] )
985  *( 1.0 + rPoint[0] - rPoint[1] ) ) / 4.0;
986  // Secondary nodes
987  case 4:
988  return (( 1.0 -rPoint[0]*rPoint[0] )
989  *( 1.0 - rPoint[1] ) ) / 2.0;
990  case 5:
991  return (( 1.0 + rPoint[0] )
992  *( 1.0 - rPoint[1]*rPoint[1] ) ) / 2.0 ;
993  case 6:
994  return (( 1.0 -rPoint[0]*rPoint[0] )
995  *( 1.0 + rPoint[1] ) ) / 2.0 ;
996  case 7:
997  return (( 1.0 -rPoint[0] )
998  *( 1.0 - rPoint[1]*rPoint[1] ) ) / 2.0 ;
999  default:
1000  KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
1001  }
1002 
1003  return 0;
1004  }
1005 
1015  Vector &rResult,
1016  const CoordinatesArrayType& rCoordinates
1017  ) const override
1018  {
1019  if(rResult.size() != 8)
1020  {
1021  rResult.resize(8,false);
1022  }
1023 
1024  // Primary nodes
1025  rResult[0] = -(( 1.0 - rCoordinates[0] )*( 1.0 - rCoordinates[1] )
1026  *( 1.0 + rCoordinates[0] + rCoordinates[1] ) ) / 4.0;
1027  rResult[1] = -(( 1.0 + rCoordinates[0] )*( 1.0 - rCoordinates[1] )
1028  *( 1.0 - rCoordinates[0] + rCoordinates[1] ) ) / 4.0;
1029  rResult[2] = -(( 1.0 + rCoordinates[0] )*( 1.0 + rCoordinates[1] )
1030  *( 1.0 - rCoordinates[0] - rCoordinates[1] ) ) / 4.0;
1031  rResult[3] = -(( 1.0 - rCoordinates[0] )*( 1.0 + rCoordinates[1] )
1032  *( 1.0 + rCoordinates[0] - rCoordinates[1] ) ) / 4.0;
1033 
1034  // Secondary nodes
1035  rResult[4] = (( 1.0 -rCoordinates[0]*rCoordinates[0] )
1036  *( 1.0 - rCoordinates[1] ) ) / 2.0;
1037  rResult[5] = (( 1.0 + rCoordinates[0] )
1038  *( 1.0 - rCoordinates[1]*rCoordinates[1] ) ) / 2.0 ;
1039  rResult[6] = (( 1.0 -rCoordinates[0]*rCoordinates[0] )
1040  *( 1.0 + rCoordinates[1] ) ) / 2.0 ;
1041  rResult[7] = (( 1.0 -rCoordinates[0] )
1042  *( 1.0 - rCoordinates[1]*rCoordinates[1] ) ) / 2.0 ;
1043 
1044  return rResult;
1045  }
1065  ShapeFunctionsGradientsType& rResult,
1066  IntegrationMethod ThisMethod ) const override
1067  {
1068  const unsigned int integration_points_number =
1069  msGeometryData.IntegrationPointsNumber( ThisMethod );
1070 
1071  if ( integration_points_number == 0 )
1072  {
1073  KRATOS_ERROR << "This integration method is not supported" << *this << std::endl;
1074  }
1075 
1076  //workaround by riccardo
1077  if ( rResult.size() != integration_points_number )
1078  {
1079  // KLUDGE: While there is a bug in ublas
1080  // vector resize, I have to put this beside resizing!!
1081  ShapeFunctionsGradientsType temp( integration_points_number );
1082  rResult.swap( temp );
1083  }
1084 
1085  //calculating the local gradients
1087  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1088 
1089  //getting the inverse jacobian matrices
1090  JacobiansType temp( integration_points_number );
1091 
1092  JacobiansType invJ = InverseOfJacobian( temp, ThisMethod );
1093 
1094  //loop over all integration points
1095  for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
1096  {
1097  rResult[pnt].resize( 4, 2, false );
1098 
1099  for ( unsigned int i = 0; i < 4; i++ )
1100  {
1101  for ( unsigned int j = 0; j < 2; j++ )
1102  {
1103  rResult[pnt]( i, j ) =
1104  ( locG[pnt]( i, 0 ) * invJ[pnt]( j, 0 ) )
1105  + ( locG[pnt]( i, 1 ) * invJ[pnt]( j, 1 ) );
1106  }
1107  }
1108  }//end of loop over integration points
1109  }
1110 
1121  std::string Info() const override
1122  {
1123  return "2 dimensional quadrilateral with eight nodes in 2D space";
1124  }
1125 
1132  void PrintInfo( std::ostream& rOStream ) const override
1133  {
1134  rOStream << "2 dimensional quadrilateral with eight nodes in 2D space";
1135  }
1136 
1146  void PrintData( std::ostream& rOStream ) const override
1147  {
1148  // Base Geometry class PrintData call
1149  BaseType::PrintData( rOStream );
1150  std::cout << std::endl;
1151 
1152  // If the geometry has valid points, calculate and output its data
1153  if (this->AllPointsAreValid()) {
1154  Matrix jacobian;
1155  this->Jacobian( jacobian, PointType() );
1156  rOStream << " Jacobian in the origin\t : " << jacobian;
1157  }
1158  }
1159 
1165  IntegrationMethod ThisMethod )
1166  {
1167  ShapeFunctionsGradientsType localGradients
1168  = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1169  const int integration_points_number
1170  = msGeometryData.IntegrationPointsNumber( ThisMethod );
1171  ShapeFunctionsGradientsType Result( integration_points_number );
1172 
1173  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1174  {
1175  Result[pnt] = localGradients[pnt];
1176  }
1177 
1178  return Result;
1179  }
1180 
1186  {
1187  IntegrationMethod ThisMethod = msGeometryData.DefaultIntegrationMethod();
1188  ShapeFunctionsGradientsType localGradients
1189  = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1190  const int integration_points_number
1191  = msGeometryData.IntegrationPointsNumber( ThisMethod );
1192  ShapeFunctionsGradientsType Result( integration_points_number );
1193 
1194  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1195  {
1196  Result[pnt] = localGradients[pnt];
1197  }
1198 
1199  return Result;
1200  }
1201 
1210  Matrix& rResult,
1211  const CoordinatesArrayType& rPoint
1212  ) const override
1213  {
1214  // Setting up result matrix
1215  rResult.resize( 8, 2 , false);
1216  noalias( rResult ) = ZeroMatrix( 8, 2 );
1217 
1218  // Primary nodes
1219  rResult( 0, 0 ) = - ((-1.0 + rPoint[1])*( 2.0 * rPoint[0] + rPoint[1])) / 4.0;
1220  rResult( 0, 1 ) = - ((-1.0 + rPoint[0])*( 2.0 * rPoint[1] + rPoint[0])) / 4.0;
1221 
1222  rResult( 1, 0 ) = ((-1.0 + rPoint[1])*(-2.0 * rPoint[0] + rPoint[1])) / 4.0;
1223  rResult( 1, 1 ) = (( 1.0 + rPoint[0])*( 2.0 * rPoint[1] - rPoint[0])) / 4.0;
1224 
1225  rResult( 2, 0 ) = (( 1.0 + rPoint[1])*( 2.0 * rPoint[0] + rPoint[1])) / 4.0;
1226  rResult( 2, 1 ) = (( 1.0 + rPoint[0])*( 2.0 * rPoint[1] + rPoint[0])) / 4.0;
1227 
1228  rResult( 3, 0 ) = - (( 1.0 + rPoint[1])*(-2.0 * rPoint[0] + rPoint[1])) / 4.0;
1229  rResult( 3, 1 ) = - ((-1.0 + rPoint[0])*( 2.0 * rPoint[1] - rPoint[0])) / 4.0;
1230 
1231  // Secondary nodes
1232  rResult( 4, 0 ) = rPoint[0] * (-1.0 + rPoint[1]);
1233  rResult( 4, 1 ) = ((1.0 + rPoint[0]) * (-1.0 + rPoint[0])) / 2.0;
1234 
1235  rResult( 5, 0 ) = - ((1.0 + rPoint[1]) * (-1.0 + rPoint[1])) / 2.0;
1236  rResult( 5, 1 ) = - rPoint[1] * ( 1.0 + rPoint[0]);
1237 
1238  rResult( 6, 0 ) = - rPoint[0] * ( 1.0 + rPoint[1]);
1239  rResult( 6, 1 ) = - ((1.0 + rPoint[0]) * (-1.0 + rPoint[0])) / 2.0;
1240 
1241  rResult( 7, 0 ) = ((1.0 + rPoint[1]) * (-1.0 + rPoint[1])) / 2.0;
1242  rResult( 7, 1 ) = rPoint[1] * (-1.0 + rPoint[0]);
1243 
1244  return( rResult );
1245  }
1246 
1252  Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
1253  {
1254  rResult.resize( 8, 2 , false);
1255  noalias( rResult ) = ZeroMatrix( 8, 2 );
1256  rResult( 0, 0 ) = -1.0;
1257  rResult( 0, 1 ) = -1.0;
1258  rResult( 1, 0 ) = 1.0;
1259  rResult( 1, 1 ) = -1.0;
1260  rResult( 2, 0 ) = 1.0;
1261  rResult( 2, 1 ) = 1.0;
1262  rResult( 3, 0 ) = -1.0;
1263  rResult( 3, 1 ) = 1.0;
1264  rResult( 4, 0 ) = 0.0;
1265  rResult( 4, 1 ) = -1.0;
1266  rResult( 5, 0 ) = 1.0;
1267  rResult( 5, 1 ) = 0.0;
1268  rResult( 6, 0 ) = 0.0;
1269  rResult( 6, 1 ) = 1.0;
1270  rResult( 7, 0 ) = -1.0;
1271  rResult( 7, 1 ) = 0.0;
1272  return rResult;
1273  }
1274 
1283  virtual Matrix& ShapeFunctionsGradients( Matrix& rResult, PointType& rPoint )
1284  {
1285  rResult.resize( 8, 2 , false);
1286  noalias( rResult ) = ZeroMatrix( 8, 2 );
1287 
1288  const CoordinatesArrayType& Coords = rPoint.Coordinates();
1289 
1290  rResult = ShapeFunctionsLocalGradients(rResult, Coords);
1291 
1292  return rResult;
1293  }
1294 
1302  {
1303  if ( rResult.size() != this->PointsNumber() )
1304  {
1305  // KLUDGE: While there is a bug in ublas
1306  // vector resize, I have to put this beside resizing!!
1308  rResult.swap( temp );
1309  }
1310 
1311  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
1312  {
1313  rResult[i].resize( 2, 2 , false);
1314  noalias( rResult[i] ) = ZeroMatrix( 2, 2 );
1315  }
1316 
1317  rResult[0]( 0, 0 ) = ( 4.0 - 4 * rPoint[1] ) / 8.0;
1318 
1319  rResult[0]( 0, 1 ) = ( -2.0 ) * ( 1.0 + 2.0 * rPoint[0] + rPoint[1] - 1.0 ) / 8.0
1320  + (( -1.0 + rPoint[1] ) * ( -2.0 ) ) / 8.0;
1321  rResult[0]( 1, 0 ) = ( -2.0 ) * ( 1.0 + rPoint[0] + 2.0 * rPoint[1] - 1.0 ) / 8.0
1322  + (( -1.0 + rPoint[0] ) * ( -2.0 ) ) / 8.0;
1323  rResult[0]( 1, 1 ) = (( -1.0 + rPoint[0] ) * ( -2.0 ) * ( 2.0 ) ) / 8.0;
1324 
1325  rResult[1]( 0, 0 ) = -(( -1.0 + rPoint[1] ) * ( -2.0 ) * ( -2.0 ) ) / 8.0;
1326  rResult[1]( 0, 1 ) = ( 2.0 ) * ( 1.0 - 2.0 * rPoint[0] + rPoint[1] - 1.0 ) / 8.0
1327  - (( -1.0 + rPoint[1] ) * ( -2.0 ) ) / 8.0;
1328  rResult[1]( 1, 0 ) = ( -1.0 + rPoint[0] - 2.0 * rPoint[1] + 1.0 ) * ( -2.0 ) / 8.0
1329  + (( 1.0 + rPoint[0] ) * ( -2.0 ) ) / 8.0;
1330  rResult[1]( 1, 1 ) = (( 1.0 + rPoint[0] ) * ( -2.0 ) * ( -2.0 ) ) / 8.0;
1331 
1332  rResult[2]( 0, 0 ) = -(( 1.0 + rPoint[1] ) * ( 2.0 ) * ( -2.0 ) ) / 8.0;
1333  rResult[2]( 0, 1 ) = -(( 1.0 ) * ( -1.0 + 2.0 * rPoint[0] + rPoint[1] + 1.0 ) * ( -2.0 ) ) / 8.0
1334  - (( 1.0 + rPoint[1] ) * ( -2.0 ) ) / 8.0;
1335  rResult[2]( 1, 0 ) = -(( 1.0 ) * ( -1.0 + rPoint[0] + 2.0 * rPoint[1] + 1.0 ) * ( -2.0 ) ) / 8.0
1336  - (( 1.0 + rPoint[0] ) * ( -2.0 ) ) / 8.0;
1337  rResult[2]( 1, 1 ) = -(( 1.0 + rPoint[0] ) * ( 2.0 ) * ( -2.0 ) ) / 8.0;
1338 
1339  rResult[3]( 0, 0 ) = (( 1.0 + rPoint[1] ) * ( -2.0 ) * ( -2.0 ) ) / 8.0;
1340  rResult[3]( 0, 1 ) = (( 1.0 ) * ( -1.0 - 2.0 * rPoint[0] + rPoint[1] + 1.0 ) * ( -2.0 ) ) / 8.0
1341  + (( 1.0 + rPoint[1] ) * ( -2.0 ) ) / 8.0;
1342  rResult[3]( 1, 0 ) = -(( -2.0 ) * ( 1.0 + rPoint[0] - 2.0 * rPoint[1] - 1.0 ) ) / 8.0
1343  - (( -1.0 + rPoint[0] ) * ( -2.0 ) ) / 8.0;
1344  rResult[3]( 1, 1 ) = -(( -1.0 + rPoint[0] ) * ( -2.0 ) * ( -2.0 ) ) / 8.0;
1345 
1346  rResult[4]( 0, 0 ) = -(( -1.0 + rPoint[1] ) * ( -2.0 ) ) / 2.0;
1347  rResult[4]( 0, 1 ) = -( rPoint[0] * ( -2.0 ) ) / 2.0;
1348  rResult[4]( 1, 0 ) = -(( 2.0 * rPoint[0] ) * ( -2.0 ) ) / 4.0;
1349  rResult[4]( 1, 1 ) = 0.0;
1350 
1351  rResult[5]( 0, 0 ) = 0.0;
1352  rResult[5]( 0, 1 ) = (( 2.0 * rPoint[1] ) * ( -2.0 ) ) / 4.0;
1353  rResult[5]( 1, 0 ) = ( rPoint[1] * ( -2.0 ) ) / 2.0;
1354  rResult[5]( 1, 1 ) = (( 1.0 + rPoint[0] ) * ( -2.0 ) ) / 2.0;
1355 
1356  rResult[6]( 0, 0 ) = (( 1.0 + rPoint[1] ) * ( -2.0 ) ) / 2.0;
1357  rResult[6]( 0, 1 ) = ( rPoint[0] * ( -2.0 ) ) / 2.0;
1358  rResult[6]( 1, 0 ) = (( 2.0 * rPoint[0] ) * ( -2.0 ) ) / 4.0;
1359  rResult[6]( 1, 1 ) = 0.0;
1360 
1361  rResult[7]( 0, 0 ) = 0.0;
1362  rResult[7]( 0, 1 ) = -(( 2.0 * rPoint[1] ) * ( -2.0 ) ) / 4.0;
1363  rResult[7]( 1, 0 ) = -( rPoint[1] * ( -2.0 ) ) / 2.0;
1364  rResult[7]( 1, 1 ) = -(( -1.0 + rPoint[0] ) * ( -2.0 ) ) / 2.0;
1365 
1366  return rResult;
1367  }
1368 
1378  const CoordinatesArrayType& rPoint
1379  ) const override
1380  {
1381 
1382  if ( rResult.size() != this->PointsNumber() )
1383  {
1384  // KLUDGE: While there is a bug in
1385  // ublas vector resize, I have to put this beside resizing!!
1386 // ShapeFunctionsGradientsType
1388  rResult.swap( temp );
1389  }
1390 
1391  for ( IndexType i = 0; i < rResult.size(); i++ )
1392  {
1394  rResult[i].swap( temp );
1395  }
1396 
1397  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
1398  {
1399  for ( int j = 0; j < 2; j++ )
1400  {
1401  rResult[i][j].resize( 2, 2 , false);
1402  noalias( rResult[i][j] ) = ZeroMatrix( 2, 2);
1403  }
1404  }
1405 
1406  rResult[0][0]( 0, 0 ) = 0.0;
1407 
1408  rResult[0][0]( 0, 1 ) = -0.5;
1409  rResult[0][0]( 1, 0 ) = -0.5;
1410  rResult[0][0]( 1, 1 ) = -0.5;
1411  rResult[0][1]( 0, 0 ) = -0.5;
1412  rResult[0][1]( 0, 1 ) = -0.5;
1413  rResult[0][1]( 1, 0 ) = -0.5;
1414  rResult[0][1]( 1, 1 ) = 0.0;
1415  rResult[1][0]( 0, 0 ) = 0.0;
1416  rResult[1][0]( 0, 1 ) = -0.5;
1417  rResult[1][0]( 1, 0 ) = -0.5;
1418  rResult[1][0]( 1, 1 ) = 0.5;
1419  rResult[1][1]( 0, 0 ) = -0.5;
1420  rResult[1][1]( 0, 1 ) = 0.5;
1421  rResult[1][1]( 1, 0 ) = 0.5;
1422  rResult[1][1]( 1, 1 ) = 0.0;
1423  rResult[2][0]( 0, 0 ) = 0.0;
1424  rResult[2][0]( 0, 1 ) = 0.5;
1425  rResult[2][0]( 1, 0 ) = 0.5;
1426  rResult[2][0]( 1, 1 ) = 0.5;
1427  rResult[2][1]( 0, 0 ) = 0.5;
1428  rResult[2][1]( 0, 1 ) = 0.5;
1429  rResult[2][1]( 1, 0 ) = 0.5;
1430  rResult[2][1]( 1, 1 ) = 0.0;
1431  rResult[3][0]( 0, 0 ) = 0.0;
1432  rResult[3][0]( 0, 1 ) = 0.5;
1433  rResult[3][0]( 1, 0 ) = 0.5;
1434  rResult[3][0]( 1, 1 ) = -0.5;
1435  rResult[3][1]( 0, 0 ) = 0.5;
1436  rResult[3][1]( 0, 1 ) = -0.5;
1437  rResult[3][1]( 1, 0 ) = -0.5;
1438  rResult[3][1]( 1, 1 ) = 0.0;
1439  rResult[4][0]( 0, 0 ) = 0.0;
1440  rResult[4][0]( 0, 1 ) = 1.0;
1441  rResult[4][0]( 1, 0 ) = 1.0;
1442  rResult[4][0]( 1, 1 ) = 0.0;
1443  rResult[4][1]( 0, 0 ) = 1.0;
1444  rResult[4][1]( 0, 1 ) = 0.0;
1445  rResult[4][1]( 1, 0 ) = 0.0;
1446  rResult[4][1]( 1, 1 ) = 0.0;
1447  rResult[5][0]( 0, 0 ) = 0.0;
1448  rResult[5][0]( 0, 1 ) = 0.0;
1449  rResult[5][0]( 1, 0 ) = 0.0;
1450  rResult[5][0]( 1, 1 ) = -1.0;
1451  rResult[5][1]( 0, 0 ) = 0.0;
1452  rResult[5][1]( 0, 1 ) = -1.0;
1453  rResult[5][1]( 1, 0 ) = 1.0;
1454  rResult[5][1]( 1, 1 ) = 0.0;
1455  rResult[6][0]( 0, 0 ) = 0.0;
1456  rResult[6][0]( 0, 1 ) = -1.0;
1457  rResult[6][0]( 1, 0 ) = -1.0;
1458  rResult[6][0]( 1, 1 ) = 0.0;
1459  rResult[6][1]( 0, 0 ) = -1.0;
1460  rResult[6][1]( 0, 1 ) = 0.0;
1461  rResult[6][1]( 1, 0 ) = 0.0;
1462  rResult[6][1]( 1, 1 ) = 0.0;
1463  rResult[7][0]( 0, 0 ) = 0.0;
1464  rResult[7][0]( 0, 1 ) = 0.0;
1465  rResult[7][0]( 1, 0 ) = 0.0;
1466  rResult[7][0]( 1, 1 ) = 1.0;
1467  rResult[7][1]( 0, 0 ) = 0.0;
1468  rResult[7][1]( 0, 1 ) = 1.0;
1469  rResult[7][1]( 1, 0 ) = -1.0;
1470  rResult[7][1]( 1, 0 ) = 0.0;
1471 
1472  return rResult;
1473  }
1474 
1475 protected:
1476 
1481 private:
1482 
1487  static const GeometryData msGeometryData;
1488 
1489  static const GeometryDimension msGeometryDimension;
1490 
1491 
1495 
1496  friend class Serializer;
1497 
1498  void save( Serializer& rSerializer ) const override
1499  {
1501  }
1502 
1503  void load( Serializer& rSerializer ) override
1504  {
1506  }
1507 
1508  Quadrilateral2D8(): BaseType( PointsArrayType(), &msGeometryData ) {}
1509 
1525  static Matrix CalculateShapeFunctionsIntegrationPointsValues( typename BaseType::IntegrationMethod ThisMethod )
1526  {
1527  IntegrationPointsContainerType all_integration_points = AllIntegrationPoints();
1528  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1529 
1530  // Number of integration points
1531  const unsigned int integration_points_number = integration_points.size();
1532 
1533  // Number of nodes in current geometry
1534  const unsigned int points_number = 8;
1535 
1536  // Setting up return matrix
1537  Matrix shape_function_values( integration_points_number, points_number );
1538 
1539  // Loop over all integration points
1540  for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
1541  {
1542  // Primary nodes
1543  row( shape_function_values, pnt )[0] =
1544  -(( 1.0 - integration_points[pnt].X() )
1545  * ( 1.0 - integration_points[pnt].Y() )
1546  * ( 1.0 + integration_points[pnt].X()
1547  + integration_points[pnt].Y() ) ) / 4.0;
1548  row( shape_function_values, pnt )[1] =
1549  -(( 1.0 + integration_points[pnt].X() )
1550  * ( 1.0 - integration_points[pnt].Y() ) * ( 1.0
1551  - integration_points[pnt].X()
1552  + integration_points[pnt].Y() ) ) / 4.0;
1553  row( shape_function_values, pnt )[2] =
1554  -(( 1.0 + integration_points[pnt].X() )
1555  * ( 1.0 + integration_points[pnt].Y() ) * ( 1.0
1556  - integration_points[pnt].X()
1557  - integration_points[pnt].Y() ) ) / 4.0;
1558  row( shape_function_values, pnt )[3] =
1559  -(( 1.0 - integration_points[pnt].X() ) * ( 1.0
1560  + integration_points[pnt].Y() ) * ( 1.0 ) * ( 1.0
1561  + integration_points[pnt].X()
1562  - integration_points[pnt].Y() ) ) / 4.0;
1563 
1564  // Secondary nodes
1565  row( shape_function_values, pnt )[4] =
1566  (( 1.0 - integration_points[pnt].X()
1567  * integration_points[pnt].X() )
1568  * ( 1.0 - integration_points[pnt].Y() ) ) / 2.0;
1569  row( shape_function_values, pnt )[5] =
1570  (( 1.0 + integration_points[pnt].X() )
1571  * ( 1.0 - integration_points[pnt].Y()
1572  * integration_points[pnt].Y() ) ) / 2.0 ;
1573  row( shape_function_values, pnt )[6] =
1574  (( 1.0 - integration_points[pnt].X()
1575  * integration_points[pnt].X() )
1576  * ( 1.0 + integration_points[pnt].Y() ) ) / 2.0 ;
1577  row( shape_function_values, pnt )[7] =
1578  (( 1.0 - integration_points[pnt].X() )
1579  * ( 1.0 - integration_points[pnt].Y()
1580  * integration_points[pnt].Y() ) ) / 2.0 ;
1581  }
1582 
1583  return shape_function_values;
1584  }
1585 
1598  CalculateShapeFunctionsIntegrationPointsLocalGradients(typename BaseType::IntegrationMethod ThisMethod )
1599  {
1600  IntegrationPointsContainerType all_integration_points = AllIntegrationPoints();
1601  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1602  // Number of integration points
1603  const unsigned int integration_points_number = integration_points.size();
1604  ShapeFunctionsGradientsType d_shape_f_values( integration_points_number );
1605  // Initialising container
1606  //std::fill(d_shape_f_values.begin(), d_shape_f_values.end(), Matrix(4,2));
1607 
1608  // Loop over all integration points
1609  for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
1610  {
1611  Matrix result = ZeroMatrix( 8, 2 );
1612 
1613  // Primary nodes
1614  result( 0, 0 ) = - ((-1.0 + integration_points[pnt].Y())*( 2.0 * integration_points[pnt].X() + integration_points[pnt].Y())) / 4.0;
1615  result( 0, 1 ) = - ((-1.0 + integration_points[pnt].X())*( 2.0 * integration_points[pnt].Y() + integration_points[pnt].X())) / 4.0;
1616 
1617  result( 1, 0 ) = ((-1.0 + integration_points[pnt].Y())*(-2.0 * integration_points[pnt].X() + integration_points[pnt].Y())) / 4.0;
1618  result( 1, 1 ) = (( 1.0 + integration_points[pnt].X())*( 2.0 * integration_points[pnt].Y() - integration_points[pnt].X())) / 4.0;
1619 
1620  result( 2, 0 ) = (( 1.0 + integration_points[pnt].Y())*( 2.0 * integration_points[pnt].X() + integration_points[pnt].Y())) / 4.0;
1621  result( 2, 1 ) = (( 1.0 + integration_points[pnt].X())*( 2.0 * integration_points[pnt].Y() + integration_points[pnt].X())) / 4.0;
1622 
1623  result( 3, 0 ) = - (( 1.0 + integration_points[pnt].Y())*(-2.0 * integration_points[pnt].X() + integration_points[pnt].Y())) / 4.0;
1624  result( 3, 1 ) = - ((-1.0 + integration_points[pnt].X())*( 2.0 * integration_points[pnt].Y() - integration_points[pnt].X())) / 4.0;
1625 
1626  // Secondary nodes
1627  result( 4, 0 ) = integration_points[pnt].X() * (-1.0 + integration_points[pnt].Y());
1628  result( 4, 1 ) = ((1.0 + integration_points[pnt].X()) * (-1.0 + integration_points[pnt].X())) / 2.0;
1629 
1630  result( 5, 0 ) = - ((1.0 + integration_points[pnt].Y()) * (-1.0 + integration_points[pnt].Y())) / 2.0;
1631  result( 5, 1 ) = - integration_points[pnt].Y() * ( 1.0 + integration_points[pnt].X());
1632 
1633  result( 6, 0 ) = - integration_points[pnt].X() * ( 1.0 + integration_points[pnt].Y());
1634  result( 6, 1 ) = - ((1.0 + integration_points[pnt].X()) * (-1.0 + integration_points[pnt].X())) / 2.0;
1635 
1636  result( 7, 0 ) = ((1.0 + integration_points[pnt].Y()) * (-1.0 + integration_points[pnt].Y())) / 2.0;
1637  result( 7, 1 ) = integration_points[pnt].Y() * (-1.0 + integration_points[pnt].X());
1638 
1639  d_shape_f_values[pnt] = result;
1640  }
1641 
1642  return d_shape_f_values;
1643  }
1644 
1648  static const IntegrationPointsContainerType AllIntegrationPoints()
1649  {
1650  IntegrationPointsContainerType integration_points =
1651  {
1652  {
1653  //Quadrature< QuadrilateralGaussLegendreIntegrationPointsNeu1,
1654  Quadrature < QuadrilateralGaussLegendreIntegrationPoints1,
1655  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1656  Quadrature < QuadrilateralGaussLegendreIntegrationPoints2,
1657  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1658  Quadrature < QuadrilateralGaussLegendreIntegrationPoints3,
1659  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1660  Quadrature < QuadrilateralGaussLegendreIntegrationPoints4,
1661  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1662  Quadrature < QuadrilateralGaussLegendreIntegrationPoints5,
1663  2, IntegrationPoint<3> >::GenerateIntegrationPoints()
1664  }
1665  };
1666  return integration_points;
1667  }
1668 
1672  static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
1673  {
1674  ShapeFunctionsValuesContainerType shape_functions_values =
1675  {
1676  {
1677  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1679  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1681  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1683  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1685  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1687  }
1688  };
1689  return shape_functions_values;
1690  }
1691 
1696  AllShapeFunctionsLocalGradients()
1697  {
1698  ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
1699  {
1700  {
1701  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(
1703  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(
1705  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(
1707  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(
1709  Quadrilateral2D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(
1711  }
1712  };
1713  return shape_functions_local_gradients;
1714  }
1715 
1720  template<class TOtherPointType> friend class Quadrilateral2D8;
1721 
1725 }; // Class Geometry
1726 
1733 template< class TPointType > inline std::istream& operator >> (
1734  std::istream& rIStream,
1736 
1740 template<class TPointType> inline std::ostream& operator << (
1741  std::ostream& rOStream,
1742  const Quadrilateral2D8<TPointType>& rThis )
1743 {
1744  rThis.PrintInfo( rOStream );
1745  rOStream << std::endl;
1746  rThis.PrintData( rOStream );
1747  return rOStream;
1748 }
1749 
1750 template<class TPointType>
1751 const GeometryData Quadrilateral2D8<TPointType>::msGeometryData(
1754  Quadrilateral2D8<TPointType>::AllIntegrationPoints(),
1755  Quadrilateral2D8<TPointType>::AllShapeFunctionsValues(),
1756  AllShapeFunctionsLocalGradients()
1757 );
1758 
1759 template<class TPointType>
1760 const GeometryDimension Quadrilateral2D8<TPointType>::msGeometryDimension(2, 2);
1761 
1762 } // namespace Kratos.
Definition: geometry_data.h:60
KratosGeometryType
Definition: geometry_data.h:110
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
std::vector< IntegrationPointType > IntegrationPointsArrayType
Definition: geometry.h:161
DataValueContainer & GetData()
Definition: geometry.h:591
SizeType LocalSpaceDimension() const
Definition: geometry.h:1300
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
virtual CoordinatesArrayType & PointLocalCoordinates(CoordinatesArrayType &rResult, const CoordinatesArrayType &rPoint) const
Returns the local coordinates of a given arbitrary point.
Definition: geometry.h:1854
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
SizeType WorkingSpaceDimension() const
Definition: geometry.h:1287
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 ComputeArea2DGeometry(const Geometry< TPointType > &rGeometry)
This method calculates and returns the volume of the geometry from a 3D geometry.
Definition: integration_utilities.h:107
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 2D line geometry with quadratic shape functions.
Definition: line_2d_3.h:63
This class defines the node.
Definition: node.h:65
CoordinatesArrayType const & Coordinates() const
Definition: point.h:215
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 eight node 2D quadrilateral geometry with quadratic shape functions.
Definition: quadrilateral_2d_8.h:49
double Length() const override
Definition: quadrilateral_2d_8.h:403
double DomainSize() const override
Definition: quadrilateral_2d_8.h:452
virtual Matrix & ShapeFunctionsGradients(Matrix &rResult, PointType &rPoint)
Definition: quadrilateral_2d_8.h:1283
Matrix & PointsLocalCoordinates(Matrix &rResult) const override
Definition: quadrilateral_2d_8.h:1252
GeometryData::KratosGeometryFamily GetGeometryFamily() const override
Definition: quadrilateral_2d_8.h:287
BaseType::JacobiansType JacobiansType
Definition: quadrilateral_2d_8.h:153
void PrintInfo(std::ostream &rOStream) const override
Definition: quadrilateral_2d_8.h:1132
Quadrilateral2D8 & operator=(const Quadrilateral2D8 &rOther)
Definition: quadrilateral_2d_8.h:312
BaseType::Pointer Create(const IndexType NewGeometryId, PointsArrayType const &rThisPoints) const override
Creates a new geometry pointer.
Definition: quadrilateral_2d_8.h:348
~Quadrilateral2D8() override
Definition: quadrilateral_2d_8.h:285
std::string Info() const override
Definition: quadrilateral_2d_8.h:1121
ShapeFunctionsThirdDerivativesType & ShapeFunctionsThirdDerivatives(ShapeFunctionsThirdDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:1376
Quadrilateral2D8(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)
Definition: quadrilateral_2d_8.h:210
BaseType::IntegrationPointType IntegrationPointType
Definition: quadrilateral_2d_8.h:117
Line2D3< TPointType > EdgeType
Definition: quadrilateral_2d_8.h:62
virtual ShapeFunctionsGradientsType ShapeFunctionsLocalGradients(IntegrationMethod ThisMethod)
Definition: quadrilateral_2d_8.h:1164
Quadrilateral2D8(Quadrilateral2D8 const &rOther)
Definition: quadrilateral_2d_8.h:260
BaseType::PointsArrayType PointsArrayType
Definition: quadrilateral_2d_8.h:110
BaseType::SizeType SizeType
Definition: quadrilateral_2d_8.h:104
Quadrilateral2D8(Quadrilateral2D8< TOtherPointType > const &rOther)
Definition: quadrilateral_2d_8.h:277
BaseType::IndexType IndexType
Definition: quadrilateral_2d_8.h:97
Matrix & Jacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:670
BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: quadrilateral_2d_8.h:146
Vector & ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType &rCoordinates) const override
Definition: quadrilateral_2d_8.h:1014
GeometryData::KratosGeometryType GetGeometryType() const override
Definition: quadrilateral_2d_8.h:292
double DeterminantOfJacobian(IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:743
void PrintData(std::ostream &rOStream) const override
Definition: quadrilateral_2d_8.h:1146
SizeType PointsNumberInDirection(IndexType LocalDirectionIndex) const override
Returns number of points per direction.
Definition: quadrilateral_2d_8.h:377
Matrix & ShapeFunctionsLocalGradients(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:1209
double ShapeFunctionValue(IndexType ShapeFunctionIndex, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:968
GeometriesArrayType GenerateEdges() const override
This method gives you all edges of this geometry.
Definition: quadrilateral_2d_8.h:942
double Area() const override
Definition: quadrilateral_2d_8.h:422
double DeterminantOfJacobian(const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:767
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_2d_8.h:465
KRATOS_CLASS_POINTER_DEFINITION(Quadrilateral2D8)
TPointType PointType
Definition: quadrilateral_2d_8.h:83
BaseType::CoordinatesArrayType CoordinatesArrayType
Definition: quadrilateral_2d_8.h:187
BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: quadrilateral_2d_8.h:160
Geometry< TPointType > BaseType
Definition: quadrilateral_2d_8.h:57
JacobiansType & Jacobian(JacobiansType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:505
JacobiansType & Jacobian(JacobiansType &rResult, IntegrationMethod ThisMethod, Matrix &DeltaPosition) const override
Definition: quadrilateral_2d_8.h:563
Quadrilateral2D8 & operator=(Quadrilateral2D8< TOtherPointType > const &rOther)
Definition: quadrilateral_2d_8.h:331
Quadrilateral2D8(const PointType &Point1, const PointType &Point2, const PointType &Point3, const PointType &Point4, const PointType &Point5, const PointType &Point6, const PointType &Point7, const PointType &Point8)
Definition: quadrilateral_2d_8.h:193
void ShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:1064
BaseType::NormalType NormalType
Definition: quadrilateral_2d_8.h:182
ShapeFunctionsSecondDerivativesType & ShapeFunctionsSecondDerivatives(ShapeFunctionsSecondDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:1301
JacobiansType & InverseOfJacobian(JacobiansType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:794
virtual ShapeFunctionsGradientsType ShapeFunctionsLocalGradients()
Definition: quadrilateral_2d_8.h:1185
Quadrilateral2D8(const PointsArrayType &ThisPoints)
Definition: quadrilateral_2d_8.h:226
SizeType EdgesNumber() const override
This method gives you number of all edges of this geometry.
Definition: quadrilateral_2d_8.h:929
GeometryData::IntegrationMethod IntegrationMethod
Definition: quadrilateral_2d_8.h:72
BaseType::IntegrationPointsContainerType IntegrationPointsContainerType
Definition: quadrilateral_2d_8.h:132
BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType
Definition: quadrilateral_2d_8.h:139
BaseType::Pointer Create(const IndexType NewGeometryId, const BaseType &rGeometry) const override
Creates a new geometry pointer.
Definition: quadrilateral_2d_8.h:362
Matrix & InverseOfJacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:839
BaseType::GeometriesArrayType GeometriesArrayType
Definition: quadrilateral_2d_8.h:78
Quadrilateral2D8(const std::string &rGeometryName, const PointsArrayType &rThisPoints)
Constructor with Geometry Name.
Definition: quadrilateral_2d_8.h:243
BaseType::ShapeFunctionsThirdDerivativesType ShapeFunctionsThirdDerivativesType
Definition: quadrilateral_2d_8.h:177
Vector & DeterminantOfJacobian(Vector &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:705
BaseType::IntegrationPointsArrayType IntegrationPointsArrayType
Definition: quadrilateral_2d_8.h:125
BaseType::ShapeFunctionsSecondDerivativesType ShapeFunctionsSecondDerivativesType
Definition: quadrilateral_2d_8.h:168
Matrix & Jacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_2d_8.h:624
Quadrilateral2D8(const IndexType GeometryId, const PointsArrayType &rThisPoints)
Constructor with Geometry Id.
Definition: quadrilateral_2d_8.h:234
Matrix & InverseOfJacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_2d_8.h:884
friend class Quadrilateral2D8
Definition: quadrilateral_2d_8.h:1720
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
#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
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
const GeometryData Quadrilateral2D8< TPointType >::msGeometryData & msGeometryDimension(), Quadrilateral2D8< TPointType >::AllShapeFunctionsValues(), AllShapeFunctionsLocalGradients()
Definition: brep_curve.h:483
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 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
def load(f)
Definition: ode_solve.py:307
int j
Definition: quadrature.py:648
float temp
Definition: rotating_cone.py:85
J
Definition: sensitivityMatrix.py:58
integer i
Definition: TensorModule.f:17