40 #ifndef KRATOS_WALL_CONDITION_H
41 #define KRATOS_WALL_CONDITION_H
96 template<
unsigned int TDim,
unsigned int TNumNodes = TDim >
158 GeometryType::Pointer pGeometry):
170 GeometryType::Pointer pGeometry,
171 PropertiesType::Pointer pProperties):
210 return Kratos::make_intrusive<WallCondition>(NewId,
GetGeometry().
Create(ThisNodes), pProperties);
214 Condition::Pointer
Create(
IndexType NewId, Condition::GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties)
const override
216 return Kratos::make_intrusive<WallCondition>(NewId, pGeom, pProperties);
231 pNewCondition->SetData(this->
GetData());
232 pNewCondition->SetFlags(this->
GetFlags());
234 return pNewCondition;
257 const ProcessInfo& r_process_info = rCurrentProcessInfo;
258 unsigned int step = r_process_info[FRACTIONAL_STEP];
262 const SizeType LocalSize = TDim * TNumNodes;
264 if (rLeftHandSideMatrix.size1() != LocalSize)
265 rLeftHandSideMatrix.
resize(LocalSize,LocalSize,
false);
266 if (rRightHandSideVector.size() != LocalSize)
267 rRightHandSideVector.
resize(LocalSize,
false);
274 this->
ApplyWallLaw(rLeftHandSideMatrix,rRightHandSideVector);
278 if(this->
Is(INTERFACE) &&
step==5 )
281 const double N = 1.0 /
static_cast<double>(TNumNodes);
284 const double Area =
norm_2(rNormal);
286 if (rLeftHandSideMatrix.size1() != TNumNodes)
287 rLeftHandSideMatrix.
resize(TNumNodes,TNumNodes,
false);
288 if (rRightHandSideVector.size() != TNumNodes)
289 rRightHandSideVector.
resize(TNumNodes,
false);
294 const double dt = r_process_info[DELTA_TIME];
295 const double equivalent_structural_density = r_process_info[DENSITY];
296 const double diag_term =
dt*Area*
N/( equivalent_structural_density );
298 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
300 rLeftHandSideMatrix(iNode,iNode) = diag_term;
305 if (rLeftHandSideMatrix.size1() != 0)
306 rLeftHandSideMatrix.
resize(0,0,
false);
308 if (rRightHandSideVector.size() != 0)
309 rRightHandSideVector.
resize(0,
false);
333 if(this->
GetGeometry()[
i].SolutionStepsDataHas(VELOCITY) ==
false)
335 if(this->
GetGeometry()[
i].SolutionStepsDataHas(MESH_VELOCITY) ==
false)
337 if(this->
GetGeometry()[
i].SolutionStepsDataHas(NORMAL) ==
false)
339 if(this->
GetGeometry()[
i].HasDofFor(VELOCITY_X) ==
false ||
340 this->
GetGeometry()[
i].HasDofFor(VELOCITY_Y) ==
false ||
358 const ProcessInfo& rCurrentProcessInfo)
const override;
367 const ProcessInfo& CurrentProcessInfo)
const override;
376 int Step = 0)
const override
378 const SizeType LocalSize = TDim * TNumNodes;
379 unsigned int LocalIndex = 0;
381 if (Values.size() != LocalSize)
382 Values.
resize(LocalSize,
false);
384 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
387 for (
unsigned int d = 0;
d < TDim; ++
d)
388 Values[LocalIndex++] = rVelocity[
d];
422 if (rVariable == NORMAL)
435 rValues[0] = const_this->
GetValue(rVariable);
443 std::vector<double>& rValues,
454 rValues[0] = const_this->
GetValue(rVariable);
465 rValues[0] = const_this->
GetValue(rVariable);
471 std::vector<Vector>& rValues,
476 rValues[0] = const_this->
GetValue(rVariable);
482 std::vector<Matrix>& rValues,
487 rValues[0] = const_this->
GetValue(rVariable);
507 std::string
Info()
const override
509 std::stringstream buffer;
517 rOStream <<
"WallCondition" << TDim <<
"D #" << this->
Id();
565 const unsigned int BlockSize = TDim;
566 const double NodalFactor = 1.0 /
double(TDim);
568 double area = NodalFactor * rGeometry.
DomainSize();
571 for(
unsigned int itNode = 0; itNode < rGeometry.
PointsNumber(); ++itNode)
573 const NodeType& rConstNode = rGeometry[itNode];
574 const double y = rConstNode.
GetValue(Y_WALL);
575 if(
y > 0.0 && rConstNode.
Is(SLIP) )
578 const array_1d<double,3>& VelMesh = rGeometry[itNode].FastGetSolutionStepValue(MESH_VELOCITY);
581 const double Ikappa = 1.0/0.41;
582 const double B = 5.2;
583 const double limit_yplus = 10.9931899;
585 const double rho = rGeometry[itNode].FastGetSolutionStepValue(DENSITY);
586 const double nu = rGeometry[itNode].FastGetSolutionStepValue(VISCOSITY);
588 double wall_vel = 0.0;
589 for (
size_t d = 0;
d < TDim;
d++)
591 wall_vel += Vel[
d]*Vel[
d];
593 wall_vel = sqrt(wall_vel);
595 if (wall_vel > 1
e-12)
599 double utau = sqrt(wall_vel *
nu /
y);
600 double yplus =
y * utau /
nu;
603 if (yplus > limit_yplus)
610 unsigned int iter = 0;
612 const double tol = 1
e-6;
613 double uplus = Ikappa * log(yplus) +
B;
615 while(iter < 100 && fabs(dx) >
tol * utau)
618 double f = utau * uplus - wall_vel;
619 double df = uplus + Ikappa;
624 yplus =
y * utau /
nu;
625 uplus = Ikappa * log(yplus) +
B;
630 std::cout <<
"WARNING: wall condition Newton-Raphson did not converge. Residual is " << dx << std::endl;
634 const double Tmp = area * utau * utau *
rho / wall_vel;
635 for (
unsigned int d = 0;
d < TDim;
d++)
637 unsigned int k = itNode*BlockSize+
d;
638 rLocalVector[
k] -= Vel[
d] * Tmp;
639 rLocalMatrix(
k,
k) += Tmp;
691 void save(
Serializer& rSerializer)
const override
743 template<
unsigned int TDim,
unsigned int TNumNodes >
751 template<
unsigned int TDim,
unsigned int TNumNodes >
756 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
PropertiesType::Pointer pGetProperties()
returns the pointer to the property of the condition. Does not throw an error, to allow copying of co...
Definition: condition.h:964
std::size_t IndexType
Definition: flags.h:74
bool Is(Flags const &rOther) const
Definition: flags.h:274
GeometryType::Pointer pGetGeometry()
Returns the pointer to the geometry.
Definition: geometrical_object.h:140
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometrical_object.h:248
Flags & GetFlags()
Returns the flags of the object.
Definition: geometrical_object.h:176
GeometryType & GetGeometry()
Returns the reference of the geometry.
Definition: geometrical_object.h:158
DataValueContainer & GetData()
Definition: geometrical_object.h:212
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
virtual double DomainSize() const
This method calculate and return length, area or volume of this geometry depending to it's dimension.
Definition: geometry.h:1371
This object defines an indexed object.
Definition: indexed_object.h:54
IndexType Id() const
Definition: indexed_object.h:107
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 & 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
A sorted associative container similar to an STL set, but uses a vector to store pointers to its data...
Definition: pointer_vector_set.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
Implements a wall condition for the monolithic formulation.
Definition: wall_condition.h:98
int Check(const ProcessInfo &rCurrentProcessInfo) const override
Check that all data required by this condition is available and reasonable.
Definition: wall_condition.h:316
std::vector< std::size_t > EquationIdVectorType
Definition: wall_condition.h:122
void CalculateLeftHandSide(MatrixType &rLeftHandSideMatrix, const ProcessInfo &rCurrentProcessInfo) override
Definition: wall_condition.h:238
KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(WallCondition)
Pointer definition of WallCondition.
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: wall_condition.h:521
void ApplyWallLaw(MatrixType &rLocalMatrix, VectorType &rLocalVector)
Commpute the wall stress and add corresponding terms to the system contributions.
Definition: wall_condition.h:561
void CalculateOnIntegrationPoints(const Variable< array_1d< double, 6 > > &rVariable, std::vector< array_1d< double, 6 > > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: wall_condition.h:458
void GetValuesVector(Vector &Values, int Step=0) const override
Returns VELOCITY_X, VELOCITY_Y, (VELOCITY_Z,) for each node.
Definition: wall_condition.h:375
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: wall_condition.h:515
WallCondition(IndexType NewId, GeometryType::Pointer pGeometry)
Constructor using Geometry.
Definition: wall_condition.h:157
std::size_t IndexType
Definition: wall_condition.h:118
void CalculateOnIntegrationPoints(const Variable< Matrix > &rVariable, std::vector< Matrix > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: wall_condition.h:480
Geometry< NodeType > GeometryType
Definition: wall_condition.h:110
void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const override
Provides the global indices for each one of this element's local rows.
Vector VectorType
Definition: wall_condition.h:114
void CalculateNormal(array_1d< double, 3 > &An)
std::vector< Dof< double >::Pointer > DofsVectorType
Definition: wall_condition.h:124
WallCondition(IndexType NewId=0)
Default constructor.
Definition: wall_condition.h:136
std::size_t SizeType
Definition: wall_condition.h:120
WallCondition & operator=(WallCondition const &rOther)
Copy constructor.
Definition: wall_condition.h:191
void CalculateOnIntegrationPoints(const Variable< double > &rVariable, std::vector< double > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: wall_condition.h:441
WallCondition(IndexType NewId, const NodesArrayType &ThisNodes)
Constructor using an array of nodes.
Definition: wall_condition.h:146
PointerVectorSet< Dof< double >, IndexedObject > DofsArrayType
Definition: wall_condition.h:126
Geometry< NodeType >::PointsArrayType NodesArrayType
Definition: wall_condition.h:112
Condition::Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const override
Create a new WallCondition object.
Definition: wall_condition.h:208
Matrix MatrixType
Definition: wall_condition.h:116
WallCondition(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
Constructor using Properties.
Definition: wall_condition.h:169
Condition::Pointer Clone(IndexType NewId, NodesArrayType const &rThisNodes) const override
Definition: wall_condition.h:227
Properties PropertiesType
Definition: wall_condition.h:108
Node NodeType
Definition: wall_condition.h:106
WallCondition(WallCondition const &rOther)
Copy constructor.
Definition: wall_condition.h:177
~WallCondition() override
Destructor.
Definition: wall_condition.h:183
std::string Info() const override
Turn back information as a string.
Definition: wall_condition.h:507
void CalculateOnIntegrationPoints(const Variable< Vector > &rVariable, std::vector< Vector > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: wall_condition.h:469
Condition::Pointer Create(IndexType NewId, Condition::GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override
Definition: wall_condition.h:214
void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Calculate wall stress term for all nodes with SLIP set.
Definition: wall_condition.h:253
void CalculateOnIntegrationPoints(const Variable< array_1d< double, 3 > > &rVariable, std::vector< array_1d< double, 3 > > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: wall_condition.h:416
void ApplyNeumannCondition(MatrixType &rLocalMatrix, VectorType &rLocalVector)
Apply boundary terms to allow imposing a pressure (normal stress), with a correction to prevent inflo...
Definition: wall_condition.cpp:227
void GetDofList(DofsVectorType &ConditionDofList, const ProcessInfo &CurrentProcessInfo) const override
Returns a list of the element's Dofs.
#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
dt
Definition: DEM_benchmarks.py:173
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
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
TABLE_NUMBER_ANGULAR_VELOCITY TABLE_NUMBER_MOMENT I33 BEAM_INERTIA_ROT_UNIT_LENGHT_Y KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BEAM_INERTIA_ROT_UNIT_LENGHT_Z) typedef std double
Definition: DEM_application_variables.h:182
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
int step
Definition: face_heat.py:88
y
Other simbols definition.
Definition: generate_axisymmetric_navier_stokes_element.py:54
f
Definition: generate_convection_diffusion_explicit_element.py:112
rho
Definition: generate_droplet_dynamics.py:86
int tol
Definition: hinsberg_optimization.py:138
nu
Definition: isotropic_damage_automatic_differentiation.py:135
def load(f)
Definition: ode_solve.py:307
int d
Definition: ode_solve.py:397
int k
Definition: quadrature.py:595
N
Definition: sensitivityMatrix.py:29
B
Definition: sensitivityMatrix.py:76
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31