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_3d_9.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_3d_3.h"
27 
28 namespace Kratos
29 {
47 template<class TPointType> class Quadrilateral3D9 : public Geometry<TPointType>
48 {
49 public:
58 
63 
68 
73 
79 
83  typedef TPointType PointType;
84 
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 
177 
178 
183 
188 // Quadrilateral3D9( const PointType& Point1, const PointType& Point2,
189 // const PointType& Point3, const PointType& Point4,
190 // const PointType& Point5, const PointType& Point6,
191 // const PointType& Point7, const PointType& Point8,
192 // const PointType& Point9 )
193 // : BaseType( PointsArrayType(), &msGeometryData )
194 // {
195 // this->Points().push_back( typename PointType::Pointer( new PointType( Point1 ) ) );
196 // this->Points().push_back( typename PointType::Pointer( new PointType( Point2 ) ) );
197 // this->Points().push_back( typename PointType::Pointer( new PointType( Point3 ) ) );
198 // this->Points().push_back( typename PointType::Pointer( new PointType( Point4 ) ) );
199 // this->Points().push_back( typename PointType::Pointer( new PointType( Point5 ) ) );
200 // this->Points().push_back( typename PointType::Pointer( new PointType( Point6 ) ) );
201 // this->Points().push_back( typename PointType::Pointer( new PointType( Point7 ) ) );
202 // this->Points().push_back( typename PointType::Pointer( new PointType( Point8 ) ) );
203 // this->Points().push_back( typename PointType::Pointer( new PointType( Point9 ) ) );
204 // }
205 
206  Quadrilateral3D9( typename PointType::Pointer pPoint1,
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 )
215  : BaseType( PointsArrayType(), &msGeometryData )
216  {
217  this->Points().push_back( pPoint1 );
218  this->Points().push_back( pPoint2 );
219  this->Points().push_back( pPoint3 );
220  this->Points().push_back( pPoint4 );
221  this->Points().push_back( pPoint5 );
222  this->Points().push_back( pPoint6 );
223  this->Points().push_back( pPoint7 );
224  this->Points().push_back( pPoint8 );
225  this->Points().push_back( pPoint9 );
226  }
227 
228  Quadrilateral3D9( const PointsArrayType& ThisPoints )
229  : BaseType( ThisPoints, &msGeometryData )
230  {
231  KRATOS_ERROR_IF( this->PointsNumber() != 9 ) << "Invalid points number. Expected 9, given " << this->PointsNumber() << std::endl;
232  }
233 
236  const IndexType GeometryId,
237  const PointsArrayType& rThisPoints
238  ) : BaseType(GeometryId, rThisPoints, &msGeometryData)
239  {
240  KRATOS_ERROR_IF( this->PointsNumber() != 9 ) << "Invalid points number. Expected 9, given " << this->PointsNumber() << std::endl;
241  }
242 
245  const std::string& rGeometryName,
246  const PointsArrayType& rThisPoints
247  ) : BaseType(rGeometryName, rThisPoints, &msGeometryData)
248  {
249  KRATOS_ERROR_IF(this->PointsNumber() != 9) << "Invalid points number. Expected 9, given " << this->PointsNumber() << std::endl;
250  }
251 
262  : BaseType( rOther )
263  {
264  }
265 
278  template<class TOtherPointType> Quadrilateral3D9(
279  Quadrilateral3D9<TOtherPointType> const& rOther )
280  : BaseType( rOther )
281  {
282  }
283 
287  ~Quadrilateral3D9() override {}
288 
290  {
292  }
293 
295  {
297  }
298 
315  {
316  BaseType::operator=( rOther );
317 
318  return *this;
319  }
320 
332  template<class TOtherPointType>
334  {
335  BaseType::operator=( rOther );
336 
337  return *this;
338  }
339 
343 
350  typename BaseType::Pointer Create(
351  const IndexType NewGeometryId,
352  PointsArrayType const& rThisPoints
353  ) const override
354  {
355  return typename BaseType::Pointer( new Quadrilateral3D9(NewGeometryId, rThisPoints ) );
356  }
357 
364  typename BaseType::Pointer Create(
365  const IndexType NewGeometryId,
366  const BaseType& rGeometry
367  ) const override
368  {
369  auto p_geometry = typename BaseType::Pointer( new Quadrilateral3D9( NewGeometryId, rGeometry.Points() ) );
370  p_geometry->SetData(rGeometry.GetData());
371  return p_geometry;
372  }
373 
379  SizeType PointsNumberInDirection(IndexType LocalDirectionIndex) const override
380  {
381  if ((LocalDirectionIndex == 0) || (LocalDirectionIndex == 1)) {
382  return 3;
383  }
384  KRATOS_ERROR << "Possible direction index reaches from 0-1. Given direction index: "
385  << LocalDirectionIndex << std::endl;
386  }
387 
404  double Length() const override
405  {
406  Vector d = this->Points()[2] - this->Points()[0];
407  return( std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ) );
408  }
409 
422  double Area() const override
423  {
424  const IntegrationMethod integration_method = msGeometryData.DefaultIntegrationMethod();
425  return IntegrationUtilities::ComputeDomainSize(*this, integration_method);
426  }
427 
435  double Volume() const override
436  {
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;
438  return Area();
439  // TODO: Replace in June 2023
440  // KRATOS_ERROR << "Quadrilateral3D9:: Method not well defined. Replace with DomainSize() instead." << std::endl;
441  // return 0.0;
442  }
443 
455  double DomainSize() const override
456  {
457  return Area();
458  }
459 
468  bool IsInside(
469  const CoordinatesArrayType& rPoint,
470  CoordinatesArrayType& rResult,
471  const double Tolerance = std::numeric_limits<double>::epsilon()
472  ) const override
473  {
474  PointLocalCoordinates( rResult, rPoint );
475 
476  if ( std::abs( rResult[0] ) <= (1.0 + Tolerance) )
477  {
478  if ( std::abs( rResult[1] ) <= (1.0 + Tolerance) )
479  {
480  return true;
481  }
482  }
483 
484  return false;
485  }
486 
494  CoordinatesArrayType& rResult,
495  const CoordinatesArrayType& rPoint
496  ) const override
497  {
498  const double tol = 1.0e-8;
499  const int maxiter = 1000;
500  //check orientation of surface
501  std::vector< unsigned int> orientation( 3 );
502 
503  double dummy = this->GetPoint( 0 ).X();
504 
505  if ( std::abs( this->GetPoint( 1 ).X() - dummy ) <= tol && std::abs( this->GetPoint( 2 ).X() - dummy ) <= tol && std::abs( this->GetPoint( 3 ).X() - dummy ) <= tol )
506  orientation[0] = 0;
507 
508  dummy = this->GetPoint( 0 ).Y();
509 
510  if ( std::abs( this->GetPoint( 1 ).Y() - dummy ) <= tol && std::abs( this->GetPoint( 2 ).Y() - dummy ) <= tol && std::abs( this->GetPoint( 3 ).Y() - dummy ) <= tol )
511  orientation[0] = 1;
512 
513  dummy = this->GetPoint( 0 ).Z();
514 
515  if ( std::abs( this->GetPoint( 1 ).Z() - dummy ) <= tol && std::abs( this->GetPoint( 2 ).Z() - dummy ) <= tol && std::abs( this->GetPoint( 3 ).Z() - dummy ) <= tol )
516  orientation[0] = 2;
517 
518  switch ( orientation[0] )
519  {
520  case 0:
521  orientation[0] = 1;
522  orientation[1] = 2;
523  orientation[2] = 0;
524  break;
525  case 1:
526  orientation[0] = 0;
527  orientation[1] = 2;
528  orientation[2] = 1;
529  break;
530  case 2:
531  orientation[0] = 0;
532  orientation[1] = 1;
533  orientation[2] = 2;
534  break;
535  default:
536  orientation[0] = 0;
537  orientation[1] = 1;
538  orientation[2] = 2;
539  }
540 
541  Matrix J = ZeroMatrix( 2, 2 );
542 
543  Matrix invJ = ZeroMatrix( 2, 2 );
544 
545  if ( rResult.size() != 2 )
546  rResult.resize( 2, false );
547 
548  //starting with xi = 0
549  rResult = ZeroVector( 3 );
550 
551  Vector DeltaXi = ZeroVector( 3 );
552 
553  CoordinatesArrayType CurrentGlobalCoords( ZeroVector( 3 ) );
554 
555  //Newton iteration:
556  for ( int k = 0; k < maxiter; k++ )
557  {
558  CurrentGlobalCoords = ZeroVector( 3 );
559  this->GlobalCoordinates( CurrentGlobalCoords, rResult );
560  noalias( CurrentGlobalCoords ) = rPoint - CurrentGlobalCoords;
561 
562  //Caluclate Inverse of Jacobian
563  noalias(J) = ZeroMatrix( 2, 2 );
564 
565  //derivatives of shape functions
566  Matrix shape_functions_gradients;
567  shape_functions_gradients = ShapeFunctionsLocalGradients(
568  shape_functions_gradients, rResult );
569 
570  //Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
571  //loop over all nodes
572  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
573  {
574  Point dummyPoint = this->GetPoint( i );
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 ) );
579  }
580 
581  //deteminant of Jacobian
582  double det_j = J( 0, 0 ) * J( 1, 1 ) - J( 0, 1 ) * J( 1, 0 );
583 
584  //filling matrix
585  invJ( 0, 0 ) = ( J( 1, 1 ) ) / ( det_j );
586 
587  invJ( 1, 0 ) = -( J( 1, 0 ) ) / ( det_j );
588 
589  invJ( 0, 1 ) = -( J( 0, 1 ) ) / ( det_j );
590 
591  invJ( 1, 1 ) = ( J( 0, 0 ) ) / ( det_j );
592 
593  DeltaXi( 0 ) = invJ( 0, 0 ) * CurrentGlobalCoords( orientation[0] ) + invJ( 0, 1 ) * CurrentGlobalCoords( orientation[1] );
594 
595  DeltaXi( 1 ) = invJ( 1, 0 ) * CurrentGlobalCoords( orientation[0] ) + invJ( 1, 1 ) * CurrentGlobalCoords( orientation[1] );
596 
597  noalias( rResult ) += DeltaXi;
598 
599  if ( MathUtils<double>::Norm3( DeltaXi ) > 30 )
600  {
601  break;
602  }
603 
604  if ( MathUtils<double>::Norm3( DeltaXi ) < tol )
605  {
606  if ( !( std::abs( CurrentGlobalCoords( orientation[2] ) ) <= tol ) )
607  rResult( 0 ) = 2.0;
608 
609  break;
610  }
611  }
612 
613  return( rResult );
614  }
615 
616 
639  IntegrationMethod ThisMethod ) const override
640  {
641  //getting derivatives of shape functions
642  ShapeFunctionsGradientsType shape_functions_gradients =
643  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
644  //getting values of shape functions
645  Matrix shape_functions_values =
646  CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
647  //workaround by riccardo...
648 
649  if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
650  {
651  // KLUDGE: While there is a bug in ublas
652  // vector resize, I have to put this beside resizing!!
653  JacobiansType temp( this->IntegrationPointsNumber( ThisMethod ) );
654  rResult.swap( temp );
655  }
656 
657  //loop over all integration points
658  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
659  {
660  //defining single jacobian matrix
661  Matrix jacobian = ZeroMatrix( 3, 2 );
662  //loop over all nodes
663 
664  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
665  {
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 ) );
672  }
673 
674  rResult[pnt] = jacobian;
675  }//end of loop over all integration points
676 
677  return rResult;
678  }
679 
698  IntegrationMethod ThisMethod,
699  Matrix & DeltaPosition ) const override
700  {
701  //getting derivatives of shape functions
702  ShapeFunctionsGradientsType shape_functions_gradients =
703  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
704  //getting values of shape functions
705  Matrix shape_functions_values =
706  CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
707  //workaround by riccardo...
708 
709  if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
710  {
711  // KLUDGE: While there is a bug in ublas
712  // vector resize, I have to put this beside resizing!!
713  JacobiansType temp( this->IntegrationPointsNumber( ThisMethod ) );
714  rResult.swap( temp );
715  }
716 
717  //loop over all integration points
718  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
719  {
720  //defining single jacobian matrix
721  Matrix jacobian = ZeroMatrix( 3, 2 );
722  //loop over all nodes
723 
724  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
725  {
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 ) );
732  }
733 
734  rResult[pnt] = jacobian;
735  }//end of loop over all integration points
736 
737  return rResult;
738  }
739 
761  Matrix& Jacobian( Matrix& rResult,
762  IndexType IntegrationPointIndex,
763  IntegrationMethod ThisMethod ) const override
764  {
765  // Setting up size of jacobian matrix
766  if (rResult.size1() != 3 || rResult.size2() != 2)
767  rResult.resize( 3, 2 , false);
768  rResult.clear();
769  //derivatives of shape functions
770  ShapeFunctionsGradientsType shape_functions_gradients =
771  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
772  Matrix ShapeFunctionsGradientInIntegrationPoint = shape_functions_gradients( IntegrationPointIndex );
773  //values of shape functions in integration points
774  DenseVector<double> ShapeFunctionsValuesInIntegrationPoint = ZeroVector( 9 );
775  /*vector<double>*/
776  ShapeFunctionsValuesInIntegrationPoint = row(
777  CalculateShapeFunctionsIntegrationPointsValues( ThisMethod ), IntegrationPointIndex );
778 
779  //Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
780  //loop over all nodes
781 
782  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
783  {
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 ) );
790  }
791 
792  return rResult;
793  }
794 
809  Matrix& Jacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
810  {
811  // Setting up size of jacobian matrix
812  if (rResult.size1() != 3 || rResult.size2() != 2)
813  rResult.resize( 3, 2 , false);
814  rResult.clear();
815  //derivatives of shape functions
816  Matrix shape_functions_gradients;
817  shape_functions_gradients = ShapeFunctionsLocalGradients(
818  shape_functions_gradients, rPoint );
819  //Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
820 
821  //loop over all nodes
822  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
823  {
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 ) );
831  }
832 
833  return rResult;
834  }
835 
839 
851  SizeType EdgesNumber() const override
852  {
853  return 4;
854  }
855 
865  {
867  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 0 ), this->pGetPoint( 1 ), this->pGetPoint( 4 ) ) );
868  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 1 ), this->pGetPoint( 2 ), this->pGetPoint( 5 ) ) );
869  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 2 ), this->pGetPoint( 3 ), this->pGetPoint( 6 ) ) );
870  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 3 ), this->pGetPoint( 0 ), this->pGetPoint( 7 ) ) );
871  return edges;
872  }
873 
890  double ShapeFunctionValue( IndexType ShapeFunctionIndex,
891  const CoordinatesArrayType& rPoint ) const override
892  {
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];
899 
900  switch ( ShapeFunctionIndex )
901  {
902  case 0:
903  return( fx1*fy1 );
904  case 1:
905  return( fx2*fy1 );
906  case 2:
907  return( fx2*fy2 );
908  case 3:
909  return( fx1*fy2 );
910  case 4:
911  return( fx3*fy1 );
912  case 5:
913  return( fx2*fy3 );
914  case 6:
915  return( fx3*fy2 );
916  case 7:
917  return( fx1*fy3 );
918  case 8:
919  return( fx3*fy3 );
920  default:
921  KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
922  }
923 
924  return 0;
925  }
926 
938  Vector& ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType& rCoordinates) const override
939  {
940  if(rResult.size() != 9) rResult.resize(9,false);
941 
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];
948 
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 ;
958 
959  return rResult;
960  }
961 
963 
975  std::string Info() const override
976  {
977  return "2 dimensional quadrilateral with nine nodes in 3D space";
978  }
979 
986  void PrintInfo( std::ostream& rOStream ) const override
987  {
988  rOStream << "2 dimensional quadrilateral with nine nodes in 3D space";
989  }
990 
1000  void PrintData( std::ostream& rOStream ) const override
1001  {
1002  // Base Geometry class PrintData call
1003  BaseType::PrintData( rOStream );
1004  std::cout << std::endl;
1005 
1006  // If the geometry has valid points, calculate and output its data
1007  if (this->AllPointsAreValid()) {
1008  Matrix jacobian;
1009  this->Jacobian( jacobian, PointType() );
1010  rOStream << " Jacobian in the origin\t : " << jacobian;
1011  }
1012  }
1013 
1019  IntegrationMethod ThisMethod )
1020  {
1021  ShapeFunctionsGradientsType localGradients
1022  = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1023  const int integration_points_number
1024  = msGeometryData.IntegrationPointsNumber( ThisMethod );
1025  ShapeFunctionsGradientsType Result( integration_points_number );
1026 
1027  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1028  {
1029  Result[pnt] = localGradients[pnt];
1030  }
1031 
1032  return Result;
1033  }
1034 
1040  {
1041  IntegrationMethod ThisMethod = msGeometryData.DefaultIntegrationMethod();
1042  ShapeFunctionsGradientsType localGradients
1043  = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
1044  const int integration_points_number
1045  = msGeometryData.IntegrationPointsNumber( ThisMethod );
1046  ShapeFunctionsGradientsType Result( integration_points_number );
1047 
1048  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1049  {
1050  Result[pnt] = localGradients[pnt];
1051  }
1052 
1053  return Result;
1054  }
1055 
1065  const CoordinatesArrayType& rPoint ) const override
1066  {
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];
1073 
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];
1080 
1081  rResult.resize( 9, 2, false );
1082  noalias( rResult ) = ZeroMatrix( 9, 2 );
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;
1101 
1102  return rResult;
1103  }
1104 
1110  Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
1111  {
1112  rResult.resize( 9, 2, false );
1113  noalias( rResult ) = ZeroMatrix( 9, 2 );
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;
1132  return rResult;
1133  }
1134 
1143  virtual Matrix& ShapeFunctionsGradients( Matrix& rResult, PointType& rPoint )
1144  {
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();
1151 
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();
1158 
1159  rResult.resize( 9, 2, false );
1160  noalias( rResult ) = ZeroMatrix( 9, 2 );
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;
1179 
1180  return rResult;
1181  }
1182 
1190  {
1191  if ( rResult.size() != this->PointsNumber() )
1192  {
1193  // KLUDGE: While there is a bug in ublas vector resize, I have to put this beside resizing!!
1195  rResult.swap( temp );
1196  }
1197 
1198  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
1199  {
1200  rResult[i].resize( 2, 2, false );
1201  noalias( rResult[i] ) = ZeroMatrix( 2, 2 );
1202  }
1203 
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];
1210 
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];
1217 
1218  double hx1 = 1.0;
1219  double hx2 = 1.0;
1220  double hx3 = -2.0;
1221  double hy1 = 1.0;
1222  double hy2 = 1.0;
1223  double hy3 = -2.0;
1224 
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;
1229 
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;
1234 
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;
1239 
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;
1244 
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;
1249 
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;
1254 
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;
1259 
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;
1264 
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;
1269 
1270  return rResult;
1271  }
1272 
1273 
1281  {
1282 
1283  if ( rResult.size() != this->PointsNumber() )
1284  {
1285  // KLUDGE: While there is a bug in
1286  // ublas vector resize, I have to put this beside resizing!!
1287 // ShapeFunctionsGradientsType
1289  rResult.swap( temp );
1290  }
1291 
1292  for ( IndexType i = 0; i < rResult.size(); i++ )
1293  {
1295  rResult[i].swap( temp );
1296  }
1297 
1298  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
1299  {
1300  for ( unsigned int j = 0; j < 2; j++ )
1301  {
1302  rResult[i][j].resize( 2, 2, false );
1303  noalias( rResult[i][j] ) = ZeroMatrix( 2, 2 );
1304  }
1305  }
1306 
1307 // double fx1 = 0.5*(rPoint[0]-1)*rPoint[0];
1308 // double fx2 = 0.5*(rPoint[0]+1)*rPoint[0];
1309 // double fx3 = 1-rPoint[0]*rPoint[0];
1310 // double fy1 = 0.5*(rPoint[1]-1)*rPoint[1];
1311 // double fy2 = 0.5*(rPoint[1]+1)*rPoint[1];
1312 // double fy3 = 1-rPoint[1]*rPoint[1];
1313 
1314  double gx1 = 0.5 * ( 2 * rPoint[0] - 1 );
1315 
1316  double gx2 = 0.5 * ( 2 * rPoint[0] + 1 );
1317 
1318  double gx3 = -2.0 * rPoint[0];
1319 
1320  double gy1 = 0.5 * ( 2 * rPoint[1] - 1 );
1321 
1322  double gy2 = 0.5 * ( 2 * rPoint[1] + 1 );
1323 
1324  double gy3 = -2.0 * rPoint[1];
1325 
1326  double hx1 = 1.0;
1327 
1328  double hx2 = 1.0;
1329 
1330  double hx3 = -2.0;
1331 
1332  double hy1 = 1.0;
1333 
1334  double hy2 = 1.0;
1335 
1336  double hy3 = -2.0;
1337 
1338  rResult[0][0]( 0, 0 ) = 0.0;
1339 
1340  rResult[0][0]( 0, 1 ) = hx1 * gy1;
1341 
1342  rResult[0][0]( 1, 0 ) = hx1 * gy1;
1343 
1344  rResult[0][0]( 1, 1 ) = gx1 * hy1;
1345 
1346  rResult[0][1]( 0, 0 ) = hx1 * gy1;
1347 
1348  rResult[0][1]( 0, 1 ) = gx1 * hy1;
1349 
1350  rResult[0][1]( 1, 0 ) = gx1 * hy1;
1351 
1352  rResult[0][1]( 1, 1 ) = 0.0;
1353 
1354  rResult[1][0]( 0, 0 ) = 0.0;
1355 
1356  rResult[1][0]( 0, 1 ) = hx2 * gy1;
1357 
1358  rResult[1][0]( 1, 0 ) = hx2 * gy1;
1359 
1360  rResult[1][0]( 1, 1 ) = gx2 * hy1;
1361 
1362  rResult[1][1]( 0, 0 ) = hx2 * gy1;
1363 
1364  rResult[1][1]( 0, 1 ) = gx2 * hy1;
1365 
1366  rResult[1][1]( 1, 0 ) = gx2 * hy1;
1367 
1368  rResult[1][1]( 1, 1 ) = 0.0;
1369 
1370  rResult[2][0]( 0, 0 ) = 0.0;
1371 
1372  rResult[2][0]( 0, 1 ) = hx2 * gy2;
1373 
1374  rResult[2][0]( 1, 0 ) = hx2 * gy2;
1375 
1376  rResult[2][0]( 1, 1 ) = gx2 * hy2;
1377 
1378  rResult[2][1]( 0, 0 ) = hx2 * gy2;
1379 
1380  rResult[2][1]( 0, 1 ) = gx2 * hy2;
1381 
1382  rResult[2][1]( 1, 0 ) = gx2 * hy2;
1383 
1384  rResult[2][1]( 1, 1 ) = 0.0;
1385 
1386  rResult[3][0]( 0, 0 ) = 0.0;
1387 
1388  rResult[3][0]( 0, 1 ) = hx1 * gy2;
1389 
1390  rResult[3][0]( 1, 0 ) = hx1 * gy2;
1391 
1392  rResult[3][0]( 1, 1 ) = gx1 * hy2;
1393 
1394  rResult[3][1]( 0, 0 ) = hx1 * gy2;
1395 
1396  rResult[3][1]( 0, 1 ) = gx1 * hy2;
1397 
1398  rResult[3][1]( 1, 0 ) = gx1 * hy2;
1399 
1400  rResult[3][1]( 1, 1 ) = 0.0;
1401 
1402  rResult[4][0]( 0, 0 ) = 0.0;
1403 
1404  rResult[4][0]( 0, 1 ) = hx3 * gy1;
1405 
1406  rResult[4][0]( 1, 0 ) = hx3 * gy1;
1407 
1408  rResult[4][0]( 1, 1 ) = gx3 * hy1;
1409 
1410  rResult[4][1]( 0, 0 ) = hx3 * gy1;
1411 
1412  rResult[4][1]( 0, 1 ) = gx3 * hy1;
1413 
1414  rResult[4][1]( 1, 0 ) = gx3 * hy1;
1415 
1416  rResult[4][1]( 1, 1 ) = 0.0;
1417 
1418  rResult[5][0]( 0, 0 ) = 0.0;
1419 
1420  rResult[5][0]( 0, 1 ) = hx2 * gy3;
1421 
1422  rResult[5][0]( 1, 0 ) = hx2 * gy3;
1423 
1424  rResult[5][0]( 1, 1 ) = gx2 * hy3;
1425 
1426  rResult[5][1]( 0, 0 ) = hx2 * gy3;
1427 
1428  rResult[5][1]( 0, 1 ) = gx2 * hy3;
1429 
1430  rResult[5][1]( 1, 0 ) = gx2 * hy3;
1431 
1432  rResult[5][1]( 1, 1 ) = 0.0;
1433 
1434  rResult[6][0]( 0, 0 ) = 0.0;
1435 
1436  rResult[6][0]( 0, 1 ) = hx3 * gy2;
1437 
1438  rResult[6][0]( 1, 0 ) = hx3 * gy2;
1439 
1440  rResult[6][0]( 1, 1 ) = gx3 * hy2;
1441 
1442  rResult[6][1]( 0, 0 ) = hx3 * gy2;
1443 
1444  rResult[6][1]( 0, 1 ) = gx3 * hy2;
1445 
1446  rResult[6][1]( 1, 0 ) = gx3 * hy2;
1447 
1448  rResult[6][1]( 1, 1 ) = 0.0;
1449 
1450  rResult[7][0]( 0, 0 ) = 0.0;
1451 
1452  rResult[7][0]( 0, 1 ) = hx1 * gy3;
1453 
1454  rResult[7][0]( 1, 0 ) = hx1 * gy3;
1455 
1456  rResult[7][0]( 1, 1 ) = gx1 * hy3;
1457 
1458  rResult[7][1]( 0, 0 ) = hx1 * gy3;
1459 
1460  rResult[7][1]( 0, 1 ) = gx1 * hy3;
1461 
1462  rResult[7][1]( 1, 0 ) = gx1 * hy3;
1463 
1464  rResult[7][1]( 1, 1 ) = 0.0;
1465 
1466  rResult[8][0]( 0, 0 ) = 0.0;
1467 
1468  rResult[8][0]( 0, 1 ) = hx3 * gy3;
1469 
1470  rResult[8][0]( 1, 0 ) = hx3 * gy3;
1471 
1472  rResult[8][0]( 1, 1 ) = gx3 * hy3;
1473 
1474  rResult[8][1]( 0, 0 ) = hx3 * gy3;
1475 
1476  rResult[8][1]( 0, 1 ) = gx3 * hy3;
1477 
1478  rResult[8][1]( 1, 0 ) = gx3 * hy3;
1479 
1480  rResult[8][1]( 1, 1 ) = 0.0;
1481 
1482 
1483  return rResult;
1484  }
1485 
1486 protected:
1487 
1488 
1494 private:
1495 
1498 
1499  static const GeometryData msGeometryData;
1500 
1501  static const GeometryDimension msGeometryDimension;
1502 
1503 
1507 
1508  friend class Serializer;
1509 
1510  void save( Serializer& rSerializer ) const override
1511  {
1513  }
1514 
1515  void load( Serializer& rSerializer ) override
1516  {
1518  }
1519 
1520  Quadrilateral3D9(): BaseType( PointsArrayType(), &msGeometryData ) {}
1521 
1540  static Matrix CalculateShapeFunctionsIntegrationPointsValues(
1541  typename BaseType::IntegrationMethod ThisMethod )
1542  {
1543  IntegrationPointsContainerType all_integration_points = AllIntegrationPoints();
1544  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1545  //number of integration points
1546  const int integration_points_number = integration_points.size();
1547  //number of nodes in current geometry
1548  const int points_number = 9;
1549  //setting up return matrix
1550  Matrix shape_function_values( integration_points_number, points_number );
1551  //loop over all integration points
1552 
1553  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1554  {
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();
1561 
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 );
1571  }
1572 
1573  return shape_function_values;
1574  }
1575 
1588  static ShapeFunctionsGradientsType CalculateShapeFunctionsIntegrationPointsLocalGradients(
1589  typename BaseType::IntegrationMethod ThisMethod )
1590  {
1591  IntegrationPointsContainerType all_integration_points = AllIntegrationPoints();
1592  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1593  //number of integration points
1594  const int integration_points_number = integration_points.size();
1595  ShapeFunctionsGradientsType d_shape_f_values( integration_points_number );
1596  //initialising container
1597  //std::fill(d_shape_f_values.begin(), d_shape_f_values.end(), Matrix(4,2));
1598  //loop over all integration points
1599 
1600  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1601  {
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();
1608 
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();
1615 
1616  Matrix result( 9, 2 );
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;
1635 
1636  d_shape_f_values[pnt] = result;
1637  }
1638 
1639  return d_shape_f_values;
1640  }
1641 
1645  static const IntegrationPointsContainerType AllIntegrationPoints()
1646  {
1647  IntegrationPointsContainerType integration_points =
1648  {
1649  {
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()
1660  }
1661  };
1662  return integration_points;
1663  }
1664 
1668  static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
1669  {
1670  ShapeFunctionsValuesContainerType shape_functions_values =
1671  {
1672  {
1673  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1675  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1677  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1679  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1681  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1683  }
1684  };
1685  return shape_functions_values;
1686  }
1687 
1691  static const ShapeFunctionsLocalGradientsContainerType AllShapeFunctionsLocalGradients()
1692  {
1693  ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
1694  {
1695  {
1696  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1698  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1700  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1702  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1704  Quadrilateral3D9<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients
1706  }
1707  };
1708  return shape_functions_local_gradients;
1709  }
1710 
1715  template<class TOtherPointType> friend class Quadrilateral3D9;
1716 
1720 }; // Class Quadrilateral3D9
1721 
1729 template< class TPointType > inline std::istream& operator >> (
1730  std::istream& rIStream,
1732 
1736 template< class TPointType > inline std::ostream& operator << (
1737  std::ostream& rOStream,
1738  const Quadrilateral3D9<TPointType>& rThis )
1739 {
1740  rThis.PrintInfo( rOStream );
1741  rOStream << std::endl;
1742  rThis.PrintData( rOStream );
1743  return rOStream;
1744 }
1745 
1746 template<class TPointType>
1747 const GeometryData Quadrilateral3D9<TPointType>::msGeometryData(
1750  Quadrilateral3D9<TPointType>::AllIntegrationPoints(),
1751  Quadrilateral3D9<TPointType>::AllShapeFunctionsValues(),
1752  AllShapeFunctionsLocalGradients()
1753 );
1754 
1755 template<class TPointType>
1756 const GeometryDimension Quadrilateral3D9<TPointType>::msGeometryDimension(3, 2);
1757 
1758 } // 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
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