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.
fs_generalized_wall_condition.h
Go to the documentation of this file.
1 /*
2  ==============================================================================
3  Kratos Fluid Dynamics Application
4  Kratos
5  A General Purpose Software for Multi-Physics Finite Element Analysis
6  Version 1.0 (Released on march 05, 2007).
7 
8  Copyright 2007
9  Pooyan Dadvand, Riccardo Rossi
10  pooyan@cimne.upc.edu
11  rrossi@cimne.upc.edu
12  CIMNE (International Center for Numerical Methods in Engineering),
13  Gran Capita' s/n, 08034 Barcelona, Spain
14 
15  Permission is hereby granted, free of charge, to any person obtaining
16  a copy of this software and associated documentation files (the
17  "Software"), to deal in the Software without restriction, including
18  without limitation the rights to use, copy, modify, merge, publish,
19  distribute, sublicense and/or sell copies of the Software, and to
20  permit persons to whom the Software is furnished to do so, subject to
21  the following condition:
22 
23  Distribution of this code for any commercial purpose is permissible
24  ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNER.
25 
26  The above copyright notice and this permission notice shall be
27  included in all copies or substantial portions of the Software.
28 
29  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
30  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
31  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
32  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
33  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
34  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
35  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
36 
37  ==============================================================================
38  */
39 
40 #ifndef KRATOS_FS_GENERALIZED_WALL_CONDITION_H
41 #define KRATOS_FS_GENERALIZED_WALL_CONDITION_H
42 
43 // System includes
44 #include <iostream>
45 #include <algorithm>
46 
47 // External includes
48 
49 // Project includes
50 #include "includes/define.h"
51 #include "includes/serializer.h"
52 #include "includes/condition.h"
53 #include "includes/process_info.h"
54 #include "includes/kratos_flags.h"
56 #include "includes/cfd_variables.h"
57 
58 // Application includes
60 
61 namespace Kratos {
64 
67 
71 
75 
79 
83 
85 
98 template<unsigned int TDim, unsigned int TNumNodes = TDim>
100 {
101 public:
104 
107 
108  typedef Node NodeType;
109 
111 
113 
115 
117 
119 
120  typedef Element::WeakPointer ElementWeakPointerType;
121 
122  typedef Element::Pointer ElementPointerType;
123 
124  typedef std::size_t IndexType;
125 
126  typedef std::size_t SizeType;
127 
128  typedef std::vector<std::size_t> EquationIdVectorType;
129 
130  typedef std::vector< Dof<double>::Pointer > DofsVectorType;
131 
133 
136 
139 
143 
145 
148  FSGeneralizedWallCondition(IndexType NewId = 0) : Condition(NewId), mInitializeWasPerformed(false), mpElement()
149  {
150  }
151 
153 
158  : Condition(NewId, ThisNodes), mInitializeWasPerformed(false), mpElement()
159  {
160  }
161 
163 
167  FSGeneralizedWallCondition(IndexType NewId, GeometryType::Pointer pGeometry)
168  : Condition(NewId, pGeometry), mInitializeWasPerformed(false), mpElement()
169  {
170  }
171 
173 
178  FSGeneralizedWallCondition(IndexType NewId, GeometryType::Pointer pGeometry,
179  PropertiesType::Pointer pProperties)
180  : Condition(NewId, pGeometry, pProperties), mInitializeWasPerformed(false), mpElement()
181  {
182  }
183 
186  mInitializeWasPerformed(rOther.mInitializeWasPerformed), mpElement(rOther.mpElement)
187  {
188  }
189 
192  {}
193 
197 
200  {
201  Condition::operator=(rOther);
202  mInitializeWasPerformed = rOther.mInitializeWasPerformed;
203  mpElement = rOther.mpElement;
204  return *this;
205  }
206 
210 
212 
217  Condition::Pointer Create(
218  IndexType NewId,
219  NodesArrayType const& ThisNodes,
220  PropertiesType::Pointer pProperties) const override
221  {
222  return Kratos::make_intrusive<FSGeneralizedWallCondition>(NewId,GetGeometry().Create(ThisNodes), pProperties);
223  }
224 
226 
231  Condition::Pointer Create(
232  IndexType NewId,
233  GeometryType::Pointer pGeom,
234  PropertiesType::Pointer pProperties) const override
235  {
236  return Kratos::make_intrusive<FSGeneralizedWallCondition>(NewId, pGeom, pProperties);
237  }
238 
240  void Initialize(const ProcessInfo &rCurrentProcessInfo) override
241  {
242  KRATOS_TRY;
243 
244  if (this->Is(SLIP))
245  {
246  const array_1d<double, 3> &rNormal = this->GetValue(NORMAL);
247  KRATOS_ERROR_IF(norm_2(rNormal) == 0.0) << "NORMAL must be calculated before using this " << this->Info() << "\n";
248  }
249 
250  if (mInitializeWasPerformed)
251  {
252  return;
253  }
254 
255  mInitializeWasPerformed = true;
256 
257  KRATOS_ERROR_IF(this->GetValue(NEIGHBOUR_ELEMENTS).size() == 0) << this->Info() << " cannot find parent element\n";
258 
259  double EdgeLength;
260  array_1d<double, 3> Edge;
261  mpElement = this->GetValue(NEIGHBOUR_ELEMENTS)(0);
262  GeometryType &rElemGeom = mpElement->GetGeometry();
263 
264  Edge = rElemGeom[1].Coordinates() - rElemGeom[0].Coordinates();
265  mMinEdgeLength = Edge[0] * Edge[0];
266  for (SizeType d = 1; d < TDim; d++)
267  {
268  mMinEdgeLength += Edge[d] * Edge[d];
269  }
270 
271  for (SizeType j = 2; j < rElemGeom.PointsNumber(); j++)
272  {
273  for (SizeType k = 0; k < j; k++)
274  {
275  Edge = rElemGeom[j].Coordinates() - rElemGeom[k].Coordinates();
276  EdgeLength = Edge[0] * Edge[0];
277 
278  for (SizeType d = 1; d < TDim; d++)
279  {
280  EdgeLength += Edge[d] * Edge[d];
281  }
282 
283  mMinEdgeLength = (EdgeLength < mMinEdgeLength) ? EdgeLength : mMinEdgeLength;
284  }
285  }
286  mMinEdgeLength = sqrt(mMinEdgeLength);
287  return;
288 
289  KRATOS_CATCH("");
290  }
291 
293  MatrixType& rLeftHandSideMatrix,
294  const ProcessInfo& rCurrentProcessInfo) override
295  {
296  VectorType RHS;
297  this->CalculateLocalSystem(rLeftHandSideMatrix, RHS, rCurrentProcessInfo);
298  }
299 
301 
306  void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
307  VectorType& rRightHandSideVector,
308  const ProcessInfo& rCurrentProcessInfo) override
309  {
310  KRATOS_TRY;
311 
312  if (mInitializeWasPerformed == false)
313  {
314  Initialize(rCurrentProcessInfo);
315  }
316 
317  if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1)
318  {
319  // Initialize local contributions
320  const SizeType LocalSize = TDim * TNumNodes;
321 
322  if (rLeftHandSideMatrix.size1() != LocalSize)
323  {
324  rLeftHandSideMatrix.resize(LocalSize, LocalSize);
325  }
326  if (rRightHandSideVector.size() != LocalSize)
327  {
328  rRightHandSideVector.resize(LocalSize);
329  }
330 
331  noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
332  noalias(rRightHandSideVector) = ZeroVector(LocalSize);
333 
334  if (this->Is(SLIP))
335  this->ApplyWallLaw(rLeftHandSideMatrix, rRightHandSideVector);
336  }
337  else if (rCurrentProcessInfo[FRACTIONAL_STEP] == 5)
338  {
339  // add IAC penalty to local pressure system
340  const SizeType LocalSize = TNumNodes;
341 
342  if (rLeftHandSideMatrix.size1() != LocalSize)
343  {
344  rLeftHandSideMatrix.resize(LocalSize, LocalSize);
345  }
346  if (rRightHandSideVector.size() != LocalSize)
347  {
348  rRightHandSideVector.resize(LocalSize);
349  }
350 
351  noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
352  noalias(rRightHandSideVector) = ZeroVector(LocalSize);
353 
354  if (this->Is(INTERFACE))
355  this->ApplyIACPenalty(rLeftHandSideMatrix, rRightHandSideVector, rCurrentProcessInfo);
356  }
357  else
358  {
359  if (rLeftHandSideMatrix.size1() != 0)
360  {
361  rLeftHandSideMatrix.resize(0,0,false);
362  }
363 
364  if (rRightHandSideVector.size() != 0)
365  {
366  rRightHandSideVector.resize(0,false);
367  }
368  }
369 
370  KRATOS_CATCH("");
371  }
372 
374  int Check(const ProcessInfo& rCurrentProcessInfo) const override
375  {
376  KRATOS_TRY;
377 
378  int Check = Condition::Check(rCurrentProcessInfo); // Checks id > 0 and area >= 0
379 
380  if (Check != 0)
381  {
382  return Check;
383  }
384  else
385  {
386  // Check that the element's nodes contain all required SolutionStepData and Degrees of freedom
387  for(unsigned int i=0; i<this->GetGeometry().size(); ++i)
388  {
389  if(this->GetGeometry()[i].SolutionStepsDataHas(VELOCITY) == false)
390  KRATOS_THROW_ERROR(std::invalid_argument,"missing VELOCITY variable on solution step data for node ",this->GetGeometry()[i].Id());
391  if(this->GetGeometry()[i].SolutionStepsDataHas(MESH_VELOCITY) == false)
392  KRATOS_THROW_ERROR(std::invalid_argument,"missing MESH_VELOCITY variable on solution step data for node ",this->GetGeometry()[i].Id());
393  if(this->GetGeometry()[i].SolutionStepsDataHas(DENSITY) == false)
394  KRATOS_THROW_ERROR(std::invalid_argument,"missing DENSITY variable on solution step data for node ",this->GetGeometry()[i].Id());
395  if(this->GetGeometry()[i].SolutionStepsDataHas(VISCOSITY) == false)
396  KRATOS_THROW_ERROR(std::invalid_argument,"missing VISCOSITY variable on solution step data for node ",this->GetGeometry()[i].Id());
397  if(this->GetGeometry()[i].SolutionStepsDataHas(NORMAL) == false)
398  KRATOS_THROW_ERROR(std::invalid_argument,"missing NORMAL variable on solution step data for node ",this->GetGeometry()[i].Id());
399  if(this->GetGeometry()[i].HasDofFor(VELOCITY_X) == false ||
400  this->GetGeometry()[i].HasDofFor(VELOCITY_Y) == false ||
401  this->GetGeometry()[i].HasDofFor(VELOCITY_Z) == false)
402  KRATOS_THROW_ERROR(std::invalid_argument,"missing VELOCITY component degree of freedom on node ",this->GetGeometry()[i].Id());
403  if(this->GetGeometry()[i].HasDofFor(PRESSURE) == false)
404  KRATOS_THROW_ERROR(std::invalid_argument,"missing PRESSURE degree of freedom on node ",this->GetGeometry()[i].Id());
405  }
406 
407  return Check;
408  }
409 
410  KRATOS_CATCH("");
411  }
412 
414 
418  void EquationIdVector(EquationIdVectorType& rResult,
419  const ProcessInfo& rCurrentProcessInfo) const override;
420 
422 
426  void GetDofList(DofsVectorType& rConditionDofList,
427  const ProcessInfo& rCurrentProcessInfo) const override;
428 
430 
434  void GetValuesVector(Vector& Values, int Step = 0) const override
435  {
436  const SizeType LocalSize = TDim * TNumNodes;
437  unsigned int LocalIndex = 0;
438 
439  if (Values.size() != LocalSize)
440  {
441  Values.resize(LocalSize, false);
442  }
443 
444  for (unsigned int i = 0; i < TNumNodes; ++i)
445  {
446  const array_1d<double,3>& rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY, Step);
447  for (unsigned int d = 0; d < TDim; ++d)
448  {
449  Values[LocalIndex++] = rVelocity[d];
450  }
451  }
452  }
453 
457 
461 
465 
467  std::string Info() const override
468  {
469  std::stringstream buffer;
470  buffer << "FSGeneralizedWallCondition" << TDim << "D";
471  return buffer.str();
472  }
473 
475  void PrintInfo(std::ostream& rOStream) const override
476  { rOStream << "FSGeneralizedWallCondition";}
477 
479  void PrintData(std::ostream& rOStream) const override
480  {}
481 
485 
487 
488 protected:
491 
495 
499 
503 
505  {
506  return mpElement->shared_from_this();
507  }
508 
509  template< class TVariableType >
510  void EvaluateInPoint(TVariableType& rResult,
512  const ShapeFunctionsType& rShapeFunc)
513  {
514  GeometryType& rGeom = this->GetGeometry();
515 
516  rResult = rShapeFunc[0] * rGeom[0].FastGetSolutionStepValue(Var);
517 
518  for(SizeType i = 1; i < TNumNodes; i++)
519  {
520  rResult += rShapeFunc[i] * rGeom[i].FastGetSolutionStepValue(Var);
521  }
522  }
523 
526  {
527  GeometryType& rElemGeom = pGetElement()->GetGeometry();
528  const SizeType NumNodes = rElemGeom.PointsNumber();
530  Vector DetJ;
531  rElemGeom.ShapeFunctionsIntegrationPointsGradients(DN_DX, DetJ,
533  ShapeFunctionDerivativesType& rDN_DX = DN_DX[0];
534 
535  const double& pres = rElemGeom[0].FastGetSolutionStepValue(PRESSURE,1);
536  for (SizeType d = 0; d < TDim; ++d)
537  {
538  rResult[d] = rDN_DX(0,d) * pres;
539  }
540 
541  for (SizeType i = 1; i < NumNodes; i++)
542  {
543  const double& pres = rElemGeom[i].FastGetSolutionStepValue(PRESSURE,1);
544  for (SizeType d = 0; d < TDim; ++d)
545  {
546  rResult[d] += rDN_DX(i,d) * pres;
547  }
548  }
549  }
550 
552 
558  void CalculateWallParameters(double& rWallHeight, array_1d<double,3>& rWallVel, double& rWallGradP, double& rArea);
559 
561 
567  double EvaluateWallFunctionResidual(const double& rWallHeight, const double& rWallVel, const double& rWallStress, const double& rWallGradP)
568  {
569  const ShapeFunctionsType& N = row(this->GetGeometry().ShapeFunctionsValues(GeometryData::IntegrationMethod::GI_GAUSS_1),0);
570  double rho, nu, func1, func2, Vel1, Vel2, Vel12, YPlus1, YPlus2, sign1, sign2;
571  EvaluateInPoint(rho, DENSITY, N);
572  EvaluateInPoint(nu, VISCOSITY, N);
573  Vel1 = sqrt(fabs(rWallStress) / rho);
574  Vel2 = pow(nu * fabs(rWallGradP) / rho, 0.333333);
575  Vel12 = Vel1 + Vel2;
576 
577  if (Vel12 == 0.0)
578  {
579  Vel12 = 1.0;
580  }
581 
582  YPlus1 = Vel1 * rWallHeight / nu;
583  YPlus2 = Vel2 * rWallHeight / nu;
584 
585  // evaluate func1(YPlus1)
586  func1 = 0.0;
587  if (YPlus1 <= 5)
588  {
589  func1 = (1.0 + (0.01 - 0.0029 * YPlus1) * YPlus1) * YPlus1;
590  }
591  else if (YPlus1 <= 30)
592  {
593  func1 = -0.872 + (1.465 + (-0.0702 + (0.00166 - 0.00001495 * YPlus1) * YPlus1) * YPlus1) * YPlus1;
594  }
595  else if (YPlus1 <= 140)
596  {
597  func1 = 8.6 + (0.1864 + (-0.002006 + (0.00001144 - 0.00000002551 * YPlus1) * YPlus1) * YPlus1) * YPlus1;
598  }
599  else
600  {
601  func1 = 2.439 * log(YPlus1) + 5.0;
602  }
603 
604  // evaluate func2(YPlus2)
605  func2 = 0.0;
606  if (YPlus2 <= 4)
607  {
608  func2 = (0.5 - 0.00731 * YPlus2) * YPlus2 * YPlus2;
609  }
610  else if (YPlus2 <= 15)
611  {
612  func2 = -15.138 + (8.4688 + (-0.81976 + (0.037292 - 0.00063866 * YPlus2) * YPlus2) * YPlus2) * YPlus2;
613  }
614  else if (YPlus2 <= 30)
615  {
616  func2 = 11.925 + (0.934 + (-0.027805 + (0.00046262 - 0.0000031442 * YPlus2) * YPlus2) * YPlus2) * YPlus2;
617  }
618  else
619  {
620  func2 = 5.0 * log(YPlus2) + 8.0;
621  }
622 
623  sign1 = (rWallStress >= 0.0) ? 1.0 : -1.0;
624  sign2 = (rWallGradP >= 0.0) ? 1.0 : -1.0;
625  return (rWallVel - sign1 * Vel1 * func1 - sign2 * Vel2 * func2) / Vel12;
626  }
627 
629 
634  double CalculateWallStress(const double& rWallHeight, const double& rWallVel, const double& rWallGradP)
635  {
636  KRATOS_TRY;
637 
638  const SizeType MaxIter = 50;
639  const double Eps = 1.0e-08;
640  double A =-0.1;
641  double B = 0.1;
642  double C, Xm, P, Q, R, S, Min1, Min2, ResA, ResB, ResC, Tol;
643  double D=0.0;
644  double E=0.0;
645 
646  for (SizeType k=0; k < 10; ++k)
647  {
648  A *= 10.0;
649  B *= 10.0;
650  ResA = EvaluateWallFunctionResidual(rWallHeight,rWallVel,A,rWallGradP);
651  ResB = EvaluateWallFunctionResidual(rWallHeight,rWallVel,B,rWallGradP);
652 
653  if (ResB * ResA <= 0.0)
654  {
655  break;
656  }
657  }
658 
659  if (ResB * ResA > 0.0)
660  {
661  return 0.;
662  //KRATOS_THROW_ERROR(std::logic_error, "Maximum wall stress search width exceeded","");
663  }
664 
665  // root finding with Brent's algorithm
666  C = B;
667  ResC = ResB;
668 
669  for (SizeType i=1; i <= MaxIter; ++i)
670  {
671  if (ResB * ResC > 0.0)
672  {
673  C = A;
674  ResC = ResA;
675  D = B - A;
676  E = D;
677  }
678 
679  if (fabs(ResC) < fabs(ResB))
680  {
681  A = B;
682  B = C;
683  C = A;
684  ResA = ResB;
685  ResB = ResC;
686  ResC = ResA;
687  }
688 
689  Tol = 2.0 * Eps * fabs(B) + 0.5e-10;
690  Xm = 0.5 * (C - B);
691 
692  if (fabs(Xm) <= Tol || ResB == 0.0)
693  {
694  return B;
695  }
696 
697  if (fabs(E) >= Tol && fabs(ResA) > fabs(ResB))
698  {
699  S = ResB / ResA;
700 
701  if (A == C)
702  {
703  P = 2.0 * Xm * S;
704  Q = 1.0 - S;
705  }
706  else
707  {
708  Q = ResA / ResC;
709  R = ResB / ResC;
710  P = S * (2.0 * Xm * Q * (Q - R) - (B - A) * (R - 1.0));
711  Q = (Q - 1.0) * (R - 1.0) * (S - 1.0);
712  }
713 
714  if (P > 0.0)
715  {
716  Q = -Q;
717  }
718 
719  P = fabs(P);
720  Min1 = 3.0 * Xm * Q - fabs(Tol * Q);
721  Min2 = fabs(E * Q);
722 
723  if (2.0 * P < std::min(Min1, Min2))
724  {
725  E = D;
726  D = P / Q;
727  }
728  else
729  {
730  D = Xm;
731  E = D;
732  }
733  }
734  else
735  {
736  D = Xm;
737  E = D;
738  }
739 
740  A = B;
741  ResA = ResB;
742 
743  if (fabs(D) > Tol)
744  {
745  B += D;
746  }
747  else
748  {
749  if (Xm >= 0.0)
750  {
751  B += fabs(Tol);
752  }
753  else
754  {
755  B -= fabs(Tol);
756  }
757  }
758 
759  ResB = EvaluateWallFunctionResidual(rWallHeight,rWallVel,B,rWallGradP);
760  }
761 
762  return 0.0;
763 
764  KRATOS_CATCH("");
765  }
766 
768 
771  bool IsCorner(double AngleCosine)
772  {
773  GeometryType& rGeometry = this->GetGeometry();
774  const array_1d<double, 3>& rNormal = this->GetValue(NORMAL);
775  const double Area = norm_2(rNormal);
776 
777  for (SizeType iNode = 0; iNode < rGeometry.PointsNumber(); ++iNode)
778  {
779  const array_1d<double, 3>& rNodeNormal = rGeometry[iNode].FastGetSolutionStepValue(NORMAL);
780  if (inner_prod(rNormal, rNodeNormal) < AngleCosine * Area * norm_2(rNodeNormal))
781  return true;
782  }
783  return false;
784  }
785 
787 
791  void ApplyWallLaw(MatrixType& rLocalMatrix, VectorType& rLocalVector)
792  {
793  const unsigned int BlockSize = TDim;
794  double WallHeight, Area, WallVelMag, tmp, WallStress, WallGradP, WallForce;
795  array_1d<double,3> WallVel;
796  GeometryType& rGeometry = this->GetGeometry();
797 
798  CalculateWallParameters(WallHeight, WallVel, WallGradP, Area);
799  WallVelMag = norm_2(WallVel);
800 
801  if (IsCorner(0.966) == false)
802  {
803  WallStress = CalculateWallStress(WallHeight, WallVelMag, WallGradP);
804  WallForce = (Area / static_cast<double>(TDim)) * WallStress;
805 
806  for(SizeType i=0; i < rGeometry.PointsNumber(); ++i)
807  {
808  const NodeType& rNode = rGeometry[i];
809  if(rNode.GetValue(Y_WALL) != 0.0 && rNode.Is(SLIP))
810  {
811  WallVel = rNode.FastGetSolutionStepValue(VELOCITY,1) - rNode.FastGetSolutionStepValue(MESH_VELOCITY,1);
812  tmp = norm_2(WallVel);
813  WallVel /= (tmp != 0.0) ? tmp : 1.0;
814 
815  for (unsigned int d=0; d < TDim; d++)
816  {
817  unsigned int k = i * BlockSize + d;
818  rLocalVector[k] -= WallVel[d] * WallForce;
819  }
820  }
821  }
822  }
823  }
824 
826 
831  void ApplyIACPenalty(MatrixType& rLeftHandSideMatrix,
832  VectorType& rRightHandSideVector,
833  const ProcessInfo& rCurrentProcessInfo)
834  {
835  GeometryType& rGeometry = this->GetGeometry();
836  const array_1d<double,3>& rNormal = this->GetValue(NORMAL);
837  const double Area = norm_2(rNormal);
838  const double DiagonalTerm = Area / static_cast<double>(TNumNodes)
839  / (rCurrentProcessInfo[DENSITY] * rCurrentProcessInfo[BDF_COEFFICIENTS][0]);
840 
841  for(SizeType iNode=0; iNode < rGeometry.PointsNumber(); ++iNode)
842  {
843  rLeftHandSideMatrix(iNode,iNode) += DiagonalTerm;
844  }
845  }
846 
850 
854 
858 
860 
861 private:
864 
868  bool mInitializeWasPerformed;
869  double mMinEdgeLength;
870  ElementWeakPointerType mpElement;
871 
875 
876  friend class Serializer;
877 
878  void save(Serializer& rSerializer) const override
879  {
881  }
882 
883  void load(Serializer& rSerializer) override
884  {
886  }
887 
891 
895 
899 
903 
907 
909 
910 }; // Class FSGeneralizedWallCondition
911 
913 
916 
920 
922 template<unsigned int TDim, unsigned int TNumNodes>
923 inline std::istream& operator >>(std::istream& rIStream,
925 {
926  return rIStream;
927 }
928 
930 template<unsigned int TDim, unsigned int TNumNodes>
931 inline std::ostream& operator <<(std::ostream& rOStream,
933 {
934  rThis.PrintInfo(rOStream);
935  rOStream << std::endl;
936  rThis.PrintData(rOStream);
937 
938  return rOStream;
939 }
940 
942 
944 
945 }// namespace Kratos.
946 
947 #endif // KRATOS_FS_GENERALIZED_WALL_CONDITION_H
Base class for all Conditions.
Definition: condition.h:59
std::size_t SizeType
Definition: condition.h:94
std::vector< DofType::Pointer > DofsVectorType
Definition: condition.h:100
virtual int Check(const ProcessInfo &rCurrentProcessInfo) const
Definition: condition.h:854
Condition & operator=(Condition const &rOther)
Assignment operator.
Definition: condition.h:181
Implements a generalized wall model accounting for pressure gradients.
Definition: fs_generalized_wall_condition.h:100
Element::Pointer ElementPointerType
Definition: fs_generalized_wall_condition.h:122
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: fs_generalized_wall_condition.h:475
FSGeneralizedWallCondition(IndexType NewId, const NodesArrayType &ThisNodes)
Constructor using an array of nodes.
Definition: fs_generalized_wall_condition.h:157
FSGeneralizedWallCondition(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
Constructor using Properties.
Definition: fs_generalized_wall_condition.h:178
std::vector< Dof< double >::Pointer > DofsVectorType
Definition: fs_generalized_wall_condition.h:130
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: fs_generalized_wall_condition.h:479
void GetValuesVector(Vector &Values, int Step=0) const override
Returns VELOCITY_X, VELOCITY_Y, (VELOCITY_Z) for each node.
Definition: fs_generalized_wall_condition.h:434
std::string Info() const override
Turn back information as a string.
Definition: fs_generalized_wall_condition.h:467
Node NodeType
Definition: fs_generalized_wall_condition.h:108
Element::WeakPointer ElementWeakPointerType
Definition: fs_generalized_wall_condition.h:120
Geometry< NodeType >::PointsArrayType NodesArrayType
Definition: fs_generalized_wall_condition.h:116
void ApplyWallLaw(MatrixType &rLocalMatrix, VectorType &rLocalVector)
Compute the wall stress and add corresponding terms to the system contributions.
Definition: fs_generalized_wall_condition.h:791
FSGeneralizedWallCondition(FSGeneralizedWallCondition const &rOther)
Copy constructor.
Definition: fs_generalized_wall_condition.h:185
int Check(const ProcessInfo &rCurrentProcessInfo) const override
Check that all data required by this condition is available and reasonable.
Definition: fs_generalized_wall_condition.h:374
std::size_t SizeType
Definition: fs_generalized_wall_condition.h:126
std::vector< std::size_t > EquationIdVectorType
Definition: fs_generalized_wall_condition.h:128
~FSGeneralizedWallCondition() override
Destructor.
Definition: fs_generalized_wall_condition.h:191
void CalculateWallParameters(double &rWallHeight, array_1d< double, 3 > &rWallVel, double &rWallGradP, double &rArea)
Calculate input parameters to wall model.
void EvaluateOldPressureGradientInElement(array_1d< double, TDim > &rResult)
Calculate the pressure gradient using nodal data from the last time step.
Definition: fs_generalized_wall_condition.h:525
void EvaluateInPoint(TVariableType &rResult, const Kratos::Variable< TVariableType > &Var, const ShapeFunctionsType &rShapeFunc)
Definition: fs_generalized_wall_condition.h:510
void GetDofList(DofsVectorType &rConditionDofList, const ProcessInfo &rCurrentProcessInfo) const override
Returns a list of the condition's Dofs.
bool IsCorner(double AngleCosine)
Check angles between condition normal and node normals.
Definition: fs_generalized_wall_condition.h:771
FSGeneralizedWallCondition(IndexType NewId=0)
Default constructor.
Definition: fs_generalized_wall_condition.h:148
Vector ShapeFunctionsType
Definition: fs_generalized_wall_condition.h:132
Condition::Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const override
Create a new FSGeneralizedWallCondition object.
Definition: fs_generalized_wall_condition.h:217
Geometry< NodeType >::PointType PointType
Definition: fs_generalized_wall_condition.h:114
double CalculateWallStress(const double &rWallHeight, const double &rWallVel, const double &rWallGradP)
Calculate the wall stress from the wall function.
Definition: fs_generalized_wall_condition.h:634
GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType
Type for an array of shape function gradient matrices.
Definition: fs_generalized_wall_condition.h:138
std::size_t IndexType
Definition: fs_generalized_wall_condition.h:124
void Initialize(const ProcessInfo &rCurrentProcessInfo) override
Find the condition's parent element.
Definition: fs_generalized_wall_condition.h:240
Geometry< NodeType > GeometryType
Definition: fs_generalized_wall_condition.h:112
double EvaluateWallFunctionResidual(const double &rWallHeight, const double &rWallVel, const double &rWallStress, const double &rWallGradP)
Evaluate the residual of the generalized wall function.
Definition: fs_generalized_wall_condition.h:567
FSGeneralizedWallCondition(IndexType NewId, GeometryType::Pointer pGeometry)
Constructor using Geometry.
Definition: fs_generalized_wall_condition.h:167
KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(FSGeneralizedWallCondition)
Pointer definition of FSGeneralizedWallCondition.
void CalculateLeftHandSide(MatrixType &rLeftHandSideMatrix, const ProcessInfo &rCurrentProcessInfo) override
Definition: fs_generalized_wall_condition.h:292
Geometry< NodeType >::GeometriesArrayType GeometriesArrayType
Definition: fs_generalized_wall_condition.h:118
void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Calculate wall stress term for all nodes with SLIP set.
Definition: fs_generalized_wall_condition.h:306
void ApplyIACPenalty(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo)
Apply an IAC penalty term.
Definition: fs_generalized_wall_condition.h:831
void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const override
Provides the global indices for each one of this condition's local rows.
Condition::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override
Create a new FSGeneralizedWallCondition object.
Definition: fs_generalized_wall_condition.h:231
FSGeneralizedWallCondition & operator=(FSGeneralizedWallCondition const &rOther)
Copy constructor.
Definition: fs_generalized_wall_condition.h:199
Properties PropertiesType
Definition: fs_generalized_wall_condition.h:110
ElementPointerType pGetElement()
Definition: fs_generalized_wall_condition.h:504
Matrix ShapeFunctionDerivativesType
Type for a matrix containing the shape function gradients.
Definition: fs_generalized_wall_condition.h:135
std::size_t IndexType
Definition: flags.h:74
bool Is(Flags const &rOther) const
Definition: flags.h:274
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometrical_object.h:248
GeometryType & GetGeometry()
Returns the reference of the geometry.
Definition: geometrical_object.h:158
Geometry base class.
Definition: geometry.h:71
SizeType PointsNumber() const
Definition: geometry.h:528
SizeType size() const
Definition: geometry.h:518
IndexType const & Id() const
Id of this Geometry.
Definition: geometry.h:964
void ShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType &rResult) const
Definition: geometry.h:3708
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
This class defines the node.
Definition: node.h:65
TVariableType::Type & FastGetSolutionStepValue(const TVariableType &rThisVariable)
Definition: node.h:435
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: node.h:466
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Properties encapsulates data shared by different Elements or Conditions. It can store any type of dat...
Definition: properties.h:69
The serialization consists in storing the state of an object into a storage format like data file or ...
Definition: serializer.h:123
Variable class contains all information needed to store and retrive data from a data container.
Definition: variable.h:63
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
#define KRATOS_SERIALIZE_SAVE_BASE_CLASS(Serializer, BaseType)
Definition: define.h:812
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_TRY
Definition: define.h:109
#define KRATOS_SERIALIZE_LOAD_BASE_CLASS(Serializer, BaseType)
Definition: define.h:815
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
static double min(double a, double b)
Definition: GeometryFunctions.h:71
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
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
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
AMatrix::MatrixRow< const TExpressionType > row(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression, std::size_t RowIndex)
Definition: amatrix_interface.h:649
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
rho
Definition: generate_droplet_dynamics.py:86
E
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:26
int C
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:27
S
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:39
tuple tmp
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:98
nu
Definition: isotropic_damage_automatic_differentiation.py:135
R
Definition: isotropic_damage_automatic_differentiation.py:172
tuple Q
Definition: isotropic_damage_automatic_differentiation.py:235
def load(f)
Definition: ode_solve.py:307
int d
Definition: ode_solve.py:397
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
A
Definition: sensitivityMatrix.py:70
N
Definition: sensitivityMatrix.py:29
B
Definition: sensitivityMatrix.py:76
integer i
Definition: TensorModule.f:17
integer function sign1(a)
Definition: TensorModule.f:872