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.
hexahedra_3d_8.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Riccardo Rossi
11 // Janosch Stascheit
12 // Felix Nagel
13 // contributors: Hoang Giang Bui
14 // Josep Maria Carbonell
15 // Bodhinanda Chandra
16 //
17 
18 #pragma once
19 
20 // System includes
21 
22 // External includes
23 
24 // Project includes
27 #include "utilities/geometry_utilities.h"
30 
31 namespace Kratos
32 {
54 template<class TPointType> class Hexahedra3D8 : public Geometry<TPointType>
55 {
56 public:
65 
71 
76 
81 
87 
91  typedef TPointType PointType;
92 
98  typedef typename BaseType::IndexType IndexType;
99 
100 
106  typedef typename BaseType::SizeType SizeType;
107 
113 
120 
128 
135 
142 
149 
156 
163 
170 
175 
180 
185 
186 
190  Hexahedra3D8( const PointType& Point1, const PointType& Point2,
191  const PointType& Point3, const PointType& Point4,
192  const PointType& Point5, const PointType& Point6,
193  const PointType& Point7, const PointType& Point8 )
194  : BaseType( PointsArrayType(), &msGeometryData )
195  {
196  this->Points().push_back( typename PointType::Pointer( new PointType( Point1 ) ) );
197  this->Points().push_back( typename PointType::Pointer( new PointType( Point2 ) ) );
198  this->Points().push_back( typename PointType::Pointer( new PointType( Point3 ) ) );
199  this->Points().push_back( typename PointType::Pointer( new PointType( Point4 ) ) );
200  this->Points().push_back( typename PointType::Pointer( new PointType( Point5 ) ) );
201  this->Points().push_back( typename PointType::Pointer( new PointType( Point6 ) ) );
202  this->Points().push_back( typename PointType::Pointer( new PointType( Point7 ) ) );
203  this->Points().push_back( typename PointType::Pointer( new PointType( Point8 ) ) );
204  }
205 
206  Hexahedra3D8( 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  : BaseType( PointsArrayType(), &msGeometryData )
215  {
216  this->Points().push_back( pPoint1 );
217  this->Points().push_back( pPoint2 );
218  this->Points().push_back( pPoint3 );
219  this->Points().push_back( pPoint4 );
220  this->Points().push_back( pPoint5 );
221  this->Points().push_back( pPoint6 );
222  this->Points().push_back( pPoint7 );
223  this->Points().push_back( pPoint8 );
224  }
225 
226  explicit Hexahedra3D8( const PointsArrayType& ThisPoints )
227  : BaseType( ThisPoints, &msGeometryData )
228  {
229  if ( this->PointsNumber() != 8 )
230  KRATOS_ERROR << "Invalid points number. Expected 8, given " << this->PointsNumber() << std::endl;
231  }
232 
234  explicit Hexahedra3D8(
235  const IndexType GeometryId,
236  const PointsArrayType& rThisPoints
237  ) : BaseType( GeometryId, rThisPoints, &msGeometryData )
238  {
239  KRATOS_ERROR_IF( this->PointsNumber() != 8 ) << "Invalid points number. Expected 8, given " << this->PointsNumber() << std::endl;
240  }
241 
243  explicit Hexahedra3D8(
244  const std::string& rGeometryName,
245  const PointsArrayType& rThisPoints
246  ) : BaseType( rGeometryName, rThisPoints, &msGeometryData)
247  {
248  KRATOS_ERROR_IF(this->PointsNumber() != 8) << "Invalid points number. Expected 8, given " << this->PointsNumber() << std::endl;
249  }
250 
260  Hexahedra3D8( Hexahedra3D8 const& rOther )
261  : BaseType( rOther )
262  {
263  }
264 
277  template<class TOtherPointType> explicit Hexahedra3D8( Hexahedra3D8<TOtherPointType> const& rOther )
278  : BaseType( rOther )
279  {
280  }
281 
283  ~Hexahedra3D8() override {}
284 
286  {
288  }
289 
291  {
293  }
294 
311  {
312  BaseType::operator=( rOther );
313  return *this;
314  }
315 
327  template<class TOtherPointType>
329  {
330  BaseType::operator=( rOther );
331 
332  return *this;
333  }
334 
335 
346  typename BaseType::Pointer Create(
347  const IndexType NewGeometryId,
348  PointsArrayType const& rThisPoints
349  ) const override
350  {
351  return typename BaseType::Pointer( new Hexahedra3D8( NewGeometryId, rThisPoints ) );
352  }
353 
360  typename BaseType::Pointer Create(
361  const IndexType NewGeometryId,
362  const BaseType& rGeometry
363  ) const override
364  {
365  auto p_geometry = typename BaseType::Pointer( new Hexahedra3D8( NewGeometryId, rGeometry.Points() ) );
366  p_geometry->SetData(rGeometry.GetData());
367  return p_geometry;
368  }
369 
389  double Length() const override
390  {
391  return std::sqrt( std::abs( this->DeterminantOfJacobian( PointType() ) ) );
392  }
393 
408  double Area() const override
409  {
410  return Volume();
411  }
412 
421  double Volume() const override
422  {
424  }
425 
439  double DomainSize() const override
440  {
441  return Volume();
442  }
443 
449  Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
450  {
451  if ( rResult.size1() != 8 || rResult.size2() != 3 )
452  rResult.resize( 8, 3, false );
453 
454  rResult( 0, 0 ) = -1.0;
455  rResult( 0, 1 ) = -1.0;
456  rResult( 0, 2 ) = -1.0;
457 
458  rResult( 1, 0 ) = 1.0;
459  rResult( 1, 1 ) = -1.0;
460  rResult( 1, 2 ) = -1.0;
461 
462  rResult( 2, 0 ) = 1.0;
463  rResult( 2, 1 ) = 1.0;
464  rResult( 2, 2 ) = -1.0;
465 
466  rResult( 3, 0 ) = -1.0;
467  rResult( 3, 1 ) = 1.0;
468  rResult( 3, 2 ) = -1.0;
469 
470  rResult( 4, 0 ) = -1.0;
471  rResult( 4, 1 ) = -1.0;
472  rResult( 4, 2 ) = 1.0;
473 
474  rResult( 5, 0 ) = 1.0;
475  rResult( 5, 1 ) = -1.0;
476  rResult( 5, 2 ) = 1.0;
477 
478  rResult( 6, 0 ) = 1.0;
479  rResult( 6, 1 ) = 1.0;
480  rResult( 6, 2 ) = 1.0;
481 
482  rResult( 7, 0 ) = -1.0;
483  rResult( 7, 1 ) = 1.0;
484  rResult( 7, 2 ) = 1.0;
485 
486  return rResult;
487  }
488 
497  bool IsInside(
498  const CoordinatesArrayType& rPoint,
499  CoordinatesArrayType& rResult,
500  const double Tolerance = std::numeric_limits<double>::epsilon()
501  ) const override
502  {
503  this->PointLocalCoordinates( rResult, rPoint );
504 
505  if ( std::abs( rResult[0] ) <= (1.0 + Tolerance) ) {
506  if ( std::abs( rResult[1] ) <= (1.0 + Tolerance) ) {
507  if ( std::abs( rResult[2] ) <= (1.0 + Tolerance) ) {
508  return true;
509  }
510  }
511  }
512 
513  return false;
514  }
515 
516 
523  // will be used by refinement algorithm, thus uncommented. janosch.
524  SizeType EdgesNumber() const override
525  {
526  return 12;
527  }
528 
536  {
538  typedef typename Geometry<TPointType>::Pointer EdgePointerType;
539  edges.push_back( EdgePointerType( new EdgeType(
540  this->pGetPoint( 0 ),
541  this->pGetPoint( 1 ) ) ) );
542  edges.push_back( EdgePointerType( new EdgeType(
543  this->pGetPoint( 1 ),
544  this->pGetPoint( 2 ) ) ) );
545  edges.push_back( EdgePointerType( new EdgeType(
546  this->pGetPoint( 2 ),
547  this->pGetPoint( 3 ) ) ) );
548  edges.push_back( EdgePointerType( new EdgeType(
549  this->pGetPoint( 3 ),
550  this->pGetPoint( 0 ) ) ) );
551  edges.push_back( EdgePointerType( new EdgeType(
552  this->pGetPoint( 4 ),
553  this->pGetPoint( 5 ) ) ) );
554  edges.push_back( EdgePointerType( new EdgeType(
555  this->pGetPoint( 5 ),
556  this->pGetPoint( 6 ) ) ) );
557  edges.push_back( EdgePointerType( new EdgeType(
558  this->pGetPoint( 6 ),
559  this->pGetPoint( 7 ) ) ) );
560  edges.push_back( EdgePointerType( new EdgeType(
561  this->pGetPoint( 7 ),
562  this->pGetPoint( 4 ) ) ) );
563  edges.push_back( EdgePointerType( new EdgeType(
564  this->pGetPoint( 0 ),
565  this->pGetPoint( 4 ) ) ) );
566  edges.push_back( EdgePointerType( new EdgeType(
567  this->pGetPoint( 1 ),
568  this->pGetPoint( 5 ) ) ) );
569  edges.push_back( EdgePointerType( new EdgeType(
570  this->pGetPoint( 2 ),
571  this->pGetPoint( 6 ) ) ) );
572  edges.push_back( EdgePointerType( new EdgeType(
573  this->pGetPoint( 3 ),
574  this->pGetPoint( 7 ) ) ) );
575  return edges;
576  }
577 
583  double AverageEdgeLength() const override
584  {
585  const TPointType& p0 = this->GetPoint(0);
586  const TPointType& p1 = this->GetPoint(1);
587  const TPointType& p2 = this->GetPoint(2);
588  const TPointType& p3 = this->GetPoint(3);
589  const TPointType& p4 = this->GetPoint(4);
590  const TPointType& p5 = this->GetPoint(5);
591  const TPointType& p6 = this->GetPoint(6);
592  const TPointType& p7 = this->GetPoint(7);
593  return (MathUtils<double>::Norm3(p0-p1) +
594  MathUtils<double>::Norm3(p1-p2) +
595  MathUtils<double>::Norm3(p2-p3) +
596  MathUtils<double>::Norm3(p3-p0) +
597  MathUtils<double>::Norm3(p4-p5) +
598  MathUtils<double>::Norm3(p5-p6) +
599  MathUtils<double>::Norm3(p6-p7) +
600  MathUtils<double>::Norm3(p7-p4) +
601  MathUtils<double>::Norm3(p0-p4) +
602  MathUtils<double>::Norm3(p1-p5) +
603  MathUtils<double>::Norm3(p2-p6) +
604  MathUtils<double>::Norm3(p3-p7)) /12.0;
605  }
606 
615  double MaxEdgeLength() const override {
616  const auto edges = GenerateEdges();
617  double max_edge_length = 0.0;
618  for (const auto& r_edge: edges) {
619  max_edge_length = std::max(max_edge_length, r_edge.Length());
620  }
621  return max_edge_length;
622  }
623 
632  double MinEdgeLength() const override {
633  const auto edges = GenerateEdges();
634  double min_edge_length = std::numeric_limits<double>::max();
635  for (const auto& r_edge: edges) {
636  min_edge_length = std::min(min_edge_length, r_edge.Length());
637  }
638  return min_edge_length;
639  }
640 
644 
652  SizeType FacesNumber() const override
653  {
654  return 6;
655  }
656 
666  {
668  typedef typename Geometry<TPointType>::Pointer FacePointerType;
669  faces.push_back( FacePointerType( new FaceType(
670  this->pGetPoint( 3 ),
671  this->pGetPoint( 2 ),
672  this->pGetPoint( 1 ),
673  this->pGetPoint( 0 ) ) ) );
674  faces.push_back( FacePointerType( new FaceType(
675  this->pGetPoint( 0 ),
676  this->pGetPoint( 1 ),
677  this->pGetPoint( 5 ),
678  this->pGetPoint( 4 ) ) ) );
679  faces.push_back( FacePointerType( new FaceType(
680  this->pGetPoint( 2 ),
681  this->pGetPoint( 6 ),
682  this->pGetPoint( 5 ),
683  this->pGetPoint( 1 ) ) ) );
684  faces.push_back( FacePointerType( new FaceType(
685  this->pGetPoint( 7 ),
686  this->pGetPoint( 6 ),
687  this->pGetPoint( 2 ),
688  this->pGetPoint( 3 ) ) ) );
689  faces.push_back( FacePointerType( new FaceType(
690  this->pGetPoint( 7 ),
691  this->pGetPoint( 3 ),
692  this->pGetPoint( 0 ),
693  this->pGetPoint( 4 ) ) ) );
694  faces.push_back( FacePointerType( new FaceType(
695  this->pGetPoint( 4 ),
696  this->pGetPoint( 5 ),
697  this->pGetPoint( 6 ),
698  this->pGetPoint( 7 ) ) ) );
699  return faces;
700  }
701 
702 
703  bool HasIntersection( const Point& rLowPoint, const Point& rHighPoint ) const override
704  {
705  using Quadrilateral3D4Type = Quadrilateral3D4<TPointType>;
706  // Check if faces have intersection
707  if(Quadrilateral3D4Type(this->pGetPoint(3),this->pGetPoint(2), this->pGetPoint(1), this->pGetPoint(0)).HasIntersection(rLowPoint, rHighPoint))
708  return true;
709  if(Quadrilateral3D4Type(this->pGetPoint(0),this->pGetPoint(1), this->pGetPoint(5), this->pGetPoint(4)).HasIntersection(rLowPoint, rHighPoint))
710  return true;
711  if(Quadrilateral3D4Type(this->pGetPoint(2),this->pGetPoint(6), this->pGetPoint(5), this->pGetPoint(1)).HasIntersection(rLowPoint, rHighPoint))
712  return true;
713  if(Quadrilateral3D4Type(this->pGetPoint(7),this->pGetPoint(6), this->pGetPoint(2), this->pGetPoint(3)).HasIntersection(rLowPoint, rHighPoint))
714  return true;
715  if(Quadrilateral3D4Type(this->pGetPoint(7),this->pGetPoint(3), this->pGetPoint(0), this->pGetPoint(4)).HasIntersection(rLowPoint, rHighPoint))
716  return true;
717  if(Quadrilateral3D4Type(this->pGetPoint(4),this->pGetPoint(5), this->pGetPoint(6), this->pGetPoint(7)).HasIntersection(rLowPoint, rHighPoint))
718  return true;
719 
720  CoordinatesArrayType local_coordinates;
721  // if there are no faces intersecting the box then or the box is inside the hexahedron or it does not have intersection
722  if(IsInside(rLowPoint,local_coordinates))
723  return true;
724 
725  return false;
726  }
727 
733  void ComputeSolidAngles(Vector& rSolidAngles) const override
734  {
735  if(rSolidAngles.size() != 8) {
736  rSolidAngles.resize(8, false);
737  }
738 
739  Vector dihedral_angles(24);
740  ComputeDihedralAngles(dihedral_angles);
741 
742  for (unsigned int i = 0; i < 8; ++i) {
743  rSolidAngles[i] = dihedral_angles[3*i]
744  + dihedral_angles[3*i + 1]
745  + dihedral_angles[3*i + 2]
746  - Globals::Pi;
747  }
748  }
749 
756  void ComputeDihedralAngles(Vector& rDihedralAngles) const override
757  {
758  if(rDihedralAngles.size() != 24) {
759  rDihedralAngles.resize(24, false);
760  }
761  const auto faces = this->GenerateFaces();
762  // The three faces that contain the node i can be obtained by doing:
763  // faces[faces_0[i]], faces[faces_1[i]], faces[faces_2[i]]
764  const std::array<unsigned int, 8> faces_0 = {0,0,0,0,5,5,5,5};
765  const std::array<unsigned int, 8> faces_1 = {1,1,3,3,1,1,3,3};
766  const std::array<unsigned int, 8> faces_2 = {4,2,2,4,4,2,2,4};
767 
768  array_1d<double, 3> normal_0, normal_1, normal_2;
769  double dihedral_angle_0, dihedral_angle_1, dihedral_angle_2;
770  for (unsigned int i = 0; i < 8; ++i) {
771  const TPointType& r_point_i = this->GetPoint(i);
772  noalias(normal_0) = faces[faces_0[i]].UnitNormal(r_point_i);
773  noalias(normal_1) = faces[faces_1[i]].UnitNormal(r_point_i);
774  noalias(normal_2) = faces[faces_2[i]].UnitNormal(r_point_i);
775  dihedral_angle_0 = std::acos(inner_prod(normal_0, -normal_1));
776  dihedral_angle_1 = std::acos(inner_prod(normal_0, -normal_2));
777  dihedral_angle_2 = std::acos(inner_prod(normal_2, -normal_1));
778  rDihedralAngles[i*3] = dihedral_angle_0;
779  rDihedralAngles[i*3 + 1] = dihedral_angle_1;
780  rDihedralAngles[i*3 + 2] = dihedral_angle_2;
781  }
782  }
783 
791  double MinDihedralAngle() const override {
792  Vector dihedral_angles(24);
793  ComputeDihedralAngles(dihedral_angles);
794  double min_dihedral_angle = dihedral_angles[0];
795  for (unsigned int i = 1; i < 24; i++) {
796  min_dihedral_angle = std::min(dihedral_angles[i], min_dihedral_angle);
797  }
798  return min_dihedral_angle;
799  }
800 
808  double MaxDihedralAngle() const override {
809  Vector dihedral_angles(24);
810  ComputeDihedralAngles(dihedral_angles);
811  double max_dihedral_angle = dihedral_angles[0];
812  for (unsigned int i = 1; i < 24; i++) {
813  max_dihedral_angle = std::max(dihedral_angles[i], max_dihedral_angle);
814  }
815  return max_dihedral_angle;
816  }
817 
829  double VolumeToRMSEdgeLength() const override {
830  const auto edges = GenerateEdges();
831  double sum_squared_lengths = 0.0;
832  for (const auto& r_edge : edges) {
833  const double length = r_edge.Length();
834  sum_squared_lengths += length*length;
835  }
836 
837  const double rms_edge = std::sqrt(1.0/12.0 * sum_squared_lengths);
838 
839  return Volume() / std::pow(rms_edge, 3.0);
840  }
841 
853  double ShortestToLongestEdgeQuality() const override {
854  const auto edges = GenerateEdges();
855  double min_edge_length = std::numeric_limits<double>::max();
856  double max_edge_length = -std::numeric_limits<double>::max();
857  for (const auto& r_edge: edges) {
858  min_edge_length = std::min(min_edge_length, r_edge.Length());
859  max_edge_length = std::max(max_edge_length, r_edge.Length());
860  }
861  return min_edge_length / max_edge_length;
862  }
863 
878  double ShapeFunctionValue( IndexType ShapeFunctionIndex,
879  const CoordinatesArrayType& rPoint ) const override
880  {
881  switch ( ShapeFunctionIndex )
882  {
883  case 0:
884  return( 0.125*( 1.0 - rPoint[0] )*( 1.0 - rPoint[1] )*( 1.0 - rPoint[2] ) );
885  case 1:
886  return( 0.125*( 1.0 + rPoint[0] )*( 1.0 - rPoint[1] )*( 1.0 - rPoint[2] ) );
887  case 2:
888  return( 0.125*( 1.0 + rPoint[0] )*( 1.0 + rPoint[1] )*( 1.0 - rPoint[2] ) );
889  case 3:
890  return( 0.125*( 1.0 - rPoint[0] )*( 1.0 + rPoint[1] )*( 1.0 - rPoint[2] ) );
891  case 4:
892  return( 0.125*( 1.0 - rPoint[0] )*( 1.0 - rPoint[1] )*( 1.0 + rPoint[2] ) );
893  case 5:
894  return( 0.125*( 1.0 + rPoint[0] )*( 1.0 - rPoint[1] )*( 1.0 + rPoint[2] ) );
895  case 6:
896  return( 0.125*( 1.0 + rPoint[0] )*( 1.0 + rPoint[1] )*( 1.0 + rPoint[2] ) );
897  case 7:
898  return( 0.125*( 1.0 - rPoint[0] )*( 1.0 + rPoint[1] )*( 1.0 + rPoint[2] ) );
899  default:
900  KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
901  }
902 
903  return 0;
904  }
905 
917  Vector& ShapeFunctionsValues (Vector &rResult, const CoordinatesArrayType& rCoordinates) const override
918  {
919  if(rResult.size() != 8) rResult.resize(8,false);
920  rResult[0] = 0.125*( 1.0 - rCoordinates[0] )*( 1.0 - rCoordinates[1] )*( 1.0 - rCoordinates[2] ) ;
921  rResult[1] = 0.125*( 1.0 + rCoordinates[0] )*( 1.0 - rCoordinates[1] )*( 1.0 - rCoordinates[2] ) ;
922  rResult[2] = 0.125*( 1.0 + rCoordinates[0] )*( 1.0 + rCoordinates[1] )*( 1.0 - rCoordinates[2] ) ;
923  rResult[3] = 0.125*( 1.0 - rCoordinates[0] )*( 1.0 + rCoordinates[1] )*( 1.0 - rCoordinates[2] ) ;
924  rResult[4] = 0.125*( 1.0 - rCoordinates[0] )*( 1.0 - rCoordinates[1] )*( 1.0 + rCoordinates[2] ) ;
925  rResult[5] = 0.125*( 1.0 + rCoordinates[0] )*( 1.0 - rCoordinates[1] )*( 1.0 + rCoordinates[2] ) ;
926  rResult[6] = 0.125*( 1.0 + rCoordinates[0] )*( 1.0 + rCoordinates[1] )*( 1.0 + rCoordinates[2] ) ;
927  rResult[7] = 0.125*( 1.0 - rCoordinates[0] )*( 1.0 + rCoordinates[1] )*( 1.0 + rCoordinates[2] ) ;
928  return rResult;
929  }
930 
939  Matrix& ShapeFunctionsLocalGradients( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
940  {
941  if ( rResult.size1() != 8 || rResult.size2() != 3 )
942  rResult.resize( 8, 3, false );
943 
944  rResult( 0, 0 ) = -0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 - rPoint[2] );
945  rResult( 0, 1 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[2] );
946  rResult( 0, 2 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[1] );
947  rResult( 1, 0 ) = 0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 - rPoint[2] );
948  rResult( 1, 1 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[2] );
949  rResult( 1, 2 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[1] );
950  rResult( 2, 0 ) = 0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 - rPoint[2] );
951  rResult( 2, 1 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[2] );
952  rResult( 2, 2 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[1] );
953  rResult( 3, 0 ) = -0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 - rPoint[2] );
954  rResult( 3, 1 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[2] );
955  rResult( 3, 2 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[1] );
956  rResult( 4, 0 ) = -0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 + rPoint[2] );
957  rResult( 4, 1 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[2] );
958  rResult( 4, 2 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[1] );
959  rResult( 5, 0 ) = 0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 + rPoint[2] );
960  rResult( 5, 1 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[2] );
961  rResult( 5, 2 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[1] );
962  rResult( 6, 0 ) = 0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 + rPoint[2] );
963  rResult( 6, 1 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[2] );
964  rResult( 6, 2 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[1] );
965  rResult( 7, 0 ) = -0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 + rPoint[2] );
966  rResult( 7, 1 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[2] );
967  rResult( 7, 2 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[1] );
968 
969  return rResult;
970  }
971 
972 
982  {
983  if ( rResult.size() != this->PointsNumber() )
984  {
985  // KLUDGE: While there is a bug in
986  // ublas vector resize, I have to put this beside resizing!!
988  rResult.swap( temp );
989  }
990 
991  for ( unsigned int i = 0; i < this->PointsNumber(); ++i )
992  {
993  rResult[i].resize(3, 3, false);
994  }
995 
996  rResult[0]( 0, 0 ) = 0.0;
997  rResult[0]( 0, 1 ) = 0.125 * ( 1.0 - rPoint[2] );
998  rResult[0]( 0, 2 ) = 0.125 * ( 1.0 - rPoint[1] );
999  rResult[0]( 1, 0 ) = 0.125 * ( 1.0 - rPoint[2] );
1000  rResult[0]( 1, 1 ) = 0.0;
1001  rResult[0]( 1, 2 ) = 0.125 * ( 1.0 - rPoint[0] );
1002  rResult[0]( 2, 0 ) = 0.125 * ( 1.0 - rPoint[1] );
1003  rResult[0]( 2, 1 ) = 0.125 * ( 1.0 - rPoint[0] );
1004  rResult[0]( 2, 2 ) = 0.0;
1005 
1006  rResult[1]( 0, 0 ) = 0.0;
1007  rResult[1]( 0, 1 ) = -0.125 * ( 1.0 - rPoint[2] );
1008  rResult[1]( 0, 2 ) = -0.125 * ( 1.0 - rPoint[1] );
1009  rResult[1]( 1, 0 ) = -0.125 * ( 1.0 - rPoint[2] );
1010  rResult[1]( 1, 1 ) = 0.0;
1011  rResult[1]( 1, 2 ) = 0.125 * ( 1.0 + rPoint[0] );
1012  rResult[1]( 2, 0 ) = -0.125 * ( 1.0 - rPoint[1] );
1013  rResult[1]( 2, 1 ) = 0.125 * ( 1.0 + rPoint[0] );
1014  rResult[1]( 2, 2 ) = 0.0;
1015 
1016  rResult[2]( 0, 0 ) = 0.0;
1017  rResult[2]( 0, 1 ) = 0.125 * ( 1.0 - rPoint[2] );
1018  rResult[2]( 0, 2 ) = -0.125 * ( 1.0 + rPoint[1] );
1019  rResult[2]( 1, 0 ) = 0.125 * ( 1.0 - rPoint[2] );
1020  rResult[2]( 1, 1 ) = 0.0;
1021  rResult[2]( 1, 2 ) = -0.125 * ( 1.0 + rPoint[0] );
1022  rResult[2]( 2, 0 ) = -0.125 * ( 1.0 + rPoint[1] );
1023  rResult[2]( 2, 1 ) = -0.125 * ( 1.0 + rPoint[0] );
1024  rResult[2]( 2, 2 ) = 0.0;
1025 
1026  rResult[3]( 0, 0 ) = 0.0;
1027  rResult[3]( 0, 1 ) = -0.125 * ( 1.0 - rPoint[2] );
1028  rResult[3]( 0, 2 ) = 0.125 * ( 1.0 + rPoint[1] );
1029  rResult[3]( 1, 0 ) = -0.125 * ( 1.0 - rPoint[2] );
1030  rResult[3]( 1, 1 ) = 0.0;
1031  rResult[3]( 1, 2 ) = -0.125 * ( 1.0 - rPoint[0] );
1032  rResult[3]( 2, 0 ) = 0.125 * ( 1.0 + rPoint[1] );
1033  rResult[3]( 2, 1 ) = -0.125 * ( 1.0 - rPoint[0] );
1034  rResult[3]( 2, 2 ) = 0.0;
1035 
1036  rResult[4]( 0, 0 ) = 0.0;
1037  rResult[4]( 0, 1 ) = 0.125 * ( 1.0 + rPoint[2] );
1038  rResult[4]( 0, 2 ) = -0.125 * ( 1.0 - rPoint[1] );
1039  rResult[4]( 1, 0 ) = 0.125 * ( 1.0 + rPoint[2] );
1040  rResult[4]( 1, 1 ) = 0.0;
1041  rResult[4]( 1, 2 ) = -0.125 * ( 1.0 - rPoint[0] );
1042  rResult[4]( 2, 0 ) = -0.125 * ( 1.0 - rPoint[1] );
1043  rResult[4]( 2, 1 ) = -0.125 * ( 1.0 - rPoint[0] );
1044  rResult[4]( 2, 2 ) = 0.0;
1045 
1046  rResult[5]( 0, 0 ) = 0.0;
1047  rResult[5]( 0, 1 ) = -0.125 * ( 1.0 + rPoint[2] );
1048  rResult[5]( 0, 2 ) = 0.125 * ( 1.0 - rPoint[1] );
1049  rResult[5]( 1, 0 ) = -0.125 * ( 1.0 + rPoint[2] );
1050  rResult[5]( 1, 1 ) = 0.0;
1051  rResult[5]( 1, 2 ) = -0.125 * ( 1.0 + rPoint[0] );
1052  rResult[5]( 2, 0 ) = 0.125 * ( 1.0 - rPoint[1] );
1053  rResult[5]( 2, 1 ) = -0.125 * ( 1.0 + rPoint[0] );
1054  rResult[5]( 2, 2 ) = 0.0;
1055 
1056  rResult[6]( 0, 0 ) = 0.0;
1057  rResult[6]( 0, 1 ) = 0.125 * ( 1.0 + rPoint[2] );
1058  rResult[6]( 0, 2 ) = 0.125 * ( 1.0 + rPoint[1] );
1059  rResult[6]( 1, 0 ) = 0.125 * ( 1.0 + rPoint[2] );
1060  rResult[6]( 1, 1 ) = 0.0;
1061  rResult[6]( 1, 2 ) = 0.125 * ( 1.0 + rPoint[0] );
1062  rResult[6]( 2, 0 ) = 0.125 * ( 1.0 + rPoint[1] );
1063  rResult[6]( 2, 1 ) = 0.125 * ( 1.0 + rPoint[0] );
1064  rResult[6]( 2, 2 ) = 0.0;
1065 
1066  rResult[7]( 0, 0 ) = 0.0;
1067  rResult[7]( 0, 1 ) = -0.125 * ( 1.0 + rPoint[2] );
1068  rResult[7]( 0, 2 ) = -0.125 * ( 1.0 + rPoint[1] );
1069  rResult[7]( 1, 0 ) = -0.125 * ( 1.0 + rPoint[2] );
1070  rResult[7]( 1, 1 ) = 0.0;
1071  rResult[7]( 1, 2 ) = 0.125 * ( 1.0 - rPoint[0] );
1072  rResult[7]( 2, 0 ) = -0.125 * ( 1.0 + rPoint[1] );
1073  rResult[7]( 2, 1 ) = 0.125 * ( 1.0 - rPoint[0] );
1074  rResult[7]( 2, 2 ) = 0.0;
1075 
1076  return rResult;
1077  }
1078 
1082 
1096  const CoordinatesArrayType& rPointGlobalCoordinates,
1097  const double Tolerance = std::numeric_limits<double>::epsilon()
1098  ) const override
1099  {
1100  // Point to compute distance to
1101  const Point point(rPointGlobalCoordinates);
1102 
1103  // Check if point is inside
1104  CoordinatesArrayType aux_coordinates;
1105  if (this->IsInside(rPointGlobalCoordinates, aux_coordinates, Tolerance)) {
1106  return 0.0;
1107  }
1108 
1109  // Compute distance to faces
1110  std::array<double, 6> distances;
1111  distances[0] = GeometryUtils::PointDistanceToQuadrilateral3D(this->GetPoint(3), this->GetPoint(2), this->GetPoint(1), this->GetPoint(0), point);
1112  distances[1] = GeometryUtils::PointDistanceToQuadrilateral3D(this->GetPoint(0), this->GetPoint(1), this->GetPoint(5), this->GetPoint(4), point);
1113  distances[2] = GeometryUtils::PointDistanceToQuadrilateral3D(this->GetPoint(2), this->GetPoint(6), this->GetPoint(5), this->GetPoint(1), point);
1114  distances[3] = GeometryUtils::PointDistanceToQuadrilateral3D(this->GetPoint(7), this->GetPoint(6), this->GetPoint(2), this->GetPoint(3), point);
1115  distances[4] = GeometryUtils::PointDistanceToQuadrilateral3D(this->GetPoint(7), this->GetPoint(3), this->GetPoint(0), this->GetPoint(4), point);
1116  distances[5] = GeometryUtils::PointDistanceToQuadrilateral3D(this->GetPoint(4), this->GetPoint(5), this->GetPoint(6), this->GetPoint(7), point);
1117  const auto min = std::min_element(distances.begin(), distances.end());
1118  return *min;
1119  }
1120 
1124 
1132  std::string Info() const override
1133  {
1134  return "3 dimensional hexahedra with eight nodes in 3D space";
1135  }
1136 
1144  void PrintInfo( std::ostream& rOStream ) const override
1145  {
1146  rOStream << "3 dimensional hexahedra with eight nodes in 3D space";
1147  }
1148 
1158  void PrintData( std::ostream& rOStream ) const override
1159  {
1160  // Base Geometry class PrintData call
1161  BaseType::PrintData( rOStream );
1162  std::cout << std::endl;
1163 
1164  // If the geometry has valid points, calculate and output its data
1165  if (this->AllPointsAreValid()) {
1166  Matrix jacobian;
1167  this->Jacobian( jacobian, PointType() );
1168  rOStream << " Jacobian in the origin\t : " << jacobian;
1169  }
1170  }
1171 
1172 protected:
1173 
1178 private:
1179 
1184  static const GeometryData msGeometryData;
1185 
1186  static const GeometryDimension msGeometryDimension;
1187 
1191 
1192  friend class Serializer;
1193 
1194  void save( Serializer& rSerializer ) const override
1195  {
1197  }
1198 
1199  void load( Serializer& rSerializer ) override
1200  {
1202  }
1203 
1204  Hexahedra3D8(): BaseType( PointsArrayType(), &msGeometryData ) {}
1205 
1222  {
1223  Matrix result = ZeroMatrix( 8, 3 );
1224  result( 0, 0 ) = -0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 - rPoint[2] );
1225  result( 0, 1 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[2] );
1226  result( 0, 2 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[1] );
1227  result( 1, 0 ) = 0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 - rPoint[2] );
1228  result( 1, 1 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[2] );
1229  result( 1, 2 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[1] );
1230  result( 2, 0 ) = 0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 - rPoint[2] );
1231  result( 2, 1 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[2] );
1232  result( 2, 2 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[1] );
1233  result( 3, 0 ) = -0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 - rPoint[2] );
1234  result( 3, 1 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[2] );
1235  result( 3, 2 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[1] );
1236  result( 4, 0 ) = -0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 + rPoint[2] );
1237  result( 4, 1 ) = -0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[2] );
1238  result( 4, 2 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 - rPoint[1] );
1239  result( 5, 0 ) = 0.125 * ( 1.0 - rPoint[1] ) * ( 1.0 + rPoint[2] );
1240  result( 5, 1 ) = -0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[2] );
1241  result( 5, 2 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 - rPoint[1] );
1242  result( 6, 0 ) = 0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 + rPoint[2] );
1243  result( 6, 1 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[2] );
1244  result( 6, 2 ) = 0.125 * ( 1.0 + rPoint[0] ) * ( 1.0 + rPoint[1] );
1245  result( 7, 0 ) = -0.125 * ( 1.0 + rPoint[1] ) * ( 1.0 + rPoint[2] );
1246  result( 7, 1 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[2] );
1247  result( 7, 2 ) = 0.125 * ( 1.0 - rPoint[0] ) * ( 1.0 + rPoint[1] );
1248  return result;
1249  }
1250 
1251 
1263  static Matrix CalculateShapeFunctionsIntegrationPointsValues(
1264  typename BaseType::IntegrationMethod ThisMethod )
1265  {
1266  IntegrationPointsContainerType all_integration_points = AllIntegrationPoints();
1267  IntegrationPointsArrayType& integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1268 
1269  //number of integration points
1270  const int integration_points_number = integration_points.size();
1271  //number of nodes in current geometry
1272  const int points_number = 8;
1273  //setting up return matrix
1274  Matrix shape_function_values( integration_points_number, points_number );
1275  //loop over all integration points
1276 
1277  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1278  {
1279  shape_function_values( pnt, 0 ) =
1280  0.125 * ( 1.0 - integration_points[pnt].X() )
1281  * ( 1.0 - integration_points[pnt].Y() )
1282  * ( 1.0 - integration_points[pnt].Z() );
1283  shape_function_values( pnt, 1 ) =
1284  0.125 * ( 1.0 + integration_points[pnt].X() )
1285  * ( 1.0 - integration_points[pnt].Y() )
1286  * ( 1.0 - integration_points[pnt].Z() );
1287  shape_function_values( pnt, 2 ) =
1288  0.125 * ( 1.0 + integration_points[pnt].X() )
1289  * ( 1.0 + integration_points[pnt].Y() )
1290  * ( 1.0 - integration_points[pnt].Z() );
1291  shape_function_values( pnt, 3 ) =
1292  0.125 * ( 1.0 - integration_points[pnt].X() )
1293  * ( 1.0 + integration_points[pnt].Y() )
1294  * ( 1.0 - integration_points[pnt].Z() );
1295  shape_function_values( pnt, 4 ) =
1296  0.125 * ( 1.0 - integration_points[pnt].X() )
1297  * ( 1.0 - integration_points[pnt].Y() )
1298  * ( 1.0 + integration_points[pnt].Z() );
1299  shape_function_values( pnt, 5 ) =
1300  0.125 * ( 1.0 + integration_points[pnt].X() )
1301  * ( 1.0 - integration_points[pnt].Y() )
1302  * ( 1.0 + integration_points[pnt].Z() );
1303  shape_function_values( pnt, 6 ) =
1304  0.125 * ( 1.0 + integration_points[pnt].X() )
1305  * ( 1.0 + integration_points[pnt].Y() )
1306  * ( 1.0 + integration_points[pnt].Z() );
1307  shape_function_values( pnt, 7 ) =
1308  0.125 * ( 1.0 - integration_points[pnt].X() )
1309  * ( 1.0 + integration_points[pnt].Y() )
1310  * ( 1.0 + integration_points[pnt].Z() );
1311  }
1312 
1313  return shape_function_values;
1314  }
1315 
1329  CalculateShapeFunctionsIntegrationPointsLocalGradients(
1330  typename BaseType::IntegrationMethod ThisMethod )
1331  {
1332  IntegrationPointsContainerType all_integration_points =
1333  AllIntegrationPoints();
1334  IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
1335  //number of integration points
1336  const int integration_points_number = integration_points.size();
1337  ShapeFunctionsGradientsType d_shape_f_values( integration_points_number );
1338  //initialising container
1339  //loop over all integration points
1340 
1341  for ( int pnt = 0; pnt < integration_points_number; pnt++ )
1342  {
1343  Matrix& result = d_shape_f_values[pnt];
1344  result = ZeroMatrix( 8, 3 );
1345  result( 0, 0 ) =
1346  -0.125 * ( 1.0 - integration_points[pnt].Y() )
1347  * ( 1.0 - integration_points[pnt].Z() );
1348  result( 0, 1 ) =
1349  -0.125 * ( 1.0 - integration_points[pnt].X() )
1350  * ( 1.0 - integration_points[pnt].Z() );
1351  result( 0, 2 ) =
1352  -0.125 * ( 1.0 - integration_points[pnt].X() )
1353  * ( 1.0 - integration_points[pnt].Y() );
1354  result( 1, 0 ) =
1355  0.125 * ( 1.0 - integration_points[pnt].Y() )
1356  * ( 1.0 - integration_points[pnt].Z() );
1357  result( 1, 1 ) =
1358  -0.125 * ( 1.0 + integration_points[pnt].X() )
1359  * ( 1.0 - integration_points[pnt].Z() );
1360  result( 1, 2 ) =
1361  -0.125 * ( 1.0 + integration_points[pnt].X() )
1362  * ( 1.0 - integration_points[pnt].Y() );
1363  result( 2, 0 ) =
1364  0.125 * ( 1.0 + integration_points[pnt].Y() )
1365  * ( 1.0 - integration_points[pnt].Z() );
1366  result( 2, 1 ) =
1367  0.125 * ( 1.0 + integration_points[pnt].X() )
1368  * ( 1.0 - integration_points[pnt].Z() );
1369  result( 2, 2 ) =
1370  -0.125 * ( 1.0 + integration_points[pnt].X() )
1371  * ( 1.0 + integration_points[pnt].Y() );
1372  result( 3, 0 ) =
1373  -0.125 * ( 1.0 + integration_points[pnt].Y() )
1374  * ( 1.0 - integration_points[pnt].Z() );
1375  result( 3, 1 ) =
1376  0.125 * ( 1.0 - integration_points[pnt].X() )
1377  * ( 1.0 - integration_points[pnt].Z() );
1378  result( 3, 2 ) =
1379  -0.125 * ( 1.0 - integration_points[pnt].X() )
1380  * ( 1.0 + integration_points[pnt].Y() );
1381  result( 4, 0 ) =
1382  -0.125 * ( 1.0 - integration_points[pnt].Y() )
1383  * ( 1.0 + integration_points[pnt].Z() );
1384  result( 4, 1 ) =
1385  -0.125 * ( 1.0 - integration_points[pnt].X() )
1386  * ( 1.0 + integration_points[pnt].Z() );
1387  result( 4, 2 ) =
1388  0.125 * ( 1.0 - integration_points[pnt].X() )
1389  * ( 1.0 - integration_points[pnt].Y() );
1390  result( 5, 0 ) =
1391  0.125 * ( 1.0 - integration_points[pnt].Y() )
1392  * ( 1.0 + integration_points[pnt].Z() );
1393  result( 5, 1 ) =
1394  -0.125 * ( 1.0 + integration_points[pnt].X() )
1395  * ( 1.0 + integration_points[pnt].Z() );
1396  result( 5, 2 ) =
1397  0.125 * ( 1.0 + integration_points[pnt].X() )
1398  * ( 1.0 - integration_points[pnt].Y() );
1399  result( 6, 0 ) =
1400  0.125 * ( 1.0 + integration_points[pnt].Y() )
1401  * ( 1.0 + integration_points[pnt].Z() );
1402  result( 6, 1 ) =
1403  0.125 * ( 1.0 + integration_points[pnt].X() )
1404  * ( 1.0 + integration_points[pnt].Z() );
1405  result( 6, 2 ) =
1406  0.125 * ( 1.0 + integration_points[pnt].X() )
1407  * ( 1.0 + integration_points[pnt].Y() );
1408  result( 7, 0 ) =
1409  -0.125 * ( 1.0 + integration_points[pnt].Y() )
1410  * ( 1.0 + integration_points[pnt].Z() );
1411  result( 7, 1 ) =
1412  0.125 * ( 1.0 - integration_points[pnt].X() )
1413  * ( 1.0 + integration_points[pnt].Z() );
1414  result( 7, 2 ) =
1415  0.125 * ( 1.0 - integration_points[pnt].X() )
1416  * ( 1.0 + integration_points[pnt].Y() );
1417  }
1418 
1419  return d_shape_f_values;
1420  }
1421 
1422  static const IntegrationPointsContainerType AllIntegrationPoints()
1423  {
1424  IntegrationPointsContainerType integration_points =
1425  {
1426  {
1427  Quadrature < HexahedronGaussLegendreIntegrationPoints1, 3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1428  Quadrature < HexahedronGaussLegendreIntegrationPoints2, 3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1429  Quadrature < HexahedronGaussLegendreIntegrationPoints3, 3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1430  Quadrature < HexahedronGaussLegendreIntegrationPoints4, 3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1431  Quadrature < HexahedronGaussLegendreIntegrationPoints5, 3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1432  Quadrature < HexahedronGaussLobattoIntegrationPoints1, 3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
1433  Quadrature < HexahedronGaussLobattoIntegrationPoints2, 3, IntegrationPoint<3> >::GenerateIntegrationPoints()
1434  }
1435  };
1436  return integration_points;
1437  }
1438 
1439  static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
1440  {
1441  ShapeFunctionsValuesContainerType shape_functions_values =
1442  {
1443  {
1444  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
1445  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
1446  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_3 ),
1447  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_4 ),
1448  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_5 ),
1449  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_1 ),
1450  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_2 )
1451  }
1452  };
1453  return shape_functions_values;
1454  }
1455 
1460  AllShapeFunctionsLocalGradients()
1461  {
1462  ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
1463  {
1464  {
1465  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
1466  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
1467  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_3 ),
1468  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_4 ),
1469  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_5 ),
1470  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_1 ),
1471  Hexahedra3D8<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_2 )
1472  }
1473  };
1474  return shape_functions_local_gradients;
1475  }
1476 
1477 
1482  template<class TOtherPointType> friend class Hexahedra3D8;
1483 
1484 
1489 };// Class Hexahedra3D8
1490 
1491 
1499 template<class TPointType> inline std::istream& operator >> (
1500  std::istream& rIStream, Hexahedra3D8<TPointType>& rThis );
1501 
1505 template<class TPointType> inline std::ostream& operator << (
1506  std::ostream& rOStream, const Hexahedra3D8<TPointType>& rThis )
1507 {
1508  rThis.PrintInfo( rOStream );
1509  rOStream << std::endl;
1510  rThis.PrintData( rOStream );
1511 
1512  return rOStream;
1513 }
1514 
1515 
1516 template<class TPointType> const
1517 GeometryData Hexahedra3D8<TPointType>::msGeometryData(
1520  Hexahedra3D8<TPointType>::AllIntegrationPoints(),
1521  Hexahedra3D8<TPointType>::AllShapeFunctionsValues(),
1522  AllShapeFunctionsLocalGradients()
1523 );
1524 
1525 template<class TPointType> const
1526 GeometryDimension Hexahedra3D8<TPointType>::msGeometryDimension(3, 3);
1527 
1528 }// namespace Kratos.
Definition: geometry_data.h:60
KratosGeometryType
Definition: geometry_data.h:110
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
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
Vector & DeterminantOfJacobian(Vector &rResult) const
Definition: geometry.h:3154
virtual CoordinatesArrayType & PointLocalCoordinates(CoordinatesArrayType &rResult, const CoordinatesArrayType &rPoint) const
Returns the local coordinates of a given arbitrary point.
Definition: geometry.h:1854
const PointsArrayType & Points() const
Definition: geometry.h:1768
bool AllPointsAreValid() const
Checks if the geometry points are valid Checks if the geometry points are valid from the pointer valu...
Definition: geometry.h:4093
std::array< IntegrationPointsArrayType, static_cast< int >GeometryData::IntegrationMethod::NumberOfIntegrationMethods)> IntegrationPointsContainerType
Definition: geometry.h:167
const ShapeFunctionsGradientsType & ShapeFunctionsLocalGradients() const
Definition: geometry.h:3539
JacobiansType & Jacobian(JacobiansType &rResult) const
Definition: geometry.h:2901
TPointType const & GetPoint(const int Index) const
Definition: geometry.h:1816
GeometryData::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: geometry.h:177
GeometryData::IntegrationMethod IntegrationMethod
Definition: geometry.h:122
static double PointDistanceToQuadrilateral3D(const Point &rQuadrilateralPoint1, const Point &rQuadrilateralPoint2, const Point &rQuadrilateralPoint3, const Point &rQuadrilateralPoint4, const Point &rPoint)
This function calculates the distance of a 3D point to a 3D quadrilateral.
Definition: geometry_utilities.cpp:376
An eight node hexahedra geometry with linear shape functions.
Definition: hexahedra_3d_8.h:55
GeometryData::KratosGeometryFamily GetGeometryFamily() const override
Definition: hexahedra_3d_8.h:285
double MinEdgeLength() const override
Definition: hexahedra_3d_8.h:632
Line3D2< TPointType > EdgeType
Definition: hexahedra_3d_8.h:69
double ShortestToLongestEdgeQuality() const override
Calculates the shortest to longest edge quality metric. Calculates the shortest to longest edge quali...
Definition: hexahedra_3d_8.h:853
double Length() const override
Definition: hexahedra_3d_8.h:389
Matrix MatrixType
Definition: hexahedra_3d_8.h:184
Geometry< TPointType > BaseType
Definition: hexahedra_3d_8.h:64
Quadrilateral3D4< TPointType > FaceType
Definition: hexahedra_3d_8.h:70
BaseType::CoordinatesArrayType CoordinatesArrayType
Definition: hexahedra_3d_8.h:179
double Area() const override
Definition: hexahedra_3d_8.h:408
void ComputeDihedralAngles(Vector &rDihedralAngles) const override
Implements the calculus of the 24 dihedral angles of the hexa.
Definition: hexahedra_3d_8.h:756
Hexahedra3D8(Hexahedra3D8 const &rOther)
Definition: hexahedra_3d_8.h:260
Hexahedra3D8(const IndexType GeometryId, const PointsArrayType &rThisPoints)
Constructor with Geometry Id.
Definition: hexahedra_3d_8.h:234
double Volume() const override
This method calculate and return volume of this geometry.
Definition: hexahedra_3d_8.h:421
BaseType::IntegrationPointsContainerType IntegrationPointsContainerType
Definition: hexahedra_3d_8.h:134
friend class Hexahedra3D8
Definition: hexahedra_3d_8.h:1482
BaseType::ShapeFunctionsSecondDerivativesType ShapeFunctionsSecondDerivativesType
Definition: hexahedra_3d_8.h:169
BaseType::NormalType NormalType
Definition: hexahedra_3d_8.h:174
void PrintInfo(std::ostream &rOStream) const override
Definition: hexahedra_3d_8.h:1144
BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType
Definition: hexahedra_3d_8.h:141
Hexahedra3D8(const PointType &Point1, const PointType &Point2, const PointType &Point3, const PointType &Point4, const PointType &Point5, const PointType &Point6, const PointType &Point7, const PointType &Point8)
Definition: hexahedra_3d_8.h:190
double MaxDihedralAngle() const override
Calculates the max dihedral angle quality metric.
Definition: hexahedra_3d_8.h:808
TPointType PointType
Definition: hexahedra_3d_8.h:91
BaseType::JacobiansType JacobiansType
Definition: hexahedra_3d_8.h:155
double MinDihedralAngle() const override
Calculates the min dihedral angle quality metric.
Definition: hexahedra_3d_8.h:791
double VolumeToRMSEdgeLength() const override
Calculates the volume to average edge length quality metric.
Definition: hexahedra_3d_8.h:829
Hexahedra3D8(Hexahedra3D8< TOtherPointType > const &rOther)
Definition: hexahedra_3d_8.h:277
double AverageEdgeLength() const override
Definition: hexahedra_3d_8.h:583
BaseType::Pointer Create(const IndexType NewGeometryId, PointsArrayType const &rThisPoints) const override
Creates a new geometry pointer.
Definition: hexahedra_3d_8.h:346
BaseType::IntegrationPointType IntegrationPointType
Definition: hexahedra_3d_8.h:119
BaseType::SizeType SizeType
Definition: hexahedra_3d_8.h:106
BaseType::IntegrationPointsArrayType IntegrationPointsArrayType
Definition: hexahedra_3d_8.h:127
BaseType::Pointer Create(const IndexType NewGeometryId, const BaseType &rGeometry) const override
Creates a new geometry pointer.
Definition: hexahedra_3d_8.h:360
ShapeFunctionsSecondDerivativesType & ShapeFunctionsSecondDerivatives(ShapeFunctionsSecondDerivativesType &rResult, const CoordinatesArrayType &rPoint) const override
Definition: hexahedra_3d_8.h:981
BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType
Definition: hexahedra_3d_8.h:148
Vector & ShapeFunctionsValues(Vector &rResult, const CoordinatesArrayType &rCoordinates) const override
Definition: hexahedra_3d_8.h:917
Hexahedra3D8(typename PointType::Pointer pPoint1, typename PointType::Pointer pPoint2, typename PointType::Pointer pPoint3, typename PointType::Pointer pPoint4, typename PointType::Pointer pPoint5, typename PointType::Pointer pPoint6, typename PointType::Pointer pPoint7, typename PointType::Pointer pPoint8)
Definition: hexahedra_3d_8.h:206
Matrix & ShapeFunctionsLocalGradients(Matrix &rResult, const CoordinatesArrayType &rPoint) const override
Definition: hexahedra_3d_8.h:939
BaseType::PointsArrayType PointsArrayType
Definition: hexahedra_3d_8.h:112
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: hexahedra_3d_8.h:497
double DomainSize() const override
Definition: hexahedra_3d_8.h:439
GeometriesArrayType GenerateFaces() const override
Returns all faces of the current geometry.
Definition: hexahedra_3d_8.h:665
double MaxEdgeLength() const override
Definition: hexahedra_3d_8.h:615
Hexahedra3D8(const PointsArrayType &ThisPoints)
Definition: hexahedra_3d_8.h:226
Hexahedra3D8 & operator=(Hexahedra3D8< TOtherPointType > const &rOther)
Definition: hexahedra_3d_8.h:328
GeometryData::IntegrationMethod IntegrationMethod
Definition: hexahedra_3d_8.h:80
std::string Info() const override
Definition: hexahedra_3d_8.h:1132
BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: hexahedra_3d_8.h:162
SizeType FacesNumber() const override
Returns the number of faces of the current geometry.
Definition: hexahedra_3d_8.h:652
void ComputeSolidAngles(Vector &rSolidAngles) const override
Implements the calculus of the 8 solid angles of the hexa.
Definition: hexahedra_3d_8.h:733
bool HasIntersection(const Point &rLowPoint, const Point &rHighPoint) const override
Definition: hexahedra_3d_8.h:703
GeometryData::KratosGeometryType GetGeometryType() const override
Definition: hexahedra_3d_8.h:290
void PrintData(std::ostream &rOStream) const override
Definition: hexahedra_3d_8.h:1158
GeometriesArrayType GenerateEdges() const override
Definition: hexahedra_3d_8.h:535
Matrix & PointsLocalCoordinates(Matrix &rResult) const override
Definition: hexahedra_3d_8.h:449
SizeType EdgesNumber() const override
Definition: hexahedra_3d_8.h:524
Hexahedra3D8(const std::string &rGeometryName, const PointsArrayType &rThisPoints)
Constructor with Geometry Name.
Definition: hexahedra_3d_8.h:243
~Hexahedra3D8() override
Destructor. Does nothing!!!
Definition: hexahedra_3d_8.h:283
Hexahedra3D8 & operator=(const Hexahedra3D8 &rOther)
Definition: hexahedra_3d_8.h:310
BaseType::GeometriesArrayType GeometriesArrayType
Definition: hexahedra_3d_8.h:86
KRATOS_CLASS_POINTER_DEFINITION(Hexahedra3D8)
double CalculateDistance(const CoordinatesArrayType &rPointGlobalCoordinates, const double Tolerance=std::numeric_limits< double >::epsilon()) const override
Computes the distance between an point in global coordinates and the closest point of this geometry....
Definition: hexahedra_3d_8.h:1095
BaseType::IndexType IndexType
Definition: hexahedra_3d_8.h:98
double ShapeFunctionValue(IndexType ShapeFunctionIndex, const CoordinatesArrayType &rPoint) const override
Definition: hexahedra_3d_8.h:878
Short class definition.
Definition: integration_point.h:52
static double ComputeVolume3DGeometry(const Geometry< TPointType > &rGeometry)
This method calculates and returns the volume of the geometry from a 3D geometry.
Definition: integration_utilities.h:168
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 3D line geometry with linear shape functions.
Definition: line_3d_2.h:64
Various mathematical utilitiy functions.
Definition: math_utils.h:62
This class defines the node.
Definition: node.h:65
Point class.
Definition: point.h:59
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 four node 3D quadrilateral geometry with bi-linear shape functions.
Definition: quadrilateral_3d_4.h:76
The serialization consists in storing the state of an object into a storage format like data file or ...
Definition: serializer.h:123
Short class definition.
Definition: array_1d.h:61
#define KRATOS_SERIALIZE_SAVE_BASE_CLASS(Serializer, BaseType)
Definition: define.h:812
#define KRATOS_SERIALIZE_LOAD_BASE_CLASS(Serializer, BaseType)
Definition: define.h:815
#define KRATOS_ERROR
Definition: exception.h:161
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
static double max(double a, double b)
Definition: GeometryFunctions.h:79
static double min(double a, double b)
Definition: GeometryFunctions.h:71
constexpr double Pi
Definition: global_variables.h:25
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
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
const GeometryData Hexahedra3D8< TPointType >::msGeometryData & msGeometryDimension(), Hexahedra3D8< TPointType >::AllShapeFunctionsValues(), AllShapeFunctionsLocalGradients()
Definition: brep_curve.h:483
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
def load(f)
Definition: ode_solve.py:307
float temp
Definition: rotating_cone.py:85
integer i
Definition: TensorModule.f:17