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_2d_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_2D_4_H_INCLUDED )
14 #define KRATOS_QUADRILATERAL_INTERFACE_2D_4_H_INCLUDED
15 
16 // System includes
17 
18 // External includes
19 
20 // Project includes
21 #include "geometries/line_2d_2.h"
23 
24 
25 namespace Kratos
26 {
29 
33 
37 
41 
45 
52 template<class TPointType> class QuadrilateralInterface2D4
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 
189 
190 // QuadrilateralInterface2D4( const PointType& FirstPoint,
191 // const PointType& SecondPoint,
192 // const PointType& ThirdPoint,
193 // const PointType& FourthPoint )
194 // : BaseType( PointsArrayType(), &msGeometryData )
195 // {
196 // double p2_p1_X = (ThirdPoint->X() + FourthPoint->X() - FirstPoint->X() - SecondPoint->X())*0.5;
197 // double p2_p1_Y = (ThirdPoint->Y() + FourthPoint->Y() - FirstPoint->Y() - SecondPoint->Y())*0.5;
198 // double p4_p3_X = (ThirdPoint->X() + SecondPoint->X() - FirstPoint->X() - FourthPoint->X())*0.5;
199 // double p4_p3_Y = (ThirdPoint->Y() + SecondPoint->Y() - FirstPoint->Y() - FourthPoint->Y())*0.5;
200 //
201 // double lx = std::sqrt(p4_p3_X*p4_p3_X + p4_p3_Y*p4_p3_Y);
202 // double ly = std::sqrt(p2_p1_X*p2_p1_X + p2_p1_Y*p2_p1_Y);
203 // if(lx > ly)
204 // {
205 // this->Points().push_back( typename PointType::Pointer( new PointType( FirstPoint ) ) );
206 // this->Points().push_back( typename PointType::Pointer( new PointType( SecondPoint ) ) );
207 // this->Points().push_back( typename PointType::Pointer( new PointType( ThirdPoint ) ) );
208 // this->Points().push_back( typename PointType::Pointer( new PointType( FourthPoint ) ) );
209 // }
210 // else
211 // {
212 // this->Points().push_back( typename PointType::Pointer( new PointType( FourthPoint ) ) );
213 // this->Points().push_back( typename PointType::Pointer( new PointType( FirstPoint ) ) );
214 // this->Points().push_back( typename PointType::Pointer( new PointType( SecondPoint ) ) );
215 // this->Points().push_back( typename PointType::Pointer( new PointType( ThirdPoint ) ) );
216 // }
217 // }
218 
219  QuadrilateralInterface2D4( typename PointType::Pointer pFirstPoint,
220  typename PointType::Pointer pSecondPoint,
221  typename PointType::Pointer pThirdPoint,
222  typename PointType::Pointer pFourthPoint )
223  : BaseType( PointsArrayType(), &msGeometryData )
224  {
225  double p2_p1_X = (pThirdPoint->X() + pFourthPoint->X() - pFirstPoint->X() - pSecondPoint->X())*0.5;
226  double p2_p1_Y = (pThirdPoint->Y() + pFourthPoint->Y() - pFirstPoint->Y() - pSecondPoint->Y())*0.5;
227  double p4_p3_X = (pThirdPoint->X() + pSecondPoint->X() - pFirstPoint->X() - pFourthPoint->X())*0.5;
228  double p4_p3_Y = (pThirdPoint->Y() + pSecondPoint->Y() - pFirstPoint->Y() - pFourthPoint->Y())*0.5;
229 
230  double lx = std::sqrt(p4_p3_X*p4_p3_X + p4_p3_Y*p4_p3_Y);
231  double ly = std::sqrt(p2_p1_X*p2_p1_X + p2_p1_Y*p2_p1_Y);
232 
233  if(ly < lx)
234  {
235  this->Points().push_back( pFirstPoint );
236  this->Points().push_back( pSecondPoint );
237  this->Points().push_back( pThirdPoint );
238  this->Points().push_back( pFourthPoint );
239  }
240  else
241  {
242  this->Points().push_back( pFourthPoint );
243  this->Points().push_back( pFirstPoint );
244  this->Points().push_back( pSecondPoint );
245  this->Points().push_back( pThirdPoint );
246  }
247  }
248 
250  : BaseType( ThisPoints, &msGeometryData )
251  {
252  if ( this->PointsNumber() != 4 )
253  KRATOS_ERROR << "Invalid points number. Expected 4, given " << this->PointsNumber() << std::endl;
254  }
255 
266  : BaseType( rOther )
267  {
268  }
269 
282  template<class TOtherPointType> QuadrilateralInterface2D4( QuadrilateralInterface2D4<TOtherPointType> const& rOther )
283  : BaseType( rOther )
284  {
285  }
286 
291 
293  {
295  }
296 
298  {
300  }
301 
305 
318  {
319  BaseType::operator=( rOther );
320  return *this;
321  }
322 
334  template<class TOtherPointType>
336  {
337  BaseType::operator=( rOther );
338  return *this;
339  }
340 
344 
345  typename BaseType::Pointer Create( PointsArrayType const& ThisPoints ) const override
346  {
347  return typename BaseType::Pointer( new QuadrilateralInterface2D4( ThisPoints ) );
348  }
349 
355  Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
356  {
357  noalias( rResult ) = ZeroMatrix( 4, 2 );
358  rResult( 0, 0 ) = -1.0;
359  rResult( 0, 1 ) = -1.0;
360  rResult( 1, 0 ) = 1.0;
361  rResult( 1, 1 ) = -1.0;
362  rResult( 2, 0 ) = 1.0;
363  rResult( 2, 1 ) = 1.0;
364  rResult( 3, 0 ) = -1.0;
365  rResult( 3, 1 ) = 1.0;
366  return rResult;
367  }
368 
372 
392  double Length() const override
393  {
396 
397  array_1d<double, 3> v( p1 - p0 );
398  return std::sqrt( v[0]*v[0] + v[1]*v[1] );
399  }
400 
416  double Area() const override
417  {
418  return Length();
419  }
420 
421 
436  double DomainSize() const override
437  {
438  return Length();
439  }
440 
449  bool IsInside(
450  const CoordinatesArrayType& rPoint,
451  CoordinatesArrayType& rResult,
452  const double Tolerance = std::numeric_limits<double>::epsilon()
453  ) const override
454  {
455  this->PointLocalCoordinates( rResult, rPoint );
456 
457  if ( std::abs(rResult[0]) <= (1.0+Tolerance) )
458  {
459  if ( std::abs(rResult[1]) <= (1.0+Tolerance) )
460  {
461  return true;
462  }
463  }
464 
465  return false;
466  }
467 
474  // TODO: Check if correct
476  CoordinatesArrayType& rResult,
477  const CoordinatesArrayType& rPoint
478  ) const override
479  {
480  rResult.clear();
481 
482  array_1d<double, 3> FirstPoint;
483  noalias(FirstPoint) = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
484  array_1d<double, 3> SecondPoint;
485  noalias(SecondPoint) = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 2 ));
486 
487  // Project point
488  double tol = 1e-14; // Tolerance
489 
490  // Normal
492  Normal[0] = SecondPoint[1] - FirstPoint[1];
493  Normal[1] = FirstPoint[0] - SecondPoint[0];
494  double norm = std::sqrt(Normal[0] * Normal[0] + Normal[1] * Normal[1]);
495  Normal /= norm;
496 
497  // Vector point and distance
498  array_1d<double,2> VectorPoint = ZeroVector(2);
499  VectorPoint[0] = rPoint[0] - FirstPoint[0];
500  VectorPoint[1] = rPoint[1] - FirstPoint[1];
501  double dist_proy = VectorPoint[0] * Normal[0] + VectorPoint[1] * Normal[1];
502 
503  if (dist_proy < tol)
504  {
505  double L = Length();
506 
507  double l1 = (rPoint[0] - FirstPoint[0]) * (rPoint[0] - FirstPoint[0])
508  + (rPoint[1] - FirstPoint[1]) * (rPoint[1] - FirstPoint[1]);
509  l1 = std::sqrt(l1);
510 
511  double l2 = (rPoint[0] - SecondPoint[0]) * (rPoint[0] - SecondPoint[0])
512  + (rPoint[1] - SecondPoint[1]) * (rPoint[1] - SecondPoint[1]);
513  l2 = std::sqrt(l2);
514 
515 // std::cout << "L: " << L << " l1: " << l1 << " l2: " << l2 << std::endl;
516 
517  if (l1 <= (L + tol) && l2 <= (L + tol))
518  {
519  rResult[0] = 2.0 * l1/(L + tol) - 1.0;
520  }
521  else
522  {
523  rResult[0] = 2.0; // Out of the line!!! TODO: Check if this value gives problems
524  }
525  }
526  else
527  {
528  rResult[0] = 2.0; // Out of the line!!!
529  }
530 
531  return( rResult );
532  }
533 
555  Matrix& Jacobian( Matrix& rResult,
556  IndexType IntegrationPointIndex,
557  IntegrationMethod ThisMethod ) const override
558  {
561  if(rResult.size1() != 2 || rResult.size2() != 1)
562  rResult.resize(2, 1, false);
563  rResult( 0, 0 ) = ( p1[0] - p0[0] ) * 0.5;
564  rResult( 1, 0 ) = ( p1[1] - p0[1] ) * 0.5;
565  return rResult;
566  }
567 
588  Matrix& Jacobian( Matrix& rResult,
589  IndexType IntegrationPointIndex,
590  IntegrationMethod ThisMethod,
591  const Matrix& DeltaPosition ) const override
592  {
595 
598  dp0[0] = (DeltaPosition(0,0) + DeltaPosition(3,0))*0.5;
599  dp0[1] = (DeltaPosition(0,1) + DeltaPosition(3,1))*0.5;
600  dp1[0] = (DeltaPosition(1,0) + DeltaPosition(2,0))*0.5;
601  dp1[1] = (DeltaPosition(1,1) + DeltaPosition(2,1))*0.5;
602 
603  if(rResult.size1() != 2 || rResult.size2() != 1)
604  rResult.resize(2, 1, false);
605  rResult( 0, 0 ) = ( (p1[0]-dp1[0]) - (p0[0]-dp0[0]) ) * 0.5;
606  rResult( 1, 0 ) = ( (p1[1]-dp1[1]) - (p0[1]-dp0[1]) ) * 0.5;
607 
608  return rResult;
609  }
610 
626  Matrix& Jacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
627  {
630  if(rResult.size1() != 2 || rResult.size2() != 1)
631  rResult.resize(2, 1, false);
632  rResult( 0, 0 ) = ( p1[0] - p0[0] ) * 0.5;
633  rResult( 1, 0 ) = ( p1[1] - p0[1] ) * 0.5;
634  return rResult;
635  }
636 
653  IntegrationMethod ThisMethod ) const override
654  {
655  //workaround by riccardo
656  if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
657  rResult.resize( this->IntegrationPointsNumber( ThisMethod ), false );
658 
659  //for all integration points
660  double detJ = Length()*0.5;
661  for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
662  rResult[pnt] = detJ;
663 
664  return rResult;
665  }
666 
689  double DeterminantOfJacobian( IndexType IntegrationPointIndex,
690  IntegrationMethod ThisMethod ) const override
691  {
692  return Length()*0.5;
693  }
694 
720  double DeterminantOfJacobian( const CoordinatesArrayType& rPoint ) const override
721  {
722  return Length()*0.5;
723  }
724 
747  IntegrationMethod ThisMethod ) const override
748  {
749  KRATOS_ERROR << "Jacobian is not square" << std::endl;
750  return rResult;
751  }
752 
777  IndexType IntegrationPointIndex,
778  IntegrationMethod ThisMethod ) const override
779  {
780  KRATOS_ERROR << "Jacobian is not square" << std::endl;
781  return rResult;
782  }
783 
801  const CoordinatesArrayType& rPoint ) const override
802  {
803  KRATOS_ERROR << "Jacobian is not square" << std::endl;
804  return rResult;
805  }
806 
810 
822  SizeType EdgesNumber() const override
823  {
824  return 4;
825  }
826 
836  {
838  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 0 ), this->pGetPoint( 1 ) ) );
839  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 1 ), this->pGetPoint( 2 ) ) );
840  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 2 ), this->pGetPoint( 3 ) ) );
841  edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 3 ), this->pGetPoint( 0 ) ) );
842  return edges;
843  }
844 
848 
861  double ShapeFunctionValue( IndexType ShapeFunctionIndex,
862  const CoordinatesArrayType& rPoint ) const override
863  {
864  switch ( ShapeFunctionIndex )
865  {
866  case 0:
867  return( 0.25*( 1.0 - rPoint[0] )*( 1.0 - rPoint[1] ) );
868  case 1:
869  return( 0.25*( 1.0 + rPoint[0] )*( 1.0 - rPoint[1] ) );
870  case 2:
871  return( 0.25*( 1.0 + rPoint[0] )*( 1.0 + rPoint[1] ) );
872  case 3:
873  return( 0.25*( 1.0 - rPoint[0] )*( 1.0 + rPoint[1] ) );
874  default:
875  KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
876  }
877 
878  return 0;
879  }
880 
892  Vector& ShapeFunctionsValues (Vector &rResult, const CoordinatesArrayType& rCoordinates) const override
893  {
894  if(rResult.size() != 4) rResult.resize(4,false);
895  rResult[0] = 0.25*( 1.0 - rCoordinates[0] )*( 1.0 - rCoordinates[1] );
896  rResult[1] = 0.25*( 1.0 + rCoordinates[0] )*( 1.0 - rCoordinates[1] );
897  rResult[2] = 0.25*( 1.0 + rCoordinates[0] )*( 1.0 + rCoordinates[1] );
898  rResult[3] = 0.25*( 1.0 - rCoordinates[0] )*( 1.0 + rCoordinates[1] );
899 
900  return rResult;
901  }
902 
921  IntegrationMethod ThisMethod ) const override
922  {
923  const unsigned int integration_points_number =
924  msGeometryData.IntegrationPointsNumber( ThisMethod );
925 
926  if ( integration_points_number == 0 )
927  KRATOS_ERROR << "This integration method is not supported" << *this << std::endl;
928 
929  //workaround by riccardo
930  if ( rResult.size() != integration_points_number )
931  {
932  // KLUDGE: While there is a bug in ublas
933  // vector resize, I have to put this beside resizing!!
934  ShapeFunctionsGradientsType temp( integration_points_number );
935  rResult.swap( temp );
936  }
937 
938  //calculating the local gradients
940  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
941 
942  //getting the inverse jacobian matrices
943  JacobiansType temp( integration_points_number );
944 
945  JacobiansType invJ = InverseOfJacobian( temp, ThisMethod );
946 
947  //loop over all integration points
948  for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
949  {
950  rResult[pnt].resize( 4, 2 ,false);
951 
952  for ( int i = 0; i < 4; i++ )
953  {
954  for ( int j = 0; j < 2; j++ )
955  {
956  rResult[pnt]( i, j ) =
957  ( locG[pnt]( i, 0 ) * invJ[pnt]( j, 0 ) )
958  + ( locG[pnt]( i, 1 ) * invJ[pnt]( j, 1 ) );
959  }
960  }
961  }//end of loop over integration points
962  }
963 
967 
975  std::string Info() const override
976  {
977  return "2 dimensional quadrilateral with four nodes in 2D space";
978  }
979 
986  void PrintInfo( std::ostream& rOStream ) const override
987  {
988  rOStream << "2 dimensional quadrilateral with four nodes in 2D space";
989  }
990 
1005  void PrintData( std::ostream& rOStream ) const override
1006  {
1007  // Base Geometry class PrintData call
1008  BaseType::PrintData( rOStream );
1009  std::cout << std::endl;
1010 
1011  // If the geometry has valid points, calculate and output its data
1012  if (this->AllPointsAreValid()) {
1013  Matrix jacobian;
1014  this->Jacobian( jacobian, PointType() );
1015  rOStream << " Jacobian in the origin\t : " << jacobian;
1016  }
1017  }
1018 
1028  Matrix& ShapeFunctionsLocalGradients( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
1029  {
1030  rResult.resize( 4, 2 ,false);
1031  noalias( rResult ) = ZeroMatrix( 4, 2 );
1032  rResult( 0, 0 ) = -0.25 * ( 1.0 - rPoint[1] );
1033  rResult( 0, 1 ) = -0.25 * ( 1.0 - rPoint[0] );
1034  rResult( 1, 0 ) = 0.25 * ( 1.0 - rPoint[1] );
1035  rResult( 1, 1 ) = -0.25 * ( 1.0 + rPoint[0] );
1036  rResult( 2, 0 ) = 0.25 * ( 1.0 + rPoint[1] );
1037  rResult( 2, 1 ) = 0.25 * ( 1.0 + rPoint[0] );
1038  rResult( 3, 0 ) = -0.25 * ( 1.0 + rPoint[1] );
1039  rResult( 3, 1 ) = 0.25 * ( 1.0 - rPoint[0] );
1040  return rResult;
1041  }
1042 
1050  const CoordinatesArrayType& rPoint ) const override
1051  {
1052  if ( rResult.size() != this->PointsNumber() )
1053  {
1054  // KLUDGE: While there is a bug in
1055  // ublas vector resize, I have to put this beside resizing!!
1057  rResult.swap( temp );
1058  }
1059 
1060  rResult[0].resize( 2, 2 ,false);
1061 
1062  rResult[1].resize( 2, 2 ,false);
1063  rResult[2].resize( 2, 2 ,false);
1064  rResult[3].resize( 2, 2 ,false);
1065  rResult[0]( 0, 0 ) = 0.0;
1066  rResult[0]( 0, 1 ) = 0.25;
1067  rResult[0]( 1, 0 ) = 0.25;
1068  rResult[0]( 1, 1 ) = 0.0;
1069  rResult[1]( 0, 0 ) = 0.0;
1070  rResult[1]( 0, 1 ) = -0.25;
1071  rResult[1]( 1, 0 ) = -0.25;
1072  rResult[1]( 1, 1 ) = 0.0;
1073  rResult[2]( 0, 0 ) = 0.0;
1074  rResult[2]( 0, 1 ) = 0.25;
1075  rResult[2]( 1, 0 ) = 0.25;
1076  rResult[2]( 1, 1 ) = 0.0;
1077  rResult[3]( 0, 0 ) = 0.0;
1078  rResult[3]( 0, 1 ) = -0.25;
1079  rResult[3]( 1, 0 ) = -0.25;
1080  rResult[3]( 1, 1 ) = 0.0;
1081  return rResult;
1082  }
1083 
1091  const CoordinatesArrayType& rPoint ) const override
1092  {
1093  if ( rResult.size() != this->PointsNumber() )
1094  {
1095  // KLUDGE: While there is a bug in
1096  // ublas vector resize, I have to put this beside resizing!!
1097 // ShapeFunctionsGradientsType
1099  rResult.swap( temp );
1100  }
1101 
1102  for ( IndexType i = 0; i < rResult.size(); i++ )
1103  {
1105  rResult[i].swap( temp );
1106  }
1107 
1108  for ( unsigned int i = 0; i < this->PointsNumber(); i++ )
1109  {
1110  for ( unsigned int j = 0; j < 2; j++ )
1111  {
1112  rResult[i][j].resize( 2, 2, false );
1113  noalias( rResult[i][j] ) = ZeroMatrix( 2, 2 );
1114  }
1115  }
1116 
1117 // rResult(0,0).resize( 2, 2);
1118 // rResult(0,1).resize( 2, 2);
1119 // rResult(1,0).resize( 2, 2);
1120 // rResult(1,1).resize( 2, 2);
1121 // rResult(2,0).resize( 2, 2);
1122 // rResult(2,1).resize( 2, 2);
1123 // rResult(3,0).resize( 2, 2);
1124 // rResult(3,0).resize( 2, 2);
1125 
1126  for ( int i = 0; i < 4; i++ )
1127  {
1128  rResult[i][0]( 0, 0 ) = 0.0;
1129  rResult[i][0]( 0, 1 ) = 0.0;
1130  rResult[i][0]( 1, 0 ) = 0.0;
1131  rResult[i][0]( 1, 1 ) = 0.0;
1132  rResult[i][1]( 0, 0 ) = 0.0;
1133  rResult[i][1]( 0, 1 ) = 0.0;
1134  rResult[i][1]( 1, 0 ) = 0.0;
1135  rResult[i][1]( 1, 1 ) = 0.0;
1136  }
1137 
1138  return rResult;
1139  }
1140 
1144 
1146 
1147 protected:
1152 private:
1155 
1156  static const GeometryData msGeometryData;
1157 
1158  static const GeometryDimension msGeometryDimension;
1159 
1163 
1164 
1168 
1169  friend class Serializer;
1170 
1171  void save( Serializer& rSerializer ) const override
1172  {
1174  }
1175 
1176  void load( Serializer& rSerializer ) override
1177  {
1179  }
1180 
1181  QuadrilateralInterface2D4(): BaseType( PointsArrayType(), &msGeometryData ) {}
1182 
1183 
1187 
1188 
1192 
1203  static Matrix CalculateShapeFunctionsIntegrationPointsValues(
1204  typename BaseType::IntegrationMethod ThisMethod )
1205  {
1206  IntegrationPointsContainerType all_integration_points =
1207  AllIntegrationPoints();
1208  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1209  //number of integration points
1210  const int integration_points_number = integration_points.size();
1211  //number of nodes in current geometry
1212  const int points_number = 4;
1213  //setting up return matrix
1214  Matrix shape_function_values( integration_points_number, points_number );
1215  //loop over all integration points
1216 
1217  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1218  {
1219  shape_function_values( pnt, 0 ) =
1220  0.25 * ( 1.0 - integration_points[pnt].X() )
1221  * ( 1.0 - integration_points[pnt].Y() );
1222  shape_function_values( pnt, 1 ) =
1223  0.25 * ( 1.0 + integration_points[pnt].X() )
1224  * ( 1.0 - integration_points[pnt].Y() );
1225  shape_function_values( pnt, 2 ) =
1226  0.25 * ( 1.0 + integration_points[pnt].X() )
1227  * ( 1.0 + integration_points[pnt].Y() );
1228  shape_function_values( pnt, 3 ) =
1229  0.25 * ( 1.0 - integration_points[pnt].X() )
1230  * ( 1.0 + integration_points[pnt].Y() );
1231  }
1232  return shape_function_values;
1233  }
1234 
1247  static ShapeFunctionsGradientsType CalculateShapeFunctionsIntegrationPointsLocalGradients(
1248  typename BaseType::IntegrationMethod ThisMethod )
1249  {
1250  IntegrationPointsContainerType all_integration_points =
1251  AllIntegrationPoints();
1252  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1253  //number of integration points
1254  const int integration_points_number = integration_points.size();
1255  ShapeFunctionsGradientsType d_shape_f_values( integration_points_number );
1256  //initialising container
1257  //std::fill(d_shape_f_values.begin(), d_shape_f_values.end(), Matrix(4,2));
1258  //loop over all integration points
1259 
1260  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1261  {
1262  Matrix result( 4, 2 );
1263  result( 0, 0 ) = -0.25 * ( 1.0 - integration_points[pnt].Y() );
1264  result( 0, 1 ) = -0.25 * ( 1.0 - integration_points[pnt].X() );
1265  result( 1, 0 ) = 0.25 * ( 1.0 - integration_points[pnt].Y() );
1266  result( 1, 1 ) = -0.25 * ( 1.0 + integration_points[pnt].X() );
1267  result( 2, 0 ) = 0.25 * ( 1.0 + integration_points[pnt].Y() );
1268  result( 2, 1 ) = 0.25 * ( 1.0 + integration_points[pnt].X() );
1269  result( 3, 0 ) = -0.25 * ( 1.0 + integration_points[pnt].Y() );
1270  result( 3, 1 ) = 0.25 * ( 1.0 - integration_points[pnt].X() );
1271  d_shape_f_values[pnt] = result;
1272  }
1273  return d_shape_f_values;
1274  }
1275 
1279  static const IntegrationPointsContainerType AllIntegrationPoints()
1280  {
1281  IntegrationPointsContainerType integration_points =
1282  {
1283  {
1284  Quadrature < QuadrilateralGaussLobattoIntegrationPoints1,
1285  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1286  Quadrature < QuadrilateralGaussLobattoIntegrationPoints2,
1287  2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1290  }
1291  };
1292  return integration_points;
1293  }
1294 
1298  static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
1299  {
1300  ShapeFunctionsValuesContainerType shape_functions_values =
1301  {
1302  {
1303  QuadrilateralInterface2D4<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1305  QuadrilateralInterface2D4<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
1307  Matrix(),
1308  Matrix()
1309  }
1310  };
1311  return shape_functions_values;
1312  }
1313 
1318  AllShapeFunctionsLocalGradients()
1319  {
1320  ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
1321  {
1322  {
1323  QuadrilateralInterface2D4<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
1324  QuadrilateralInterface2D4<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
1327  }
1328  };
1329  return shape_functions_local_gradients;
1330  }
1331 
1335 
1336 
1340 
1341 
1345 
1346  template<class TOtherPointType> friend class QuadrilateralInterface2D4;
1347 
1351 
1352 
1354 }; // Class Geometry
1355 
1359 
1360 
1364 
1367 template<class TPointType> inline std::istream& operator >> (
1368  std::istream& rIStream,
1373 template<class TPointType> inline std::ostream& operator << (
1374  std::ostream& rOStream,
1376 {
1377  rThis.PrintInfo( rOStream );
1378  rOStream << std::endl;
1379  rThis.PrintData( rOStream );
1380  return rOStream;
1381 }
1382 
1384 
1385 template<class TPointType> const
1386 GeometryData QuadrilateralInterface2D4<TPointType>::msGeometryData(
1389  QuadrilateralInterface2D4<TPointType>::AllIntegrationPoints(),
1390  QuadrilateralInterface2D4<TPointType>::AllShapeFunctionsValues(),
1391  AllShapeFunctionsLocalGradients()
1392 );
1393 
1394 template<class TPointType>
1396 
1397 }// namespace Kratos.
1398 
1399 #endif // KRATOS_QUADRILATERAL_INTERFACE_2D_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
Geometry & operator=(const Geometry &rOther)
Definition: geometry.h:400
std::vector< IntegrationPointType > IntegrationPointsArrayType
Definition: geometry.h:161
virtual array_1d< double, 3 > Normal(const CoordinatesArrayType &rPointLocalCoordinates) const
It returns a vector that is normal to its corresponding geometry in the given local point.
Definition: geometry.h:1543
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
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
An two node 2D line geometry with linear shape functions.
Definition: line_2d_2.h:65
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_2d_4.h:54
double DeterminantOfJacobian(const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:720
BaseType::IndexType IndexType
Definition: quadrilateral_interface_2d_4.h:97
BaseType::ShapeFunctionsThirdDerivativesType ShapeFunctionsThirdDerivativesType
Definition: quadrilateral_interface_2d_4.h:179
void PrintInfo(std::ostream &rOStream) const override
Definition: quadrilateral_interface_2d_4.h:986
double Area() const override
Definition: quadrilateral_interface_2d_4.h:416
BaseType::GeometriesArrayType GeometriesArrayType
Definition: quadrilateral_interface_2d_4.h:84
Matrix & InverseOfJacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_2d_4.h:776
GeometryData::IntegrationMethod IntegrationMethod
Definition: quadrilateral_interface_2d_4.h:78
CoordinatesArrayType & PointLocalCoordinates(CoordinatesArrayType &rResult, const CoordinatesArrayType &rPoint) const override
Returns the local coordinates of a given arbitrary point.
Definition: quadrilateral_interface_2d_4.h:475
BaseType::IntegrationPointsContainerType IntegrationPointsContainerType
Definition: quadrilateral_interface_2d_4.h:137
QuadrilateralInterface2D4(typename PointType::Pointer pFirstPoint, typename PointType::Pointer pSecondPoint, typename PointType::Pointer pThirdPoint, typename PointType::Pointer pFourthPoint)
Definition: quadrilateral_interface_2d_4.h:219
Vector & DeterminantOfJacobian(Vector &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_2d_4.h:652
BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType
Definition: quadrilateral_interface_2d_4.h:143
ShapeFunctionsSecondDerivativesType & ShapeFunctionsSecondDerivatives(ShapeFunctionsSecondDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:1049
Matrix & PointsLocalCoordinates(Matrix &rResult) const override
Definition: quadrilateral_interface_2d_4.h:355
double DomainSize() const override
Definition: quadrilateral_interface_2d_4.h:436
QuadrilateralInterface2D4 & operator=(const QuadrilateralInterface2D4 &rOther)
Definition: quadrilateral_interface_2d_4.h:317
JacobiansType & InverseOfJacobian(JacobiansType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_2d_4.h:746
BaseType::JacobiansType JacobiansType
Definition: quadrilateral_interface_2d_4.h:156
Matrix & ShapeFunctionsLocalGradients(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:1028
BaseType::PointsArrayType PointsArrayType
Definition: quadrilateral_interface_2d_4.h:110
QuadrilateralInterface2D4(QuadrilateralInterface2D4 const &rOther)
Definition: quadrilateral_interface_2d_4.h:265
void ShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType &rResult, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_2d_4.h:919
GeometriesArrayType GenerateEdges() const override
This method gives you all edges of this geometry.
Definition: quadrilateral_interface_2d_4.h:835
Matrix & InverseOfJacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:800
friend class QuadrilateralInterface2D4
Definition: quadrilateral_interface_2d_4.h:1346
BaseType::IntegrationPointType IntegrationPointType
Definition: quadrilateral_interface_2d_4.h:121
Vector & ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType &rCoordinates) const override
Definition: quadrilateral_interface_2d_4.h:892
Line2D2< TPointType > EdgeType
Definition: quadrilateral_interface_2d_4.h:68
SizeType EdgesNumber() const override
This method gives you number of all edges of this geometry.
Definition: quadrilateral_interface_2d_4.h:822
double ShapeFunctionValue(IndexType ShapeFunctionIndex, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:861
TPointType PointType
Definition: quadrilateral_interface_2d_4.h:89
std::string Info() const override
Definition: quadrilateral_interface_2d_4.h:975
virtual ~QuadrilateralInterface2D4()
Definition: quadrilateral_interface_2d_4.h:290
void PrintData(std::ostream &rOStream) const override
Definition: quadrilateral_interface_2d_4.h:1005
BaseType::CoordinatesArrayType CoordinatesArrayType
Definition: quadrilateral_interface_2d_4.h:115
BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: quadrilateral_interface_2d_4.h:163
Geometry< TPointType > BaseType
Definition: quadrilateral_interface_2d_4.h:63
double Length() const override
Definition: quadrilateral_interface_2d_4.h:392
GeometryData::KratosGeometryType GetGeometryType() const override
Definition: quadrilateral_interface_2d_4.h:297
BaseType::NormalType NormalType
Definition: quadrilateral_interface_2d_4.h:184
Matrix & Jacobian(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:626
BaseType::ShapeFunctionsSecondDerivativesType ShapeFunctionsSecondDerivativesType
Definition: quadrilateral_interface_2d_4.h:171
ShapeFunctionsThirdDerivativesType & ShapeFunctionsThirdDerivatives(ShapeFunctionsThirdDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: quadrilateral_interface_2d_4.h:1090
Matrix & Jacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_2d_4.h:555
GeometryData::KratosGeometryFamily GetGeometryFamily() const override
Definition: quadrilateral_interface_2d_4.h:292
BaseType::Pointer Create(PointsArrayType const &ThisPoints) const override
Creates a new geometry pointer.
Definition: quadrilateral_interface_2d_4.h:345
BaseType::SizeType SizeType
Definition: quadrilateral_interface_2d_4.h:104
BaseType::IntegrationPointsArrayType IntegrationPointsArrayType
Definition: quadrilateral_interface_2d_4.h:130
Matrix & Jacobian(Matrix &rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod, const Matrix &DeltaPosition) const override
Definition: quadrilateral_interface_2d_4.h:588
QuadrilateralInterface2D4(const PointsArrayType &ThisPoints)
Definition: quadrilateral_interface_2d_4.h:249
KRATOS_CLASS_POINTER_DEFINITION(QuadrilateralInterface2D4)
BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: quadrilateral_interface_2d_4.h:149
QuadrilateralInterface2D4 & operator=(QuadrilateralInterface2D4< TOtherPointType > const &rOther)
Definition: quadrilateral_interface_2d_4.h:335
QuadrilateralInterface2D4(QuadrilateralInterface2D4< TOtherPointType > const &rOther)
Definition: quadrilateral_interface_2d_4.h:282
double DeterminantOfJacobian(IndexType IntegrationPointIndex, IntegrationMethod ThisMethod) const override
Definition: quadrilateral_interface_2d_4.h:689
bool IsInside(const CoordinatesArrayType &rPoint, CoordinatesArrayType &rResult, const double Tolerance=std::numeric_limits< double >::epsilon()) const override
Definition: quadrilateral_interface_2d_4.h:449
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 clear()
Definition: array_1d.h:325
#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
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
const GeometryData QuadrilateralInterface2D4< TPointType >::msGeometryData & msGeometryDimension(), QuadrilateralInterface2D4< TPointType >::AllShapeFunctionsValues(), AllShapeFunctionsLocalGradients()
Definition: brep_curve.h:483
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
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
v
Definition: generate_convection_diffusion_explicit_element.py:114
int tol
Definition: hinsberg_optimization.py:138
def load(f)
Definition: ode_solve.py:307
int L
Definition: ode_solve.py:390
int j
Definition: quadrature.py:648
float temp
Definition: rotating_cone.py:85
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31