40 #ifndef KRATOS_FS_GENERALIZED_WALL_CONDITION_H
41 #define KRATOS_FS_GENERALIZED_WALL_CONDITION_H
98 template<
unsigned int TDim,
unsigned int TNumNodes = TDim>
158 :
Condition(NewId, ThisNodes), mInitializeWasPerformed(false), mpElement()
168 :
Condition(NewId, pGeometry), mInitializeWasPerformed(false), mpElement()
179 PropertiesType::Pointer pProperties)
180 :
Condition(NewId, pGeometry, pProperties), mInitializeWasPerformed(false), mpElement()
186 mInitializeWasPerformed(rOther.mInitializeWasPerformed), mpElement(rOther.mpElement)
202 mInitializeWasPerformed = rOther.mInitializeWasPerformed;
203 mpElement = rOther.mpElement;
220 PropertiesType::Pointer pProperties)
const override
222 return Kratos::make_intrusive<FSGeneralizedWallCondition>(NewId,
GetGeometry().
Create(ThisNodes), pProperties);
233 GeometryType::Pointer pGeom,
234 PropertiesType::Pointer pProperties)
const override
236 return Kratos::make_intrusive<FSGeneralizedWallCondition>(NewId, pGeom, pProperties);
250 if (mInitializeWasPerformed)
255 mInitializeWasPerformed =
true;
261 mpElement = this->
GetValue(NEIGHBOUR_ELEMENTS)(0);
264 Edge = rElemGeom[1].Coordinates() - rElemGeom[0].Coordinates();
265 mMinEdgeLength = Edge[0] * Edge[0];
268 mMinEdgeLength += Edge[
d] * Edge[
d];
275 Edge = rElemGeom[
j].Coordinates() - rElemGeom[
k].Coordinates();
276 EdgeLength = Edge[0] * Edge[0];
280 EdgeLength += Edge[
d] * Edge[
d];
283 mMinEdgeLength = (EdgeLength < mMinEdgeLength) ? EdgeLength : mMinEdgeLength;
286 mMinEdgeLength = sqrt(mMinEdgeLength);
312 if (mInitializeWasPerformed ==
false)
317 if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1)
320 const SizeType LocalSize = TDim * TNumNodes;
322 if (rLeftHandSideMatrix.size1() != LocalSize)
324 rLeftHandSideMatrix.
resize(LocalSize, LocalSize);
326 if (rRightHandSideVector.size() != LocalSize)
328 rRightHandSideVector.
resize(LocalSize);
335 this->
ApplyWallLaw(rLeftHandSideMatrix, rRightHandSideVector);
337 else if (rCurrentProcessInfo[FRACTIONAL_STEP] == 5)
340 const SizeType LocalSize = TNumNodes;
342 if (rLeftHandSideMatrix.size1() != LocalSize)
344 rLeftHandSideMatrix.
resize(LocalSize, LocalSize);
346 if (rRightHandSideVector.size() != LocalSize)
348 rRightHandSideVector.
resize(LocalSize);
354 if (this->
Is(INTERFACE))
355 this->
ApplyIACPenalty(rLeftHandSideMatrix, rRightHandSideVector, rCurrentProcessInfo);
359 if (rLeftHandSideMatrix.size1() != 0)
361 rLeftHandSideMatrix.
resize(0,0,
false);
364 if (rRightHandSideVector.size() != 0)
366 rRightHandSideVector.
resize(0,
false);
389 if(this->
GetGeometry()[
i].SolutionStepsDataHas(VELOCITY) ==
false)
391 if(this->
GetGeometry()[
i].SolutionStepsDataHas(MESH_VELOCITY) ==
false)
393 if(this->
GetGeometry()[
i].SolutionStepsDataHas(DENSITY) ==
false)
395 if(this->
GetGeometry()[
i].SolutionStepsDataHas(VISCOSITY) ==
false)
397 if(this->
GetGeometry()[
i].SolutionStepsDataHas(NORMAL) ==
false)
399 if(this->
GetGeometry()[
i].HasDofFor(VELOCITY_X) ==
false ||
400 this->
GetGeometry()[
i].HasDofFor(VELOCITY_Y) ==
false ||
419 const ProcessInfo& rCurrentProcessInfo)
const override;
427 const ProcessInfo& rCurrentProcessInfo)
const override;
436 const SizeType LocalSize = TDim * TNumNodes;
437 unsigned int LocalIndex = 0;
439 if (Values.size() != LocalSize)
441 Values.
resize(LocalSize,
false);
444 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
447 for (
unsigned int d = 0;
d < TDim; ++
d)
449 Values[LocalIndex++] = rVelocity[
d];
467 std::string
Info()
const override
469 std::stringstream buffer;
470 buffer <<
"FSGeneralizedWallCondition" << TDim <<
"D";
476 { rOStream <<
"FSGeneralizedWallCondition";}
506 return mpElement->shared_from_this();
509 template<
class TVariableType >
516 rResult = rShapeFunc[0] * rGeom[0].FastGetSolutionStepValue(Var);
520 rResult += rShapeFunc[
i] * rGeom[
i].FastGetSolutionStepValue(Var);
535 const double& pres = rElemGeom[0].FastGetSolutionStepValue(PRESSURE,1);
538 rResult[
d] = rDN_DX(0,
d) * pres;
543 const double& pres = rElemGeom[
i].FastGetSolutionStepValue(PRESSURE,1);
546 rResult[
d] += rDN_DX(
i,
d) * pres;
570 double rho,
nu, func1, func2, Vel1, Vel2, Vel12, YPlus1, YPlus2,
sign1, sign2;
573 Vel1 = sqrt(fabs(rWallStress) /
rho);
574 Vel2 = pow(
nu * fabs(rWallGradP) /
rho, 0.333333);
582 YPlus1 = Vel1 * rWallHeight /
nu;
583 YPlus2 = Vel2 * rWallHeight /
nu;
589 func1 = (1.0 + (0.01 - 0.0029 * YPlus1) * YPlus1) * YPlus1;
591 else if (YPlus1 <= 30)
593 func1 = -0.872 + (1.465 + (-0.0702 + (0.00166 - 0.00001495 * YPlus1) * YPlus1) * YPlus1) * YPlus1;
595 else if (YPlus1 <= 140)
597 func1 = 8.6 + (0.1864 + (-0.002006 + (0.00001144 - 0.00000002551 * YPlus1) * YPlus1) * YPlus1) * YPlus1;
601 func1 = 2.439 * log(YPlus1) + 5.0;
608 func2 = (0.5 - 0.00731 * YPlus2) * YPlus2 * YPlus2;
610 else if (YPlus2 <= 15)
612 func2 = -15.138 + (8.4688 + (-0.81976 + (0.037292 - 0.00063866 * YPlus2) * YPlus2) * YPlus2) * YPlus2;
614 else if (YPlus2 <= 30)
616 func2 = 11.925 + (0.934 + (-0.027805 + (0.00046262 - 0.0000031442 * YPlus2) * YPlus2) * YPlus2) * YPlus2;
620 func2 = 5.0 * log(YPlus2) + 8.0;
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;
639 const double Eps = 1.0e-08;
642 double C, Xm, P,
Q,
R,
S, Min1, Min2, ResA, ResB, ResC, Tol;
653 if (ResB * ResA <= 0.0)
659 if (ResB * ResA > 0.0)
671 if (ResB * ResC > 0.0)
679 if (fabs(ResC) < fabs(ResB))
689 Tol = 2.0 * Eps * fabs(
B) + 0.5e-10;
692 if (fabs(Xm) <= Tol || ResB == 0.0)
697 if (fabs(
E) >= Tol && fabs(ResA) > fabs(ResB))
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);
720 Min1 = 3.0 * Xm *
Q - fabs(Tol *
Q);
775 const double Area =
norm_2(rNormal);
779 const array_1d<double, 3>& rNodeNormal = rGeometry[iNode].FastGetSolutionStepValue(NORMAL);
780 if (
inner_prod(rNormal, rNodeNormal) < AngleCosine * Area *
norm_2(rNodeNormal))
793 const unsigned int BlockSize = TDim;
794 double WallHeight, Area, WallVelMag,
tmp, WallStress, WallGradP, WallForce;
799 WallVelMag =
norm_2(WallVel);
804 WallForce = (Area /
static_cast<double>(TDim)) * WallStress;
809 if(rNode.
GetValue(Y_WALL) != 0.0 && rNode.
Is(SLIP))
813 WallVel /= (
tmp != 0.0) ?
tmp : 1.0;
815 for (
unsigned int d=0;
d < TDim;
d++)
817 unsigned int k =
i * BlockSize +
d;
818 rLocalVector[
k] -= WallVel[
d] * WallForce;
837 const double Area =
norm_2(rNormal);
838 const double DiagonalTerm = Area /
static_cast<double>(TNumNodes)
839 / (rCurrentProcessInfo[DENSITY] * rCurrentProcessInfo[BDF_COEFFICIENTS][0]);
843 rLeftHandSideMatrix(iNode,iNode) += DiagonalTerm;
868 bool mInitializeWasPerformed;
869 double mMinEdgeLength;
878 void save(
Serializer& rSerializer)
const override
922 template<
unsigned int TDim,
unsigned int TNumNodes>
930 template<
unsigned int TDim,
unsigned int TNumNodes>
935 rOStream << std::endl;
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