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_interface_3d_4.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: Ignasi de Pouplana
11 //
12 
13 #if !defined(KRATOS_QUADRILATERAL_INTERFACE_3D_4_H_INCLUDED )
14 #define KRATOS_QUADRILATERAL_INTERFACE_3D_4_H_INCLUDED
15 
16 // System includes
17 
18 // External includes
19 
20 // Project includes
21 #include "geometries/line_3d_2.h"
23 
24 
25 namespace Kratos
26 {
29 
33 
37 
41 
45 
52 template<class TPointType> class QuadrilateralInterface3D4
53  : public Geometry<TPointType>
54 {
55 public:
59 
64 
69 
74 
79 
85 
89  typedef TPointType PointType;
90 
97  typedef typename BaseType::IndexType IndexType;
98 
104  typedef typename BaseType::SizeType SizeType;
105 
111 
116 
122 
131 
138 
144 
150 
157 
164 
172 
180 
185 
186 
190 
191 // QuadrilateralInterface3D4( const PointType& FirstPoint,
192 // const PointType& SecondPoint,
193 // const PointType& ThirdPoint,
194 // const PointType& FourthPoint )
195 // : BaseType( PointsArrayType(), &msGeometryData )
196 // {
197 // array_1d< double , 3 > vx;
198 // vx.clear();
199 // vx = SecondPoint + ThirdPoint - FirstPoint - FourthPoint;
200 // vx *= 0.5;
201 //
202 // array_1d< double , 3 > vy;
203 // vy.clear();
204 // vy = ThirdPoint + FourthPoint - FirstPoint - SecondPoint;
205 // vy *= 0.5;
206 //
207 // double lx = MathUtils<double>::Norm3(vx);
208 // double ly = MathUtils<double>::Norm3(vy);
209 //
210 // if(lx > ly)
211 // {
212 // this->Points().push_back( typename PointType::Pointer( new PointType( FirstPoint ) ) );
213 // this->Points().push_back( typename PointType::Pointer( new PointType( SecondPoint ) ) );
214 // this->Points().push_back( typename PointType::Pointer( new PointType( ThirdPoint ) ) );
215 // this->Points().push_back( typename PointType::Pointer( new PointType( FourthPoint ) ) );
216 // }
217 // else
218 // {
219 // this->Points().push_back( typename PointType::Pointer( new PointType( FourthPoint ) ) );
220 // this->Points().push_back( typename PointType::Pointer( new PointType( FirstPoint ) ) );
221 // this->Points().push_back( typename PointType::Pointer( new PointType( SecondPoint ) ) );
222 // this->Points().push_back( typename PointType::Pointer( new PointType( ThirdPoint ) ) );
223 // }
224 // }
225 
226  QuadrilateralInterface3D4( typename PointType::Pointer pFirstPoint,
227  typename PointType::Pointer pSecondPoint,
228  typename PointType::Pointer pThirdPoint,
229  typename PointType::Pointer pFourthPoint )
230  : BaseType( PointsArrayType(), &msGeometryData )
231  {
233  noalias(vx) = (*pSecondPoint + *pThirdPoint - *pFirstPoint - *pFourthPoint)*0.5;
234 
236  noalias(vy) = (*pThirdPoint + *pFourthPoint - *pFirstPoint - *pSecondPoint)*0.5;
237 
238  double lx = MathUtils<double>::Norm3(vx);
239  double ly = MathUtils<double>::Norm3(vy);
240 
241  if(ly < lx)
242  {
243  this->Points().push_back( pFirstPoint );
244  this->Points().push_back( pSecondPoint );
245  this->Points().push_back( pThirdPoint );
246  this->Points().push_back( pFourthPoint );
247  }
248  else
249  {
250  this->Points().push_back( pFourthPoint );
251  this->Points().push_back( pFirstPoint );
252  this->Points().push_back( pSecondPoint );
253  this->Points().push_back( pThirdPoint );
254  }
255  }
256 
258  : BaseType( ThisPoints, &msGeometryData )
259  {
260  if ( this->PointsNumber() != 4 )
261  KRATOS_ERROR << "Invalid points number. Expected 4, given " << this->PointsNumber() << std::endl;
262  }
263 
274  : BaseType( rOther )
275  {
276  }
277 
290  template<class TOtherPointType> QuadrilateralInterface3D4( QuadrilateralInterface3D4<TOtherPointType> const& rOther )
291  : BaseType( rOther )
292  {
293  }
294 
299 
301  {
303  }
304 
306  {
308  }
309 
313 
326  {
327  BaseType::operator=( rOther );
328  return *this;
329  }
330 
342  template<class TOtherPointType>
344  {
345  BaseType::operator=( rOther );
346  return *this;
347  }
348 
352 
353  typename BaseType::Pointer Create( PointsArrayType const& ThisPoints ) const override
354  {
355  return typename BaseType::Pointer( new QuadrilateralInterface3D4( ThisPoints ) );
356  }
357 
363  Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
364  {
365  noalias( rResult ) = ZeroMatrix( 4, 2 );
366  rResult( 0, 0 ) = -1.0;
367  rResult( 0, 1 ) = -1.0;
368  rResult( 1, 0 ) = 1.0;
369  rResult( 1, 1 ) = -1.0;
370  rResult( 2, 0 ) = 1.0;
371  rResult( 2, 1 ) = 1.0;
372  rResult( 3, 0 ) = -1.0;
373  rResult( 3, 1 ) = 1.0;
374  return rResult;
375  }
376 
380 
400  double Length() const override
401  {
404 
405  array_1d<double, 3> vx( p1 - p0 );
406 
407  return MathUtils<double>::Norm3(vx);
408  }
409 
425  double Area() const override
426  {
427  return Length();
428  }
429 
444  double DomainSize() const override
445  {
446  return Area();
447  }
448 
456  double Volume() const override
457  {
458  KRATOS_WARNING("QuadrilateralInterface3D4") << "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;
459  return Area();
460  // TODO: Replace in June 2023
461  // KRATOS_ERROR << "QuadrilateralInterface3D4:: Method not well defined. Replace with DomainSize() instead." << std::endl;
462  // return 0.0;
463  }
464 
473  bool IsInside(
474  const CoordinatesArrayType& rPoint,
475  CoordinatesArrayType& rResult,
476  const double Tolerance = std::numeric_limits<double>::epsilon()
477  ) const override
478  {
479  PointLocalCoordinates( rResult, rPoint );
480 
481  if ( (rResult[0] >= (-1.0 - Tolerance)) && (rResult[0] <= (1.0 + Tolerance)) )
482  {
483  if ( (rResult[1] >= (-1.0 - Tolerance)) && (rResult[1] <= (1.0 + Tolerance)) )
484  {
485  return true;
486  }
487  }
488 
489  return false;
490  }
491 
499  CoordinatesArrayType& rResult,
500  const CoordinatesArrayType& rPoint
501  ) const override
502  {
505  for(unsigned int i=0; i<this->size();i++) {
506  X(0,i ) = this->GetPoint( i ).X();
507  X(1,i ) = this->GetPoint( i ).Y();
508  X(2,i ) = this->GetPoint( i ).Z();
509  }
510 
511  double tol = 1.0e-8;
512  int maxiter = 1000;
513 
514  Matrix J = ZeroMatrix( 2, 2 );
515  Matrix invJ = ZeroMatrix( 2, 2 );
516 
517  //starting with xi = 0
518  rResult = ZeroVector( 3 );
519  Vector DeltaXi = ZeroVector( 2 );
520  array_1d<double,3> CurrentGlobalCoords;
521 
522 
523  //Newton iteration:
524  for ( int k = 0; k < maxiter; k++ )
525  {
526  noalias(CurrentGlobalCoords) = ZeroVector( 3 );
527  this->GlobalCoordinates( CurrentGlobalCoords, rResult );
528 
529  noalias( CurrentGlobalCoords ) = rPoint - CurrentGlobalCoords;
530 
531 
532  //derivatives of shape functions
533  Matrix shape_functions_gradients;
534  shape_functions_gradients = ShapeFunctionsLocalGradients(shape_functions_gradients, rResult );
535  noalias(DN) = prod(X,shape_functions_gradients);
536 
537  noalias(J) = prod(trans(DN),DN);
538  Vector res = prod(trans(DN),CurrentGlobalCoords);
539 
540  //deteminant of Jacobian
541  const double det_j = J( 0, 0 ) * J( 1, 1 ) - J( 0, 1 ) * J( 1, 0 );
542 
543  //filling matrix
544  invJ( 0, 0 ) = ( J( 1, 1 ) ) / ( det_j );
545  invJ( 1, 0 ) = -( J( 1, 0 ) ) / ( det_j );
546  invJ( 0, 1 ) = -( J( 0, 1 ) ) / ( det_j );
547  invJ( 1, 1 ) = ( J( 0, 0 ) ) / ( det_j );
548 
549 
550  DeltaXi( 0 ) = invJ( 0, 0 ) * res[0] + invJ( 0, 1 ) * res[1];
551  DeltaXi( 1 ) = invJ( 1, 0 ) * res[0] + invJ( 1, 1 ) * res[1];
552 
553  rResult[0] += DeltaXi[0];
554  rResult[1] += DeltaXi[1];
555  rResult[2] = 0.0;
556 
557  if ( norm_2( DeltaXi ) > 300 )
558  {
559  res[0] = 0.0;
560  res[1] = 0.0;
561  std::cout << "detJ =" << det_j << "DeltaX = " << DeltaXi << " stopping calculation and assigning the baricenter" << std::endl;
562  break;
563  //KRATOS_ERROR << "Computation of local coordinates failed at iteration" << k << std::endl;
564  }
565 
566  if ( norm_2( DeltaXi ) < tol )
567  {
568  break;
569  }
570  }
571 
572  return( rResult );
573  }
574 
594  IntegrationMethod ThisMethod,
595  Matrix & DeltaPosition ) const override
596  {
599 
602  dp0[0] = (DeltaPosition(0,0) + DeltaPosition(3,0))*0.5;
603  dp0[1] = (DeltaPosition(0,1) + DeltaPosition(3,1))*0.5;
604  dp0[2] = (DeltaPosition(0,2) + DeltaPosition(3,2))*0.5;
605  dp1[0] = (DeltaPosition(1,0) + DeltaPosition(2,0))*0.5;
606  dp1[1] = (DeltaPosition(1,1) + DeltaPosition(2,1))*0.5;
607  dp1[2] = (DeltaPosition(1,2) + DeltaPosition(2,2))*0.5;
608 
609  Matrix jacobian( 3, 1 );
610  jacobian( 0, 0 ) = ( p1[0]-dp1[0] - (p0[0]-dp0[0]) ) * 0.5; //on the Gauss points (J is constant at each element)
611  jacobian( 1, 0 ) = ( p1[1]-dp1[1] - (p0[1]-dp0[1]) ) * 0.5;
612  jacobian( 2, 0 ) = ( p1[2]-dp1[2] - (p0[2]-dp0[2]) ) * 0.5;
613 
614  if ( rResult.size() != BaseType::IntegrationPointsNumber( ThisMethod ) )
615  {
616  // KLUDGE: While there is a bug in ublas vector resize, I have to put this beside resizing!!
618  rResult.swap( temp );
619  }
620 
621  std::fill( rResult.begin(), rResult.end(), jacobian );
622 
623  return rResult;
624  }
646  Matrix& Jacobian( Matrix& rResult,
647  IndexType IntegrationPointIndex,
648  IntegrationMethod ThisMethod ) const override
649  {
652 
653  if(rResult.size1() != 3 || rResult.size2() != 1)
654  rResult.resize(3, 1, false);
655 
656  rResult( 0, 0 ) = ( p1[0] - p0[0] ) * 0.5;
657  rResult( 1, 0 ) = ( p1[1] - p0[1] ) * 0.5;
658  rResult( 2, 0 ) = ( p1[2] - p0[2] ) * 0.5;
659 
660  return rResult;
661  }
662 
678  Matrix& Jacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
679  {
682 
683  if(rResult.size1() != 3 || rResult.size2() != 1)
684  rResult.resize(3, 1, false);
685 
686  rResult( 0, 0 ) = ( p1[0] - p0[0] ) * 0.5;
687  rResult( 1, 0 ) = ( p1[1] - p0[1] ) * 0.5;
688  rResult( 2, 0 ) = ( p1[2] - p0[2] ) * 0.5;
689 
690  return rResult;
691  }
692 
709  IntegrationMethod ThisMethod ) const override
710  {
711  KRATOS_ERROR << "QuadrilateralInterface3D4::DeterminantOfJacobian" <<"Jacobian is not square" << std::endl;
712  return rResult;
713  }
714 
737  double DeterminantOfJacobian( IndexType IntegrationPointIndex,
738  IntegrationMethod ThisMethod ) const override
739  {
740  KRATOS_ERROR << "QuadrilateralInterface3D4::DeterminantOfJacobian" << "Jacobian is not square" << std::endl;
741  return 0.0;
742  }
743 
769  double DeterminantOfJacobian( const CoordinatesArrayType& rPoint ) const override
770  {
771  KRATOS_ERROR << "QuadrilateralInterface3D4::DeterminantOfJacobian" << "Jacobian is not square" << std::endl;
772  return 0.0;
773  }
774 
797  IntegrationMethod ThisMethod ) const override
798  {
799  KRATOS_ERROR << "QuadrilateralInterface3D4::DeterminantOfJacobian" << "Jacobian is not square" << std::endl;
800  return rResult;
801  }
802 
827  IndexType IntegrationPointIndex,
828  IntegrationMethod ThisMethod ) const override
829  {
830  KRATOS_ERROR << "QuadrilateralInterface3D4::DeterminantOfJacobian" << "Jacobian is not square" << std::endl;
831  return rResult;
832  }
833 
851  const CoordinatesArrayType& rPoint ) const override
852  {
853  KRATOS_ERROR << "QuadrilateralInterface3D4::DeterminantOfJacobian" << "Jacobian is not square" << std::endl;
854  return rResult;
855  }
856 
857 
861 
873  SizeType EdgesNumber() const override
874  {
875  return 4;
876  }
877 
887  {
889  typedef typename Geometry<TPointType>::Pointer EdgePointerType;
890  edges.push_back( EdgePointerType( new EdgeType( this->pGetPoint( 0 ), this->pGetPoint( 1 ) ) ) );
891  edges.push_back( EdgePointerType( new EdgeType( this->pGetPoint( 1 ), this->pGetPoint( 2 ) ) ) );
892  edges.push_back( EdgePointerType( new EdgeType( this->pGetPoint( 2 ), this->pGetPoint( 3 ) ) ) );
893  edges.push_back( EdgePointerType( new EdgeType( this->pGetPoint( 3 ), this->pGetPoint( 0 ) ) ) );
894  return edges;
895  }
896 
900 
913  double ShapeFunctionValue( IndexType ShapeFunctionIndex,
914  const CoordinatesArrayType& rPoint ) const override
915  {
916  switch ( ShapeFunctionIndex )
917  {
918  case 0:
919  return( 0.25*( 1.0 - rPoint[0] )*( 1.0 - rPoint[1] ) );
920  case 1:
921  return( 0.25*( 1.0 + rPoint[0] )*( 1.0 - rPoint[1] ) );
922  case 2:
923  return( 0.25*( 1.0 + rPoint[0] )*( 1.0 + rPoint[1] ) );
924  case 3:
925  return( 0.25*( 1.0 - rPoint[0] )*( 1.0 + rPoint[1] ) );
926  default:
927  KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
928  }
929 
930  return 0;
931  }
932 
944  Vector& ShapeFunctionsValues (Vector &rResult, const CoordinatesArrayType& rCoordinates) const override
945  {
946  if(rResult.size() != 4) rResult.resize(4,false);
947  rResult[0] = 0.25*( 1.0 - rCoordinates[0] )*( 1.0 - rCoordinates[1] );
948  rResult[1] = 0.25*( 1.0 + rCoordinates[0] )*( 1.0 - rCoordinates[1] );
949  rResult[2] = 0.25*( 1.0 + rCoordinates[0] )*( 1.0 + rCoordinates[1] );
950  rResult[3] = 0.25*( 1.0 - rCoordinates[0] )*( 1.0 + rCoordinates[1] );
951 
952  return rResult;
953  }
954 
973  IntegrationMethod ThisMethod ) const override
974  {
975  const unsigned int integration_points_number =
976  msGeometryData.IntegrationPointsNumber( ThisMethod );
977 
978  if ( integration_points_number == 0 )
979  KRATOS_ERROR << "This integration method is not supported" << *this << std::endl;
980 
981  //workaround by riccardo
982  if ( rResult.size() != integration_points_number )
983  {
984  // KLUDGE: While there is a bug in ublas
985  // vector resize, I have to put this beside resizing!!
986  ShapeFunctionsGradientsType temp( integration_points_number );
987  rResult.swap( temp );
988  }
989 
990  //calculating the local gradients
992  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
993 
994  //getting the inverse jacobian matrices
995  JacobiansType temp( integration_points_number );
996 
997  JacobiansType invJ = InverseOfJacobian( temp, ThisMethod );
998 
999  //loop over all integration points
1000  for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
1001  {
1002  rResult[pnt].resize( 4, 2, false );
1003 
1004  for ( int i = 0; i < 4; i++ )
1005  {
1006  for ( int j = 0; j < 2; j++ )
1007  {
1008  rResult[pnt]( i, j ) =
1009  ( locG[pnt]( i, 0 ) * invJ[pnt]( j, 0 ) )
1010  + ( locG[pnt]( i, 1 ) * invJ[pnt]( j, 1 ) );
1011  }
1012  }
1013  }//end of loop over integration points
1014  }
1015 
1019 
1027  std::string Info() const override
1028  {
1029  return "3 dimensional quadrilateral with four nodes in 3D space";
1030  }
1031 
1038  void PrintInfo( std::ostream& rOStream ) const override
1039  {
1040  rOStream << "3 dimensional quadrilateral with four nodes in 3D space";
1041  }
1042 
1057  void PrintData( std::ostream& rOStream ) const override
1058  {
1059  // Base Geometry class PrintData call
1060  BaseType::PrintData( rOStream );
1061  std::cout << std::endl;
1062 
1063  // If the geometry has valid points, calculate and output its data
1064  if (this->AllPointsAreValid()) {
1065  Matrix jacobian;
1066  this->Jacobian( jacobian, PointType() );
1067  rOStream << " Jacobian in the origin\t : " << jacobian;
1068  }
1069  }
1070 
1081  const CoordinatesArrayType& rPoint ) const override
1082  {
1083  rResult.resize( 4, 2, false );
1084  noalias( rResult ) = ZeroMatrix( 4, 2 );
1085  rResult( 0, 0 ) = -0.25 * ( 1.0 - rPoint[1] );
1086  rResult( 0, 1 ) = -0.25 * ( 1.0 - rPoint[0] );
1087  rResult( 1, 0 ) = 0.25 * ( 1.0 - rPoint[1] );
1088  rResult( 1, 1 ) = -0.25 * ( 1.0 + rPoint[0] );
1089  rResult( 2, 0 ) = 0.25 * ( 1.0 + rPoint[1] );
1090  rResult( 2, 1 ) = 0.25 * ( 1.0 + rPoint[0] );
1091  rResult( 3, 0 ) = -0.25 * ( 1.0 + rPoint[1] );
1092  rResult( 3, 1 ) = 0.25 * ( 1.0 - rPoint[0] );
1093  return rResult;
1094  }
1095 
1103  const CoordinatesArrayType& rPoint ) const override
1104  {
1105  if ( rResult.size() != this->PointsNumber() )
1106  {
1107  // KLUDGE: While there is a bug in
1108  // ublas vector resize, I have to put this beside resizing!!
1110  rResult.swap( temp );
1111  }
1112 
1113  rResult[0].resize( 2, 2, false );
1114 
1115  rResult[1].resize( 2, 2, false );
1116  rResult[2].resize( 2, 2, false );
1117  rResult[3].resize( 2, 2, false );
1118  rResult[0]( 0, 0 ) = 0.0;
1119  rResult[0]( 0, 1 ) = 0.25;
1120  rResult[0]( 1, 0 ) = 0.25;
1121  rResult[0]( 1, 1 ) = 0.0;
1122  rResult[1]( 0, 0 ) = 0.0;
1123  rResult[1]( 0, 1 ) = -0.25;
1124  rResult[1]( 1, 0 ) = -0.25;
1125  rResult[1]( 1, 1 ) = 0.0;
1126  rResult[2]( 0, 0 ) = 0.0;
1127  rResult[2]( 0, 1 ) = 0.25;
1128  rResult[2]( 1, 0 ) = 0.25;
1129  rResult[2]( 1, 1 ) = 0.0;
1130  rResult[3]( 0, 0 ) = 0.0;
1131  rResult[3]( 0, 1 ) = -0.25;
1132  rResult[3]( 1, 0 ) = -0.25;
1133  rResult[3]( 1, 1 ) = 0.0;
1134  return rResult;
1135  }
1136 
1144  const CoordinatesArrayType& rPoint ) const override
1145  {
1146  if ( rResult.size() != this->PointsNumber() )
1147  {
1148  // KLUDGE: While there is a bug in
1149  // ublas vector resize, I have to put this beside resizing!!
1150 // ShapeFunctionsGradientsType
1152  rResult.swap( temp );
1153  }
1154 
1155  for ( IndexType i = 0; i < rResult.size(); i++ )
1156  {
1158  rResult[i].swap( temp );
1159  }
1160 
1161  rResult[0][0].resize( 2, 2 , false);
1162 
1163  rResult[0][1].resize( 2, 2 , false);
1164  rResult[1][0].resize( 2, 2 , false);
1165  rResult[1][1].resize( 2, 2 , false);
1166  rResult[2][0].resize( 2, 2 , false);
1167  rResult[2][1].resize( 2, 2 , false);
1168  rResult[3][0].resize( 2, 2 , false);
1169  rResult[3][1].resize( 2, 2 , false);
1170 
1171  for ( int i = 0; i < 4; i++ )
1172  {
1173 
1174  rResult[i][0]( 0, 0 ) = 0.0;
1175  rResult[i][0]( 0, 1 ) = 0.0;
1176  rResult[i][0]( 1, 0 ) = 0.0;
1177  rResult[i][0]( 1, 1 ) = 0.0;
1178  rResult[i][1]( 0, 0 ) = 0.0;
1179  rResult[i][1]( 0, 1 ) = 0.0;
1180  rResult[i][1]( 1, 0 ) = 0.0;
1181  rResult[i][1]( 1, 1 ) = 0.0;
1182  }
1183 
1184  return rResult;
1185  }
1186 
1190 
1192 
1193 protected:
1198 private:
1201 
1202  static const GeometryData msGeometryData;
1203 
1204  static const GeometryDimension msGeometryDimension;
1205 
1209 
1210  friend class Serializer;
1211 
1212  void save( Serializer& rSerializer ) const override
1213  {
1215  }
1216 
1217  void load( Serializer& rSerializer ) override
1218  {
1220  }
1221 
1222  QuadrilateralInterface3D4(): BaseType( PointsArrayType(), &msGeometryData ) {}
1223 
1227 
1228 
1232 
1233 
1237 
1248  static Matrix CalculateShapeFunctionsIntegrationPointsValues(
1249  typename BaseType::IntegrationMethod ThisMethod )
1250  {
1251  IntegrationPointsContainerType all_integration_points =
1252  AllIntegrationPoints();
1253  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1254  //number of integration points
1255  const int integration_points_number = integration_points.size();
1256  //number of nodes in current geometry
1257  const int points_number = 4;
1258  //setting up return matrix
1259  Matrix shape_function_values( integration_points_number, points_number );
1260  //loop over all integration points
1261 
1262  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1263  {
1264  shape_function_values( pnt, 0 ) =
1265  0.25 * ( 1.0 - integration_points[pnt].X() )
1266  * ( 1.0 - integration_points[pnt].Y() );
1267  shape_function_values( pnt, 1 ) =
1268  0.25 * ( 1.0 + integration_points[pnt].X() )
1269  * ( 1.0 - integration_points[pnt].Y() );
1270  shape_function_values( pnt, 2 ) =
1271  0.25 * ( 1.0 + integration_points[pnt].X() )
1272  * ( 1.0 + integration_points[pnt].Y() );
1273  shape_function_values( pnt, 3 ) =
1274  0.25 * ( 1.0 - integration_points[pnt].X() )
1275  * ( 1.0 + integration_points[pnt].Y() );
1276  }
1277 
1278  return shape_function_values;
1279  }
1280 
1293  static ShapeFunctionsGradientsType CalculateShapeFunctionsIntegrationPointsLocalGradients(
1294  typename BaseType::IntegrationMethod ThisMethod )
1295  {
1296  IntegrationPointsContainerType all_integration_points =
1297  AllIntegrationPoints();
1298  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1299  //number of integration points
1300  const int integration_points_number = integration_points.size();
1301  ShapeFunctionsGradientsType d_shape_f_values( integration_points_number );
1302  //initialising container
1303  //std::fill(d_shape_f_values.begin(), d_shape_f_values.end(), Matrix(4,2));
1304  //loop over all integration points
1305 
1306  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1307  {
1308  Matrix result( 4, 2 );
1309  result( 0, 0 ) = -0.25 * ( 1.0 - integration_points[pnt].Y() );
1310  result( 0, 1 ) = -0.25 * ( 1.0 - integration_points[pnt].X() );
1311  result( 1, 0 ) = 0.25 * ( 1.0 - integration_points[pnt].Y() );
1312  result( 1, 1 ) = -0.25 * ( 1.0 + integration_points[pnt].X() );
1313  result( 2, 0 ) = 0.25 * ( 1.0 + integration_points[pnt].Y() );
1314  result( 2, 1 ) = 0.25 * ( 1.0 + integration_points[pnt].X() );
1315  result( 3, 0 ) = -0.25 * ( 1.0 + integration_points[pnt].Y() );
1316  result( 3, 1 ) = 0.25 * ( 1.0 - integration_points[pnt].X() );
1317  d_shape_f_values[pnt] = result;
1318  }
1319 
1320  return d_shape_f_values;
1321  }
1322 
1326  static const IntegrationPointsContainerType AllIntegrationPoints()
1327  {
1328  IntegrationPointsContainerType integration_points =
1329  {
1330  {
1331  Quadrature < QuadrilateralGaussLobattoIntegrationPoints1,
1332  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1333  Quadrature < QuadrilateralGaussLobattoIntegrationPoints2,
1334  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1337  }
1338  };
1339  return integration_points;
1340  }
1341 
1345  static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
1346  {
1347  ShapeFunctionsValuesContainerType shape_functions_values =
1348  {
1349  {
1350  QuadrilateralInterface3D4<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1352  QuadrilateralInterface3D4<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1354  Matrix(),
1355  Matrix()
1356  }
1357  };
1358  return shape_functions_values;
1359  }
1360 
1365  AllShapeFunctionsLocalGradients()
1366  {
1367  ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
1368  {
1369  {
1370  QuadrilateralInterface3D4<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
1371  QuadrilateralInterface3D4<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
1374  }
1375  };
1376  return shape_functions_local_gradients;
1377  }
1378 
1382 
1383 
1387 
1388 
1392 
1393  template<class TOtherPointType> friend class QuadrilateralInterface3D4;
1394 
1398 
1399 
1401 }; // Class Geometry
1402 
1406 
1407 
1411 
1414 template<class TPointType> inline std::istream& operator >> (
1415  std::istream& rIStream,
1420 template<class TPointType> inline std::ostream& operator << (
1421  std::ostream& rOStream,
1423 {
1424  rThis.PrintInfo( rOStream );
1425  rOStream << std::endl;
1426  rThis.PrintData( rOStream );
1427  return rOStream;
1428 }
1429 
1431 
1432 template<class TPointType> const
1433 GeometryData QuadrilateralInterface3D4<TPointType>::msGeometryData(
1436  QuadrilateralInterface3D4<TPointType>::AllIntegrationPoints(),
1437  QuadrilateralInterface3D4<TPointType>::AllShapeFunctionsValues(),
1438  AllShapeFunctionsLocalGradients()
1439 );
1440 
1441 template<class TPointType>
1443 
1444 }// namespace Kratos.
1445 
1446 #endif // KRATOS_QUADRILATERAL_INTERFACE_3D_4_H_INCLUDED defined
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
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
SizeType size() const
Definition: geometry.h:518
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
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
const ShapeFunctionsGradientsType & ShapeFunctionsLocalGradients() const
Definition: geometry.h:3539
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
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
iterator end()
Definition: amatrix_interface.h:243
iterator begin()
Definition: amatrix_interface.h:241
An two node 3D line geometry with linear shape functions.
Definition: line_3d_2.h:64
static double Norm3(const TVectorType &a)
Calculates the norm of vector "a" which is assumed to be of size 3.
Definition: math_utils.h:691
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
Definition: quadrilateral_interface_3d_4.h:54
CoordinatesArrayType & PointLocalCoordinates(CoordinatesArrayType &rResult, const CoordinatesArrayType &rPoint) const override
Returns the local coordinates of a given arbitrary point.
Definition: quadrilateral_interface_3d_4.h:498
BaseType::ShapeFunctionsSecondDerivativesType ShapeFunctionsSecondDerivativesType
Definition: quadrilateral_interface_3d_4.h:171
void ShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_3d_4.h:971
BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: quadrilateral_interface_3d_4.h:163
Geometry< TPointType > BaseType
Definition: quadrilateral_interface_3d_4.h:63
double DomainSize() const override
Definition: quadrilateral_interface_3d_4.h:444
double ShapeFunctionValue(IndexType ShapeFunctionIndex, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:913
BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: quadrilateral_interface_3d_4.h:149
double DeterminantOfJacobian(IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_3d_4.h:737
GeometryData::KratosGeometryFamily GetGeometryFamily() const override
Definition: quadrilateral_interface_3d_4.h:300
void PrintInfo(std::ostream &rOStream) const override
Definition: quadrilateral_interface_3d_4.h:1038
virtual ~QuadrilateralInterface3D4()
Definition: quadrilateral_interface_3d_4.h:298
BaseType::IndexType IndexType
Definition: quadrilateral_interface_3d_4.h:97
QuadrilateralInterface3D4 & operator=(const QuadrilateralInterface3D4 &rOther)
Definition: quadrilateral_interface_3d_4.h:325
SizeType EdgesNumber() const override
This method gives you number of all edges of this geometry.
Definition: quadrilateral_interface_3d_4.h:873
QuadrilateralInterface3D4(QuadrilateralInterface3D4< TOtherPointType > const &rOther)
Definition: quadrilateral_interface_3d_4.h:290
JacobiansType & Jacobian(JacobiansType &rResult, IntegrationMethod ThisMethod, Matrix &DeltaPosition) const override
Definition: quadrilateral_interface_3d_4.h:593
QuadrilateralInterface3D4(QuadrilateralInterface3D4 const &rOther)
Definition: quadrilateral_interface_3d_4.h:273
BaseType::IntegrationPointsArrayType IntegrationPointsArrayType
Definition: quadrilateral_interface_3d_4.h:130
QuadrilateralInterface3D4(const PointsArrayType &ThisPoints)
Definition: quadrilateral_interface_3d_4.h:257
GeometriesArrayType GenerateEdges() const override
This method gives you all edges of this geometry.
Definition: quadrilateral_interface_3d_4.h:886
Matrix & ShapeFunctionsLocalGradients(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:1080
BaseType::IntegrationPointsContainerType IntegrationPointsContainerType
Definition: quadrilateral_interface_3d_4.h:137
BaseType::JacobiansType JacobiansType
Definition: quadrilateral_interface_3d_4.h:156
Matrix & InverseOfJacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_3d_4.h:826
Matrix & InverseOfJacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:850
QuadrilateralInterface3D4 & operator=(QuadrilateralInterface3D4< TOtherPointType > const &rOther)
Definition: quadrilateral_interface_3d_4.h:343
Vector & DeterminantOfJacobian(Vector &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_3d_4.h:708
BaseType::IntegrationPointType IntegrationPointType
Definition: quadrilateral_interface_3d_4.h:121
double DeterminantOfJacobian(const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:769
double Volume() const override
This method calculates and returns the volume of this geometry.
Definition: quadrilateral_interface_3d_4.h:456
BaseType::NormalType NormalType
Definition: quadrilateral_interface_3d_4.h:184
Matrix & PointsLocalCoordinates(Matrix &rResult) const override
Definition: quadrilateral_interface_3d_4.h:363
friend class QuadrilateralInterface3D4
Definition: quadrilateral_interface_3d_4.h:1393
KRATOS_CLASS_POINTER_DEFINITION(QuadrilateralInterface3D4)
BaseType::CoordinatesArrayType CoordinatesArrayType
Definition: quadrilateral_interface_3d_4.h:115
JacobiansType & InverseOfJacobian(JacobiansType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_3d_4.h:796
TPointType PointType
Definition: quadrilateral_interface_3d_4.h:89
BaseType::Pointer Create(PointsArrayType const &ThisPoints) const override
Creates a new geometry pointer.
Definition: quadrilateral_interface_3d_4.h:353
GeometryData::KratosGeometryType GetGeometryType() const override
Definition: quadrilateral_interface_3d_4.h:305
std::string Info() const override
Definition: quadrilateral_interface_3d_4.h:1027
GeometryData::IntegrationMethod IntegrationMethod
Definition: quadrilateral_interface_3d_4.h:78
BaseType::ShapeFunctionsThirdDerivativesType ShapeFunctionsThirdDerivativesType
Definition: quadrilateral_interface_3d_4.h:179
Matrix & Jacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:678
BaseType::SizeType SizeType
Definition: quadrilateral_interface_3d_4.h:104
BaseType::GeometriesArrayType GeometriesArrayType
Definition: quadrilateral_interface_3d_4.h:84
bool IsInside(const CoordinatesArrayType &rPoint, CoordinatesArrayType &rResult, const double Tolerance=std::numeric_limits< double >::epsilon()) const override
Definition: quadrilateral_interface_3d_4.h:473
double Length() const override
Definition: quadrilateral_interface_3d_4.h:400
Line3D2< TPointType > EdgeType
Definition: quadrilateral_interface_3d_4.h:68
void PrintData(std::ostream &rOStream) const override
Definition: quadrilateral_interface_3d_4.h:1057
ShapeFunctionsThirdDerivativesType & ShapeFunctionsThirdDerivatives(ShapeFunctionsThirdDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:1143
Matrix & Jacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_3d_4.h:646
BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType
Definition: quadrilateral_interface_3d_4.h:143
Vector & ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType &rCoordinates) const override
Definition: quadrilateral_interface_3d_4.h:944
double Area() const override
Definition: quadrilateral_interface_3d_4.h:425
BaseType::PointsArrayType PointsArrayType
Definition: quadrilateral_interface_3d_4.h:110
ShapeFunctionsSecondDerivativesType & ShapeFunctionsSecondDerivatives(ShapeFunctionsSecondDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_3d_4.h:1102
QuadrilateralInterface3D4(typename PointType::Pointer pFirstPoint, typename PointType::Pointer pSecondPoint, typename PointType::Pointer pThirdPoint, typename PointType::Pointer pFourthPoint)
Definition: quadrilateral_interface_3d_4.h:226
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_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
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
const GeometryData QuadrilateralInterface3D4< TPointType >::msGeometryData & msGeometryDimension(), QuadrilateralInterface3D4< 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::TransposeMatrix< const T > trans(const T &TheMatrix)
Definition: amatrix_interface.h:486
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
res
Definition: generate_convection_diffusion_explicit_element.py:211
DN
Definition: generate_convection_diffusion_explicit_element.py:98
int tol
Definition: hinsberg_optimization.py:138
def load(f)
Definition: ode_solve.py:307
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
float temp
Definition: rotating_cone.py:85
J
Definition: sensitivityMatrix.py:58
integer i
Definition: TensorModule.f:17