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.
coordinate_transformation_utilities.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:
11 //
12 //
13 
14 
15 #ifndef KRATOS_COORDINATE_TRANSFORMATION_UTILITIES_H
16 #define KRATOS_COORDINATE_TRANSFORMATION_UTILITIES_H
17 
18 // system includes
19 
20 // external includes
21 #include "boost/numeric/ublas/matrix_proxy.hpp"
22 
23 // kratos includes
24 #include "includes/define.h"
25 #include "includes/node.h"
26 #include "includes/model_part.h"
27 #include "containers/variable.h"
28 #include "geometries/geometry.h"
29 
30 namespace Kratos {
31 
34 
37 
41 
45 
49 
53 
56 template<class TLocalMatrixType, class TLocalVectorType, class TValueType>
57 class CoordinateTransformationUtils {
58 public:
61 
63  KRATOS_CLASS_POINTER_DEFINITION(CoordinateTransformationUtils);
64 
65  typedef Node NodeType;
66 
67  typedef Geometry< Node > GeometryType;
68 
69 // typedef boost::numeric::ublas::matrix_row<TLocalMatrixType> LocalRowType;
70 //
71 // typedef boost::numeric::ublas::matrix_range<TLocalMatrixType> MatrixBlockType;
72 
76 
78 
82  CoordinateTransformationUtils(const unsigned int DomainSize,
83  const unsigned int NumRowsPerNode,
84  const Kratos::Flags& rSelectionFlag = SLIP):
85  mDomainSize(DomainSize),
86  mBlockSize(NumRowsPerNode),
87  mrFlag(rSelectionFlag)
88  {}
89 
91  virtual ~CoordinateTransformationUtils() {}
92 
96 
100 
110  virtual void CalculateRotationOperatorPure(
111  TLocalMatrixType& rRotationMatrix,
112  const GeometryType::PointType& rThisPoint) const
113  {
114  KRATOS_TRY
115 
116  if (mDomainSize == 2) {
117  BoundedMatrix<double, 2, 2> local_matrix;
118  this->LocalRotationOperatorPure(local_matrix, rThisPoint);
119  if (rRotationMatrix.size1() != 2 || rRotationMatrix.size2() != 2) {
120  rRotationMatrix.resize(2, 2, false);
121  }
122  noalias(rRotationMatrix) = local_matrix;
123  } else if (mDomainSize == 3) {
124  BoundedMatrix<double, 3, 3> local_matrix;
125  this->LocalRotationOperatorPure(local_matrix, rThisPoint);
126  if (rRotationMatrix.size1() != 3 || rRotationMatrix.size2() != 3) {
127  rRotationMatrix.resize(3, 3, false);
128  }
129  noalias(rRotationMatrix) = local_matrix;
130  } else {
131  KRATOS_ERROR << "Unsupported domain size [ mDomainSize = " << mDomainSize
132  << " ].\n";
133  }
134 
135  KRATOS_CATCH("");
136  }
137 
138  void LocalRotationOperatorPure(
139  BoundedMatrix<double, 3, 3>& rRot,
140  const GeometryType::PointType& rThisPoint) const
141  {
142  // Get the normal evaluated at the node
143  const array_1d<double, 3>& rNormal = rThisPoint.FastGetSolutionStepValue(NORMAL);
144 
145  double aux = rNormal[0] * rNormal[0] + rNormal[1] * rNormal[1] +
146  rNormal[2] * rNormal[2];
147  aux = sqrt(aux);
148  rRot(0, 0) = rNormal[0] / aux;
149  rRot(0, 1) = rNormal[1] / aux;
150  rRot(0, 2) = rNormal[2] / aux;
151  // Define the new coordinate system, where the first vector is aligned with the normal
152 
153  // To choose the remaining two vectors, we project the first component of the cartesian base to the tangent plane
154  array_1d<double, 3> rT1;
155  rT1(0) = 1.0;
156  rT1(1) = 0.0;
157  rT1(2) = 0.0;
158  double dot = rRot(0, 0); // this->Dot(rN,rT1);
159 
160  // It is possible that the normal is aligned with (1,0,0), resulting in
161  // norm(rT1) = 0 If this is the case, repeat the procedure using (0,1,0)
162  if (fabs(dot) > 0.99) {
163  rT1(0) = 0.0;
164  rT1(1) = 1.0;
165  rT1(2) = 0.0;
166 
167  dot = rRot(0, 1); // this->Dot(rN,rT1);
168  }
169 
170  // calculate projection and normalize
171  rT1[0] -= dot * rRot(0, 0);
172  rT1[1] -= dot * rRot(0, 1);
173  rT1[2] -= dot * rRot(0, 2);
174  Normalize(rT1);
175  rRot(1, 0) = rT1[0];
176  rRot(1, 1) = rT1[1];
177  rRot(1, 2) = rT1[2];
178 
179  // The third base component is choosen as N x T1, which is normalized by construction
180  rRot(2, 0) = rRot(0, 1) * rT1[2] - rRot(0, 2) * rT1[1];
181  rRot(2, 1) = rRot(0, 2) * rT1[0] - rRot(0, 0) * rT1[2];
182  rRot(2, 2) = rRot(0, 0) * rT1[1] - rRot(0, 1) * rT1[0];
183  }
184 
185  void LocalRotationOperatorPure(
186  BoundedMatrix<double, 2, 2>& rRot,
187  const GeometryType::PointType& rThisPoint) const
188  {
189  // Get the normal evaluated at the node
190  const array_1d<double, 3>& rNormal = rThisPoint.FastGetSolutionStepValue(NORMAL);
191 
192  double aux = rNormal[0] * rNormal[0] + rNormal[1] * rNormal[1];
193  aux = sqrt(aux);
194 
195  rRot(0, 0) = rNormal[0] / aux;
196  rRot(0, 1) = rNormal[1] / aux;
197  rRot(1, 0) = -rNormal[1] / aux;
198  rRot(1, 1) = rNormal[0] / aux;
199  }
200 
216  virtual void CalculateRotationOperatorPureShapeSensitivities(
217  TLocalMatrixType& rRotationMatrixShapeDerivative,
218  const std::size_t DerivativeNodeIndex,
219  const std::size_t DerivativeDirectionIndex,
220  const GeometryType::PointType& rThisPoint) const
221  {
222  KRATOS_TRY
223 
224  if (mDomainSize == 2) {
225  BoundedMatrix<double, 2, 2> local_matrix;
226  this->CalculateRotationOperatorPureShapeSensitivities(
227  local_matrix, DerivativeNodeIndex, DerivativeDirectionIndex, rThisPoint);
228  if (rRotationMatrixShapeDerivative.size1() != 2 ||
229  rRotationMatrixShapeDerivative.size2() != 2) {
230  rRotationMatrixShapeDerivative.resize(2, 2, false);
231  }
232  noalias(rRotationMatrixShapeDerivative) = local_matrix;
233  } else if (mDomainSize == 3) {
234  BoundedMatrix<double, 3, 3> local_matrix;
235  this->CalculateRotationOperatorPureShapeSensitivities(
236  local_matrix, DerivativeNodeIndex, DerivativeDirectionIndex, rThisPoint);
237  if (rRotationMatrixShapeDerivative.size1() != 3 ||
238  rRotationMatrixShapeDerivative.size2() != 3) {
239  rRotationMatrixShapeDerivative.resize(3, 3, false);
240  }
241  noalias(rRotationMatrixShapeDerivative) = local_matrix;
242  } else {
243  KRATOS_ERROR << "Unsupported domain size [ mDomainSize = " << mDomainSize
244  << " ].\n";
245  }
246 
247 
248  KRATOS_CATCH("");
249  }
250 
267  virtual void CalculateRotationOperatorPureShapeSensitivities(
268  BoundedMatrix<double, 2, 2>& rOutput,
269  const std::size_t DerivativeNodeIndex,
270  const std::size_t DerivativeDirectionIndex,
271  const GeometryType::PointType& rThisPoint) const
272  {
273  KRATOS_TRY
274 
275  KRATOS_ERROR_IF(!rThisPoint.SolutionStepsDataHas(NORMAL))
276  << "NORMAL is not found in node at " << rThisPoint.Coordinates() << ".";
277  KRATOS_ERROR_IF(!rThisPoint.Has(NORMAL_SHAPE_DERIVATIVE))
278  << "NORMAL_SHAPE_DERIVATIVE is not found in node [ Node.Id() = "
279  << rThisPoint.Id() << " ] at " << rThisPoint.Coordinates() << ".";
280 
281  const array_1d<double, 3>& r_nodal_normal =
282  rThisPoint.FastGetSolutionStepValue(NORMAL);
283  const double nodal_normal_magnitude = norm_2(r_nodal_normal);
284 
285  KRATOS_ERROR_IF(nodal_normal_magnitude == 0.0)
286  << "NORMAL at node " << rThisPoint.Coordinates()
287  << " is not properly initialized.";
288 
289  const Matrix& r_sensitivity_values = rThisPoint.GetValue(NORMAL_SHAPE_DERIVATIVE);
290 
291  KRATOS_DEBUG_ERROR_IF(r_sensitivity_values.size2() != 2)
292  << "NORMAL_SHAPE_DERIVATIVE is not properly initialized at node [ Node.Id() = "
293  << rThisPoint.Id() << " ] "
294  << rThisPoint.Coordinates() << " to calculate 2D rotation operator shape sensitivities. [ required number of columns = 2, available number of columns = "
295  << r_sensitivity_values.size2() << " ].";
296 
297  const std::size_t require_rows = (DerivativeNodeIndex + 1) * 2;
298  KRATOS_DEBUG_ERROR_IF(r_sensitivity_values.size1() < require_rows)
299  << "NORMAL_SHAPE_DERIVATIVE is not properly initialized at node [ Node.Id() = "
300  << rThisPoint.Id() << " ] "
301  << rThisPoint.Coordinates() << " to calculate 2D rotation operator shape sensitivities. [ required number of rows >= "
302  << require_rows
303  << ", available number of rows = " << r_sensitivity_values.size1() << " ].";
304 
305  const Vector& r_nodal_normal_derivatives =
306  row(r_sensitivity_values, DerivativeNodeIndex * 2 + DerivativeDirectionIndex);
307 
308  rOutput(0, 0) = r_nodal_normal_derivatives[0] / nodal_normal_magnitude;
309  rOutput(0, 1) = r_nodal_normal_derivatives[1] / nodal_normal_magnitude;
310  rOutput(1, 0) = -r_nodal_normal_derivatives[1] / nodal_normal_magnitude;
311  rOutput(1, 1) = r_nodal_normal_derivatives[0] / nodal_normal_magnitude;
312 
313  const double nodal_normal_magnitude_derivative =
314  (r_nodal_normal[0] * r_nodal_normal_derivatives[0] +
315  r_nodal_normal[1] * r_nodal_normal_derivatives[1]) /
316  nodal_normal_magnitude;
317 
318  const double coeff = nodal_normal_magnitude_derivative /
319  (std::pow(nodal_normal_magnitude, 2));
320 
321  rOutput(0, 0) -= r_nodal_normal[0] * coeff;
322  rOutput(0, 1) -= r_nodal_normal[1] * coeff;
323  rOutput(1, 0) -= -r_nodal_normal[1] * coeff;
324  rOutput(1, 1) -= r_nodal_normal[0] * coeff;
325 
326  KRATOS_CATCH("");
327  }
328 
345  virtual void CalculateRotationOperatorPureShapeSensitivities(
346  BoundedMatrix<double, 3, 3>& rOutput,
347  const std::size_t DerivativeNodeIndex,
348  const std::size_t DerivativeDirectionIndex,
349  const GeometryType::PointType& rThisPoint) const
350  {
351  KRATOS_TRY
352 
353  KRATOS_ERROR_IF(!rThisPoint.SolutionStepsDataHas(NORMAL))
354  << "NORMAL is not found in node at " << rThisPoint.Coordinates() << ".";
355  KRATOS_ERROR_IF(!rThisPoint.Has(NORMAL_SHAPE_DERIVATIVE))
356  << "NORMAL_SHAPE_DERIVATIVE is not found in node at "
357  << rThisPoint.Coordinates() << ".";
358 
359  const array_1d<double, 3>& r_nodal_normal =
360  rThisPoint.FastGetSolutionStepValue(NORMAL);
361  const double nodal_normal_magnitude = norm_2(r_nodal_normal);
362 
363  KRATOS_ERROR_IF(nodal_normal_magnitude == 0.0)
364  << "NORMAL at node " << rThisPoint.Coordinates()
365  << " is not properly initialized.";
366 
367  const Matrix& r_sensitivity_values = rThisPoint.GetValue(NORMAL_SHAPE_DERIVATIVE);
368 
369  KRATOS_DEBUG_ERROR_IF(r_sensitivity_values.size2() != 3)
370  << "NORMAL_SHAPE_DERIVATIVE is not properly initialized at node "
371  << rThisPoint.Coordinates() << " to calculate 3D rotation operator shape sensitivities. [ required number of columns = 3, available number of columns = "
372  << r_sensitivity_values.size2() << " ].";
373 
374  const std::size_t require_rows = (DerivativeNodeIndex + 1) * 3;
375  KRATOS_DEBUG_ERROR_IF(r_sensitivity_values.size1() < require_rows)
376  << "NORMAL_SHAPE_DERIVATIVE is not properly initialized at node "
377  << rThisPoint.Coordinates() << " to calculate 3D rotation operator shape sensitivities. [ required number of rows >= "
378  << require_rows
379  << ", available number of rows = " << r_sensitivity_values.size1() << " ].";
380 
381  const Vector& r_nodal_normal_derivative =
382  row(r_sensitivity_values, DerivativeNodeIndex * 3 + DerivativeDirectionIndex);
383 
384  const double nodal_normal_magnitude_derivative = VectorNormDerivative(nodal_normal_magnitude, r_nodal_normal, r_nodal_normal_derivative);
385  const array_1d<double, 3>& unit_normal = r_nodal_normal / nodal_normal_magnitude;
386  const array_1d<double, 3>& unit_normal_derivative = UnitVectorDerivative(nodal_normal_magnitude, nodal_normal_magnitude_derivative, r_nodal_normal, r_nodal_normal_derivative);
387 
388  rOutput(0, 0) = unit_normal_derivative[0];
389  rOutput(0, 1) = unit_normal_derivative[1];
390  rOutput(0, 2) = unit_normal_derivative[2];
391 
392  array_1d<double, 3> rT1(3, 0.0);
393  rT1[0] = 1.0;
394  double dot = unit_normal[0];
395  double dot_derivative = unit_normal_derivative[0];
396 
397  if (std::abs(dot) > 0.99) {
398  rT1[0] = 0.0;
399  rT1[1] = 1.0;
400  dot = unit_normal[1];
401  dot_derivative = unit_normal_derivative[1];
402  }
403 
404  // calculate rT1
405  noalias(rT1) -= unit_normal * dot;
406  const double rT1_norm = norm_2(rT1);
407  const array_1d<double, 3>& unit_rT1 = rT1 / rT1_norm;
408 
409  // calculate rT1 derivative
410  const array_1d<double, 3>& rT1_derivative = (unit_normal_derivative * dot + unit_normal * dot_derivative) * -1.0;
411 
412  // calculate rT1 norm derivative
413  const double rT1_norm_derivative = VectorNormDerivative(rT1_norm, rT1, rT1_derivative);
414 
415  const array_1d<double, 3>& unit_rT1_derivative =
416  UnitVectorDerivative(rT1_norm, rT1_norm_derivative, rT1, rT1_derivative);
417 
418  rOutput(1, 0) = unit_rT1_derivative[0];
419  rOutput(1, 1) = unit_rT1_derivative[1];
420  rOutput(1, 2) = unit_rT1_derivative[2];
421 
422  rOutput(2, 0) = unit_normal_derivative[1] * unit_rT1[2]
423  + unit_normal[1] * unit_rT1_derivative[2]
424  - unit_normal_derivative[2] * unit_rT1[1]
425  - unit_normal[2] * unit_rT1_derivative[1];
426 
427 
428  rOutput(2, 1) = unit_normal_derivative[2] * unit_rT1[0]
429  + unit_normal[2] * unit_rT1_derivative[0]
430  - unit_normal_derivative[0] * unit_rT1[2]
431  - unit_normal[0] * unit_rT1_derivative[2];
432 
433  rOutput(2, 2) = unit_normal_derivative[0] * unit_rT1[1]
434  + unit_normal[0] * unit_rT1_derivative[1]
435  - unit_normal_derivative[1] * unit_rT1[0]
436  - unit_normal[1] * unit_rT1_derivative[0];
437 
438  KRATOS_CATCH("");
439  }
440 
442 
447  virtual void Rotate(TLocalMatrixType& rLocalMatrix,
448  TLocalVectorType& rLocalVector,
449  GeometryType& rGeometry) const
450  {
451  if(mBlockSize != mDomainSize) //Monolithic case
452  {
453  if(mDomainSize == 2) RotateAux<2,3>(rLocalMatrix,rLocalVector,rGeometry);
454  if(mDomainSize == 3) RotateAux<3,4>(rLocalMatrix,rLocalVector,rGeometry);
455  }
456  else //fractional step case
457  {
458  if(mDomainSize == 2) RotateAuxPure<2>(rLocalMatrix,rLocalVector,rGeometry);
459  if(mDomainSize == 3) RotateAuxPure<3>(rLocalMatrix,rLocalVector,rGeometry);
460  }
461 
462  }
463 
465  virtual void Rotate(TLocalVectorType& rLocalVector,
466  GeometryType& rGeometry) const
467  {
468  //const unsigned int LocalSize = rLocalVector.size(); // We expect this to work both with elements (4 nodes) and conditions (3 nodes)
469 
470  //unsigned int Index = 0;
471 
472  if (rLocalVector.size() > 0)
473  {
474  if(mBlockSize != mDomainSize) //Monolithic case
475  {
476  for(unsigned int j = 0; j < rGeometry.PointsNumber(); ++j)
477  {
478  if( this->IsSlip(rGeometry[j]) )
479  {
480  if(mDomainSize == 3)
481  {
482  array_1d<double,4> aux,aux1;
483  BoundedMatrix<double,4,4> rRot;
484  LocalRotationOperator3D<4>(rRot,rGeometry[j]);
485 
486  for(unsigned int k=0; k<4; k++)
487  aux[k] = rLocalVector[j*mBlockSize+k];
488 
489  noalias(aux1) = prod(rRot,aux);
490 
491  for(unsigned int k=0; k<4; k++)
492  rLocalVector[j*mBlockSize+k] = aux1[k];
493  }
494  else
495  {
496  array_1d<double,3> aux,aux1;
497  BoundedMatrix<double,3,3> rRot;
498  LocalRotationOperator2D<3>(rRot,rGeometry[j]);
499 
500  for(unsigned int k=0; k<3; k++)
501  {
502  aux[k] = rLocalVector[j*mBlockSize+k];
503  }
504 
505  noalias(aux1) = prod(rRot,aux);
506 
507  for(unsigned int k=0; k<3; k++)
508  rLocalVector[j*mBlockSize+k] = aux1[k];
509  }
510  }
511  //Index += mBlockSize;
512  }
513 
514  }
515  else //fractional step case
516  {
517  for(unsigned int j = 0; j < rGeometry.PointsNumber(); ++j)
518  {
519  if( this->IsSlip(rGeometry[j]) )
520  {
521  if(mDomainSize == 3)
522  {
523  array_1d<double,3> aux,aux1;
524  BoundedMatrix<double,3,3> rRot;
525  LocalRotationOperatorPure(rRot,rGeometry[j]);
526 
527  for(unsigned int k=0; k<3; k++)
528  aux[k] = rLocalVector[j*mBlockSize+k];
529 
530  noalias(aux1) = prod(rRot,aux);
531 
532  for(unsigned int k=0; k<3; k++)
533  rLocalVector[j*mBlockSize+k] = aux1[k];
534  }
535  else
536  {
537  array_1d<double,2> aux,aux1;
538  BoundedMatrix<double,2,2> rRot;
539  LocalRotationOperatorPure(rRot,rGeometry[j]);
540 
541  for(unsigned int k=0; k<2; k++)
542  aux[k] = rLocalVector[j*mBlockSize+k];
543 
544  noalias(aux1) = prod(rRot,aux);
545 
546  for(unsigned int k=0; k<2; k++)
547  rLocalVector[j*mBlockSize+k] = aux1[k];
548  }
549  }
550  //Index += mBlockSize;
551  }
552 
553  }
554 
555  }
556 
557  }
558 
560 
565  virtual void ApplySlipCondition(TLocalMatrixType& rLocalMatrix,
566  TLocalVectorType& rLocalVector,
567  GeometryType& rGeometry) const
568  {
569  const unsigned int LocalSize = rLocalVector.size(); // We expect this to work both with elements (4 nodes) and conditions (3 nodes)
570 
571  if (LocalSize > 0)
572  {
573  for(unsigned int itNode = 0; itNode < rGeometry.PointsNumber(); ++itNode)
574  {
575  if( this->IsSlip(rGeometry[itNode]))
576  {
577  // We fix the first dof (normal velocity) for each rotated block
578  unsigned int j = itNode * mBlockSize;
579  //const double k = rLocalMatrix(j,j)+rLocalMatrix(j,j+1)+rLocalMatrix(j,j+2);
580 
581  // If the mesh is moving, we must impose v_normal = vmesh_normal
582  array_1d<double,3> VMesh = rGeometry[itNode].FastGetSolutionStepValue(MESH_VELOCITY);
583  VMesh -= rGeometry[itNode].FastGetSolutionStepValue(VELOCITY);
584  array_1d<double,3> rN = rGeometry[itNode].FastGetSolutionStepValue(NORMAL);
585  this->Normalize(rN);
586 
587  for( unsigned int i = 0; i < j; ++i)// Skip term (i,i)
588  {
589  rLocalMatrix(i,j) = 0.0;
590  rLocalMatrix(j,i) = 0.0;
591  }
592  for( unsigned int i = j+1; i < LocalSize; ++i)
593  {
594  rLocalMatrix(i,j) = 0.0;
595  rLocalMatrix(j,i) = 0.0;
596  }
597 
598  rLocalVector(j) = inner_prod(rN,VMesh);
599  rLocalMatrix(j,j) = 1.0;
600  }
601  }
602  }
603  }
604 
606  virtual void ApplySlipCondition(TLocalVectorType& rLocalVector,
607  GeometryType& rGeometry) const
608  {
609  if (rLocalVector.size() > 0)
610  {
611  for(unsigned int itNode = 0; itNode < rGeometry.PointsNumber(); ++itNode)
612  {
613  if( this->IsSlip(rGeometry[itNode]) )
614  {
615  // We fix the first dof (normal velocity) for each rotated block
616  unsigned int j = itNode * mBlockSize;
617 
618  // If the mesh is moving, we must impose v_normal = vmesh_normal
619  array_1d<double,3> VMesh = rGeometry[itNode].FastGetSolutionStepValue(MESH_VELOCITY);
620  VMesh -= rGeometry[itNode].FastGetSolutionStepValue(VELOCITY);
621  array_1d<double,3> rN = rGeometry[itNode].FastGetSolutionStepValue(NORMAL);
622  this->Normalize(rN);
623 
624  rLocalVector[j] = inner_prod(rN,VMesh);
625  }
626  }
627  }
628  }
629 
631  virtual void RotateVelocities(ModelPart& rModelPart) const
632  {
633  TLocalVectorType Vel(mDomainSize);
634  TLocalVectorType Tmp(mDomainSize);
635 
636  ModelPart::NodeIterator it_begin = rModelPart.NodesBegin();
637 #pragma omp parallel for firstprivate(Vel,Tmp)
638  for(int iii=0; iii<static_cast<int>(rModelPart.Nodes().size()); iii++)
639  {
640  ModelPart::NodeIterator itNode = it_begin+iii;
641  if( this->IsSlip(*itNode) )
642  {
643  //this->RotationOperator<TLocalMatrixType>(Rotation,);
644  if(mDomainSize == 3)
645  {
646  BoundedMatrix<double,3,3> rRot;
647  LocalRotationOperatorPure(rRot,*itNode);
648 
649  array_1d<double,3>& rVelocity = itNode->FastGetSolutionStepValue(VELOCITY);
650  for(unsigned int i = 0; i < 3; i++) Vel[i] = rVelocity[i];
651  noalias(Tmp) = prod(rRot,Vel);
652  for(unsigned int i = 0; i < 3; i++) rVelocity[i] = Tmp[i];
653  }
654  else
655  {
656  BoundedMatrix<double,2,2> rRot;
657  LocalRotationOperatorPure(rRot,*itNode);
658 
659  array_1d<double,3>& rVelocity = itNode->FastGetSolutionStepValue(VELOCITY);
660  for(unsigned int i = 0; i < 2; i++) Vel[i] = rVelocity[i];
661  noalias(Tmp) = prod(rRot,Vel);
662  for(unsigned int i = 0; i < 2; i++) rVelocity[i] = Tmp[i];
663  }
664  }
665  }
666  }
667 
669  virtual void RecoverVelocities(ModelPart& rModelPart) const
670  {
671  TLocalVectorType Vel(mDomainSize);
672  TLocalVectorType Tmp(mDomainSize);
673 
674  ModelPart::NodeIterator it_begin = rModelPart.NodesBegin();
675 #pragma omp parallel for firstprivate(Vel,Tmp)
676  for(int iii=0; iii<static_cast<int>(rModelPart.Nodes().size()); iii++)
677  {
678  ModelPart::NodeIterator itNode = it_begin+iii;
679  if( this->IsSlip(*itNode) )
680  {
681  if(mDomainSize == 3)
682  {
683  BoundedMatrix<double,3,3> rRot;
684  LocalRotationOperatorPure(rRot,*itNode);
685 
686  array_1d<double,3>& rVelocity = itNode->FastGetSolutionStepValue(VELOCITY);
687  for(unsigned int i = 0; i < 3; i++) Vel[i] = rVelocity[i];
688  noalias(Tmp) = prod(trans(rRot),Vel);
689  for(unsigned int i = 0; i < 3; i++) rVelocity[i] = Tmp[i];
690  }
691  else
692  {
693  BoundedMatrix<double,2,2> rRot;
694  LocalRotationOperatorPure(rRot,*itNode);
695 
696  array_1d<double,3>& rVelocity = itNode->FastGetSolutionStepValue(VELOCITY);
697  for(unsigned int i = 0; i < 2; i++) Vel[i] = rVelocity[i];
698  noalias(Tmp) = prod(trans(rRot),Vel);
699  for(unsigned int i = 0; i < 2; i++) rVelocity[i] = Tmp[i];
700  }
701  }
702  }
703  }
704 
708 
712 
716 
718  virtual std::string Info() const
719  {
720  std::stringstream buffer;
721  buffer << "CoordinateTransformationUtils";
722  return buffer.str();
723  }
724 
726  virtual void PrintInfo(std::ostream& rOStream) const
727  {
728  rOStream << "CoordinateTransformationUtils";
729  }
730 
732  virtual void PrintData(std::ostream& rOStream) const {}
733 
737 
739 
740 protected:
743 
747 
751 
755 
756  template<unsigned int TDim, unsigned int TBlockSize, unsigned int TSkip = 0>
757  void RotateAux(TLocalMatrixType& rLocalMatrix,
758  TLocalVectorType& rLocalVector,
759  GeometryType& rGeometry) const
760  {
761  const unsigned int LocalSize = rLocalVector.size();
762 
763  //unsigned int Index = 0;
764  int rotations_needed = 0;
765  const unsigned int NumBlocks = LocalSize / TBlockSize;
766  DenseVector<bool> NeedRotation( NumBlocks, false);
767 
768  std::vector< BoundedMatrix<double,TBlockSize,TBlockSize> > rRot(NumBlocks);
769  for(unsigned int j = 0; j < NumBlocks; ++j)
770  {
771  if( this->IsSlip(rGeometry[j]) )
772  {
773  NeedRotation[j] = true;
774  rotations_needed++;
775 
776  if constexpr (TDim == 2) LocalRotationOperator2D<TBlockSize,TSkip>(rRot[j],rGeometry[j]);
777  else LocalRotationOperator3D<TBlockSize,TSkip>(rRot[j],rGeometry[j]);
778  }
779 
780  //Index += TBlockSize;
781  }
782 
783  if(rotations_needed > 0)
784  {
785  BoundedMatrix<double,TBlockSize,TBlockSize> mat_block, tmp;
786  array_1d<double,TBlockSize> aux, aux1;
787 
788  for(unsigned int i=0; i<NumBlocks; i++)
789  {
790  if(NeedRotation[i] == true)
791  {
792  for(unsigned int j=0; j<NumBlocks; j++)
793  {
794  if(NeedRotation[j] == true)
795  {
796  ReadBlockMatrix<TBlockSize>(mat_block, rLocalMatrix, i*TBlockSize, j*TBlockSize);
797  noalias(tmp) = prod(mat_block,trans(rRot[j]));
798  noalias(mat_block) = prod(rRot[i],tmp);
799  WriteBlockMatrix<TBlockSize>(mat_block, rLocalMatrix, i*TBlockSize, j*TBlockSize);
800  }
801  else
802  {
803  ReadBlockMatrix<TBlockSize>(mat_block, rLocalMatrix, i*TBlockSize, j*TBlockSize);
804  noalias(tmp) = prod(rRot[i],mat_block);
805  WriteBlockMatrix<TBlockSize>(tmp, rLocalMatrix, i*TBlockSize, j*TBlockSize);
806  }
807  }
808 
809  for(unsigned int k=0; k<TBlockSize; k++)
810  aux[k] = rLocalVector[i*TBlockSize+k];
811 
812  noalias(aux1) = prod(rRot[i],aux);
813 
814  for(unsigned int k=0; k<TBlockSize; k++)
815  rLocalVector[i*TBlockSize+k] = aux1[k];
816 
817  }
818  else
819  {
820  for(unsigned int j=0; j<NumBlocks; j++)
821  {
822  if(NeedRotation[j] == true)
823  {
824  ReadBlockMatrix<TBlockSize>(mat_block, rLocalMatrix, i*TBlockSize, j*TBlockSize);
825  noalias(tmp) = prod(mat_block,trans(rRot[j]));
826  WriteBlockMatrix<TBlockSize>(tmp, rLocalMatrix, i*TBlockSize, j*TBlockSize);
827  }
828  }
829  }
830 
831  }
832  }
833  }
834 
835  //to be used when there is only velocity (no additional pressure or other var block)
836  template<unsigned int TDim>
837  void RotateAuxPure(TLocalMatrixType& rLocalMatrix,
838  TLocalVectorType& rLocalVector,
839  GeometryType& rGeometry) const
840  {
841  const unsigned int LocalSize = rLocalVector.size();
842 
843  //unsigned int Index = 0;
844  int rotations_needed = 0;
845  const unsigned int NumBlocks = LocalSize / mBlockSize;
846  DenseVector<bool> NeedRotation( NumBlocks, false);
847 
848  std::vector< BoundedMatrix<double,TDim,TDim> > rRot(NumBlocks);
849  for(unsigned int j = 0; j < NumBlocks; ++j)
850  {
851  if( this->IsSlip(rGeometry[j]) )
852  {
853  NeedRotation[j] = true;
854  rotations_needed++;
855 
856  LocalRotationOperatorPure(rRot[j],rGeometry[j]);
857  }
858 
859  //Index += mBlockSize;
860  }
861 
862  if(rotations_needed > 0)
863  {
864  BoundedMatrix<double,TDim,TDim> mat_block, tmp;
865  array_1d<double,TDim> aux, aux1;
866 
867  for(unsigned int i=0; i<NumBlocks; i++)
868  {
869  if(NeedRotation[i] == true)
870  {
871  for(unsigned int j=0; j<NumBlocks; j++)
872  {
873  if(NeedRotation[j] == true)
874  {
875  ReadBlockMatrix<TDim>(mat_block, rLocalMatrix, i*mBlockSize, j*mBlockSize);
876  noalias(tmp) = prod(mat_block,trans(rRot[j]));
877  noalias(mat_block) = prod(rRot[i],tmp);
878  WriteBlockMatrix<TDim>(mat_block, rLocalMatrix, i*mBlockSize, j*mBlockSize);
879  }
880  else
881  {
882  ReadBlockMatrix<TDim>(mat_block, rLocalMatrix, i*mBlockSize, j*mBlockSize);
883  noalias(tmp) = prod(rRot[i],mat_block);
884  WriteBlockMatrix<TDim>(tmp, rLocalMatrix, i*mBlockSize, j*mBlockSize);
885  }
886  }
887 
888  for(unsigned int k=0; k<TDim; k++)
889  aux[k] = rLocalVector[i*mBlockSize+k];
890 
891  noalias(aux1) = prod(rRot[i],aux);
892 
893  for(unsigned int k=0; k<TDim; k++)
894  rLocalVector[i*mBlockSize+k] = aux1[k];
895 
896  }
897  else
898  {
899  for(unsigned int j=0; j<NumBlocks; j++)
900  {
901  if(NeedRotation[j] == true)
902  {
903  ReadBlockMatrix<TDim>(mat_block, rLocalMatrix, i*mBlockSize, j*mBlockSize);
904  noalias(tmp) = prod(mat_block,trans(rRot[j]));
905  WriteBlockMatrix<TDim>(tmp, rLocalMatrix, i*mBlockSize, j*mBlockSize);
906  }
907  }
908  }
909 
910  }
911  }
912  }
913 
914  template<unsigned int TBlockSize, unsigned int TSkip = 0>
915  void LocalRotationOperator2D(
916  BoundedMatrix<double,TBlockSize,TBlockSize>& rRot,
917  GeometryType::PointType& rThisPoint) const
918  {
919  noalias(rRot) = IdentityMatrix(TBlockSize);
920 
921  // Get the normal evaluated at the node
922  const array_1d<double,3>& rNormal = rThisPoint.FastGetSolutionStepValue(NORMAL);
923 
924  double aux = rNormal[0]*rNormal[0] + rNormal[1]*rNormal[1];
925  aux = sqrt(aux);
926 
927  rRot(TSkip ,TSkip ) = rNormal[0]/aux;
928  rRot(TSkip ,TSkip+1) = rNormal[1]/aux;
929  rRot(TSkip+1,TSkip ) = -rNormal[1]/aux;
930  rRot(TSkip+1,TSkip+1) = rNormal[0]/aux;
931  }
932 
933  template<unsigned int TBlockSize, unsigned int TSkip = 0>
934  void LocalRotationOperator3D(
935  BoundedMatrix<double,TBlockSize,TBlockSize>& rRot,
936  GeometryType::PointType& rThisPoint) const
937  {
938  noalias(rRot) = IdentityMatrix(TBlockSize);
939 
940  // Get the normal evaluated at the node
941  const array_1d<double,3>& rNormal = rThisPoint.FastGetSolutionStepValue(NORMAL);
942 
943  double aux = rNormal[0]*rNormal[0] + rNormal[1]*rNormal[1] + rNormal[2]*rNormal[2];
944  aux = sqrt(aux);
945  rRot(TSkip,TSkip ) = rNormal[0]/aux;
946  rRot(TSkip,TSkip+1) = rNormal[1]/aux;
947  rRot(TSkip,TSkip+2) = rNormal[2]/aux;
948  // Define the new coordinate system, where the first vector is aligned with the normal
949 
950  // To choose the remaining two vectors, we project the first component of the cartesian base to the tangent plane
951  array_1d<double,3> rT1;
952  rT1(0) = 1.0;
953  rT1(1) = 0.0;
954  rT1(2) = 0.0;
955  double dot = rRot(TSkip,TSkip);//this->Dot(rN,rT1);
956 
957  // It is possible that the normal is aligned with (1,0,0), resulting in norm(rT1) = 0
958  // If this is the case, repeat the procedure using (0,1,0)
959  if ( fabs(dot) > 0.99 )
960  {
961  rT1(0) = 0.0;
962  rT1(1) = 1.0;
963  rT1(2) = 0.0;
964 
965  dot = rRot(TSkip,TSkip+1); //this->Dot(rN,rT1);
966  }
967 
968  // calculate projection and normalize
969  rT1[0] -= dot*rRot(TSkip,TSkip);
970  rT1[1] -= dot*rRot(TSkip,TSkip+1);
971  rT1[2] -= dot*rRot(TSkip,TSkip+2);
972  this->Normalize(rT1);
973  rRot(TSkip+1,TSkip ) = rT1[0];
974  rRot(TSkip+1,TSkip+1) = rT1[1];
975  rRot(TSkip+1,TSkip+2) = rT1[2];
976 
977  // The third base component is choosen as N x T1, which is normalized by construction
978  rRot(TSkip+2,TSkip ) = rRot(TSkip,TSkip+1)*rT1[2] - rRot(TSkip,TSkip+2)*rT1[1];
979  rRot(TSkip+2,TSkip+1) = rRot(TSkip,TSkip+2)*rT1[0] - rRot(TSkip,TSkip )*rT1[2];
980  rRot(TSkip+2,TSkip+2) = rRot(TSkip,TSkip )*rT1[1] - rRot(TSkip,TSkip+1)*rT1[0];
981  }
982 
983  bool IsSlip(const Node& rNode) const
984  {
985  return rNode.Is(mrFlag);
986  }
987 
989 
993  template< class TVectorType >
994  double Normalize(TVectorType& rThis) const
995  {
996  double Norm = 0;
997  for(typename TVectorType::iterator iComponent = rThis.begin(); iComponent < rThis.end(); ++iComponent)
998  Norm += (*iComponent)*(*iComponent);
999  Norm = sqrt(Norm);
1000  for(typename TVectorType::iterator iComponent = rThis.begin(); iComponent < rThis.end(); ++iComponent)
1001  *iComponent /= Norm;
1002  return Norm;
1003  }
1004 
1008 
1009  unsigned int GetDomainSize() const
1010  {
1011  return mDomainSize;
1012  }
1013 
1014  unsigned int GetBlockSize() const
1015  {
1016  return mBlockSize;
1017  }
1018 
1022 
1026 
1028 
1029 private:
1032 
1036 
1038  const unsigned int mDomainSize;
1039 
1041 
1043  const unsigned int mBlockSize;
1044 
1045  const Kratos::Flags& mrFlag;
1046 
1050 
1054 
1055 // /// Compute a rotation matrix to transform values from the cartesian base to one oriented with the node's normal
1056 // /**
1057 // * The normal is read from solution step data NORMAL. Use NormalCalculationUtils::CalculateOnSimplex to
1058 // * obtain and store the nodal normal from the normals of the model's conditons.
1059 // * @param rRot The rotation matrix (output)
1060 // * @param rThisPoint The point used to orient the new coordinate system.
1061 // * @see NormalCalculationUtils
1062 // */
1063 // template<class TMatrixType>
1064 // void RotationOperator(TMatrixType& rRot,
1065 // GeometryType::PointType& rThisPoint) const
1066 // {
1067 // typedef boost::numeric::ublas::matrix_row<TMatrixType> ThisRowType;
1068 // // Get the normal evaluated at the node
1069 // const array_1d<double,3>& rNormal = rThisPoint.FastGetSolutionStepValue(NORMAL);
1070 //
1071 // if(mDomainSize == 3)
1072 // {
1073 // // Define the new coordinate system, where the first vector is aligned with the normal
1074 // ThisRowType rN(rRot,0);
1075 // for( unsigned int i = 0; i < 3; ++i)
1076 // rN[i] = rNormal[i];
1077 // this->Normalize(rN);
1078 //
1079 // // To choose the remaining two vectors, we project the first component of the cartesian base to the tangent plane
1080 // ThisRowType rT1(rRot,1);
1081 // rT1(0) = 1.0;
1082 // rT1(1) = 0.0;
1083 // rT1(2) = 0.0;
1084 //
1085 // double dot = this->Dot(rN,rT1);
1086 //
1087 // // It is possible that the normal is aligned with (1,0,0), resulting in norm(rT1) = 0
1088 // // If this is the case, repeat the procedure using (0,1,0)
1089 // if ( fabs(dot) > 0.99 )
1090 // {
1091 // rT1(0) = 0.0;
1092 // rT1(1) = 1.0;
1093 // rT1(2) = 0.0;
1094 //
1095 // dot = this->Dot(rN,rT1);
1096 // }
1097 //
1098 // // calculate projection and normalize
1099 // rT1 -= dot * rN;
1100 // this->Normalize(rT1);
1101 //
1102 // // The third base component is choosen as N x T1, which is normalized by construction
1103 // ThisRowType rT2(rRot,2);
1104 // rT2(0) = rN(1)*rT1(2) - rN(2)*rT1(1);
1105 // rT2(1) = rN(2)*rT1(0) - rN(0)*rT1(2);
1106 // rT2(2) = rN(0)*rT1(1) - rN(1)*rT1(0);
1107 // }
1108 // else //if(mDomainSize == 2)
1109 // {
1110 // /* The basis for the new coordinate system is (normal,tangent)
1111 // Tangent vector is chosen (-normal_y, normal_x) so that the resulting base
1112 // is right-handed.
1113 // */
1114 // ThisRowType rN(rRot,0);
1115 // ThisRowType rT(rRot,1);
1116 //
1117 // rN[0] = rNormal[0];
1118 // rN[1] = rNormal[1];
1119 // this->Normalize(rN);
1120 // rT[0] = -rN[1];
1121 // rT[1] = rN[0];
1122 // }
1123 //
1124 // }
1125 
1126  template< class TVectorType >
1127  double Dot(const TVectorType& rV1,
1128  const TVectorType& rV2) const
1129  {
1130  double dot = 0.0;
1131  for( typename TVectorType::const_iterator iV1 = rV1.begin(),iV2 = rV2.begin(); iV1 != rV1.end(); ++iV1, ++iV2)
1132  {
1133  dot += (*iV1) * (*iV2);
1134  }
1135  return dot;
1136  }
1137 
1138  inline double VectorNormDerivative(
1139  const double ValueNorm,
1140  const array_1d<double, 3>& rValue,
1141  const array_1d<double, 3>& rValueDerivative) const
1142  {
1143  return inner_prod(rValue, rValueDerivative) / ValueNorm;
1144  }
1145 
1146  inline array_1d<double, 3> UnitVectorDerivative(
1147  const double VectorNorm,
1148  const double VectorNormDerivative,
1149  const array_1d<double, 3>& rVector,
1150  const array_1d<double, 3>& rVectorDerivative) const
1151  {
1152  return (rVectorDerivative * VectorNorm - rVector * VectorNormDerivative) /
1153  std::pow(VectorNorm, 2);
1154  }
1155 
1157 // void ApplyRotation(TLocalMatrixType& rMatrix,
1158 // const TLocalMatrixType& rRotation) const
1159 // {
1160 // // compute B = R*A*transpose(R)
1161 // const unsigned int LocalSize = rMatrix.size1();
1162 // const unsigned int NumBlocks = LocalSize / mBlockSize;
1163 // //TLocalMatrixType Tmp = ZeroMatrix(LocalSize,LocalSize);
1164 // /*
1165 // for (unsigned int iBlock = 0; iBlock < NumBlocks; iBlock++)
1166 // {
1167 // for (unsigned int jBlock = 0; jBlock < NumBlocks; jBlock++)
1168 // {
1169 // for (unsigned int i = iBlock*mBlockSize; i < (iBlock+1)*mBlockSize; i++)
1170 // {
1171 // for(unsigned int j = jBlock*mBlockSize; j < (jBlock+1)*mBlockSize; j++)
1172 // {
1173 // double& tij = Tmp(i,j);
1174 // for(unsigned int k = iBlock*mBlockSize; k < (iBlock+1)*mBlockSize; k++)
1175 // {
1176 // for(unsigned int l = jBlock*mBlockSize; l < (jBlock+1)*mBlockSize; l++)
1177 // {
1178 // tij += rRotation(i,k)*rMatrix(k,l)*rRotation(j,l);
1179 // }
1180 // }
1181 // }
1182 // }
1183 // }
1184 // }*/
1185 //
1186 // Matrix Tmp = prod(rMatrix,trans(rRotation));
1187 // noalias(rMatrix) = prod(rRotation,Tmp);
1188 //
1189 // // noalias(rMatrix) = Tmp;
1190 // }
1191 
1192  //auxiliary functions
1193  template< unsigned int TBlockSize >
1194  void ReadBlockMatrix( BoundedMatrix<double,TBlockSize, TBlockSize>& block, const Matrix& origin, const unsigned int Ibegin, const unsigned int Jbegin) const
1195  {
1196  for(unsigned int i=0; i<TBlockSize; i++)
1197  {
1198  for(unsigned int j=0; j<TBlockSize; j++)
1199  {
1200  block(i,j) = origin(Ibegin+i, Jbegin+j);
1201  }
1202  }
1203  }
1204 
1205  template< unsigned int TBlockSize >
1206  void WriteBlockMatrix( const BoundedMatrix<double,TBlockSize, TBlockSize>& block, Matrix& destination, const unsigned int Ibegin, const unsigned int Jbegin) const
1207  {
1208  for(unsigned int i=0; i<TBlockSize; i++)
1209  {
1210  for(unsigned int j=0; j<TBlockSize; j++)
1211  {
1212  destination(Ibegin+i, Jbegin+j) = block(i,j);
1213  }
1214  }
1215  }
1216 
1220 
1224 
1228 
1230  CoordinateTransformationUtils& operator=(CoordinateTransformationUtils const& rOther) {}
1231 
1233  CoordinateTransformationUtils(CoordinateTransformationUtils const& rOther) {}
1234 
1236 };
1237 
1239 
1242 
1246 
1248 template<class TLocalMatrixType, class TLocalVectorType, class TValueType>
1249 inline std::istream& operator >>(std::istream& rIStream,
1250  CoordinateTransformationUtils<TLocalMatrixType, TLocalVectorType,
1251  TValueType>& rThis) {
1252  return rIStream;
1253 }
1254 
1256 template<class TLocalMatrixType, class TLocalVectorType, class TValueType>
1257 inline std::ostream& operator <<(std::ostream& rOStream,
1258  const CoordinateTransformationUtils<TLocalMatrixType, TLocalVectorType,
1259  TValueType>& rThis) {
1260  rThis.PrintInfo(rOStream);
1261  rOStream << std::endl;
1262  rThis.PrintData(rOStream);
1263 
1264  return rOStream;
1265 }
1266 
1268 
1270 
1271 }
1272 
1273 #endif // KRATOS_COORDINATE_TRANSFORMATION_UTILITIES_H
Kratos::CoordinateTransformationUtils
Definition: coordinate_transformation_utilities.h:57
Kratos::CoordinateTransformationUtils::PrintData
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: coordinate_transformation_utilities.h:732
Kratos::CoordinateTransformationUtils::RotateVelocities
virtual void RotateVelocities(ModelPart &rModelPart) const
Transform nodal velocities to the rotated coordinates (aligned with each node's normal)
Definition: coordinate_transformation_utilities.h:631
Kratos::CoordinateTransformationUtils::GeometryType
Geometry< Node > GeometryType
Definition: coordinate_transformation_utilities.h:67
Kratos::CoordinateTransformationUtils::Normalize
double Normalize(TVectorType &rThis) const
Normalize a vector.
Definition: coordinate_transformation_utilities.h:994
Kratos::CoordinateTransformationUtils::Rotate
virtual void Rotate(TLocalVectorType &rLocalVector, GeometryType &rGeometry) const
RHS only version of Rotate.
Definition: coordinate_transformation_utilities.h:465
Kratos::CoordinateTransformationUtils::RotateAux
void RotateAux(TLocalMatrixType &rLocalMatrix, TLocalVectorType &rLocalVector, GeometryType &rGeometry) const
Definition: coordinate_transformation_utilities.h:757
Kratos::CoordinateTransformationUtils::CalculateRotationOperatorPureShapeSensitivities
virtual void CalculateRotationOperatorPureShapeSensitivities(BoundedMatrix< double, 3, 3 > &rOutput, const std::size_t DerivativeNodeIndex, const std::size_t DerivativeDirectionIndex, const GeometryType::PointType &rThisPoint) const
Calculate 3d rotation nodal matrix shape sensitivities.
Definition: coordinate_transformation_utilities.h:345
Kratos::CoordinateTransformationUtils::RotateAuxPure
void RotateAuxPure(TLocalMatrixType &rLocalMatrix, TLocalVectorType &rLocalVector, GeometryType &rGeometry) const
Definition: coordinate_transformation_utilities.h:837
Kratos::CoordinateTransformationUtils::Info
virtual std::string Info() const
Turn back information as a string.
Definition: coordinate_transformation_utilities.h:718
Kratos::CoordinateTransformationUtils::CalculateRotationOperatorPure
virtual void CalculateRotationOperatorPure(TLocalMatrixType &rRotationMatrix, const GeometryType::PointType &rThisPoint) const
Calculates rotation operator for given point.
Definition: coordinate_transformation_utilities.h:110
Kratos::CoordinateTransformationUtils::CoordinateTransformationUtils
CoordinateTransformationUtils(const unsigned int DomainSize, const unsigned int NumRowsPerNode, const Kratos::Flags &rSelectionFlag=SLIP)
Constructor.
Definition: coordinate_transformation_utilities.h:82
Kratos::CoordinateTransformationUtils::NodeType
Node NodeType
Definition: coordinate_transformation_utilities.h:65
Kratos::CoordinateTransformationUtils::~CoordinateTransformationUtils
virtual ~CoordinateTransformationUtils()
Destructor.
Definition: coordinate_transformation_utilities.h:91
Kratos::CoordinateTransformationUtils::RecoverVelocities
virtual void RecoverVelocities(ModelPart &rModelPart) const
Transform nodal velocities from the rotated system to the original one.
Definition: coordinate_transformation_utilities.h:669
Kratos::CoordinateTransformationUtils::LocalRotationOperator2D
void LocalRotationOperator2D(BoundedMatrix< double, TBlockSize, TBlockSize > &rRot, GeometryType::PointType &rThisPoint) const
Definition: coordinate_transformation_utilities.h:915
Kratos::CoordinateTransformationUtils::ApplySlipCondition
virtual void ApplySlipCondition(TLocalMatrixType &rLocalMatrix, TLocalVectorType &rLocalVector, GeometryType &rGeometry) const
Apply slip boundary conditions to the rotated local contributions.
Definition: coordinate_transformation_utilities.h:565
Kratos::CoordinateTransformationUtils::IsSlip
bool IsSlip(const Node &rNode) const
Definition: coordinate_transformation_utilities.h:983
Kratos::CoordinateTransformationUtils::ApplySlipCondition
virtual void ApplySlipCondition(TLocalVectorType &rLocalVector, GeometryType &rGeometry) const
RHS only version of ApplySlipCondition.
Definition: coordinate_transformation_utilities.h:606
Kratos::CoordinateTransformationUtils::GetBlockSize
unsigned int GetBlockSize() const
Definition: coordinate_transformation_utilities.h:1014
Kratos::CoordinateTransformationUtils::PrintInfo
virtual void PrintInfo(std::ostream &rOStream) const
Print information about this object.
Definition: coordinate_transformation_utilities.h:726
Kratos::CoordinateTransformationUtils::CalculateRotationOperatorPureShapeSensitivities
virtual void CalculateRotationOperatorPureShapeSensitivities(BoundedMatrix< double, 2, 2 > &rOutput, const std::size_t DerivativeNodeIndex, const std::size_t DerivativeDirectionIndex, const GeometryType::PointType &rThisPoint) const
Calculate 2d rotation nodal matrix shape sensitivities.
Definition: coordinate_transformation_utilities.h:267
Kratos::CoordinateTransformationUtils::LocalRotationOperatorPure
void LocalRotationOperatorPure(BoundedMatrix< double, 3, 3 > &rRot, const GeometryType::PointType &rThisPoint) const
Definition: coordinate_transformation_utilities.h:138
Kratos::CoordinateTransformationUtils::CalculateRotationOperatorPureShapeSensitivities
virtual void CalculateRotationOperatorPureShapeSensitivities(TLocalMatrixType &rRotationMatrixShapeDerivative, const std::size_t DerivativeNodeIndex, const std::size_t DerivativeDirectionIndex, const GeometryType::PointType &rThisPoint) const
Calculates rotation nodal matrix shape sensitivities.
Definition: coordinate_transformation_utilities.h:216
Kratos::CoordinateTransformationUtils::Rotate
virtual void Rotate(TLocalMatrixType &rLocalMatrix, TLocalVectorType &rLocalVector, GeometryType &rGeometry) const
Rotate the local system contributions so that they are oriented with each node's normal.
Definition: coordinate_transformation_utilities.h:447
Kratos::CoordinateTransformationUtils::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(CoordinateTransformationUtils)
Pointer definition of CoordinateTransformationUtils.
Kratos::CoordinateTransformationUtils::LocalRotationOperator3D
void LocalRotationOperator3D(BoundedMatrix< double, TBlockSize, TBlockSize > &rRot, GeometryType::PointType &rThisPoint) const
Definition: coordinate_transformation_utilities.h:934
Kratos::CoordinateTransformationUtils::LocalRotationOperatorPure
void LocalRotationOperatorPure(BoundedMatrix< double, 2, 2 > &rRot, const GeometryType::PointType &rThisPoint) const
Definition: coordinate_transformation_utilities.h:185
Kratos::CoordinateTransformationUtils::GetDomainSize
unsigned int GetDomainSize() const
Definition: coordinate_transformation_utilities.h:1009
Kratos::Flags
Definition: flags.h:58
Kratos::Flags::Is
bool Is(Flags const &rOther) const
Definition: flags.h:274
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::PointsNumber
SizeType PointsNumber() const
Definition: geometry.h:528
Kratos::Internals::Matrix
Definition: amatrix_interface.h:41
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::NodeIterator
MeshType::NodeIterator NodeIterator
Definition: model_part.h:134
Kratos::ModelPart::Nodes
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::Node::FastGetSolutionStepValue
TVariableType::Type & FastGetSolutionStepValue(const TVariableType &rThisVariable)
Definition: node.h:435
Kratos::Node::Id
IndexType Id() const
Definition: node.h:262
Kratos::Node::SolutionStepsDataHas
bool SolutionStepsDataHas(const VariableData &rThisVariable) const
Definition: node.h:427
Kratos::Node::GetValue
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: node.h:466
Kratos::Node::Has
bool Has(const Variable< TDataType > &rThisVariable) const
Definition: node.h:498
Kratos::Point::Coordinates
CoordinatesArrayType const & Coordinates() const
Definition: point.h:215
Kratos::array_1d< double, 3 >
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
KRATOS_ERROR
#define KRATOS_ERROR
Definition: exception.h:161
KRATOS_ERROR_IF
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
KRATOS_DEBUG_ERROR_IF
#define KRATOS_DEBUG_ERROR_IF(conditional)
Definition: exception.h:171
model_part.h
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::norm_2
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
Kratos::IdentityMatrix
AMatrix::IdentityMatrix< double > IdentityMatrix
Definition: amatrix_interface.h:564
Kratos::Matrix
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
Kratos::prod
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
Kratos::inner_prod
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
Kratos::operator>>
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::row
AMatrix::MatrixRow< const TExpressionType > row(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression, std::size_t RowIndex)
Definition: amatrix_interface.h:649
Kratos::trans
AMatrix::TransposeMatrix< const T > trans(const T &TheMatrix)
Definition: amatrix_interface.h:486
Kratos::operator<<
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
bombardelli_test.coeff
list coeff
Definition: bombardelli_test.py:41
generate_total_lagrangian_mixed_volumetric_strain_element.tmp
tuple tmp
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:98
isotropic_damage_automatic_differentiation.aux1
float aux1
Definition: isotropic_damage_automatic_differentiation.py:239
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
sand_production_post_process_tool.Norm
def Norm(vector)
Definition: sand_production_post_process_tool.py:10
tensors::i
integer i
Definition: TensorModule.f:17
Generated on Sat Feb 3 2024 02:48:02 for KratosMultiphysics by doxygen 1.9.1