48 #ifndef MONOLITHIC_DEM_COUPLED_H
49 #define MONOLITHIC_DEM_COUPLED_H
64 #include "utilities/geometry_utilities.h"
134 template<
unsigned int TDim,
135 unsigned int TNumNodes = TDim + 1 >
226 Element(NewId, pGeometry, pProperties)
252 PropertiesType::Pointer pProperties)
const override
273 if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1) {
274 const unsigned int LocalSize = (TDim + 1) * TNumNodes;
277 if (rLeftHandSideMatrix.size1() != LocalSize)
278 rLeftHandSideMatrix.
resize(LocalSize, LocalSize,
false);
283 this->CalculateRightHandSide(rRightHandSideVector, rCurrentProcessInfo);
287 const unsigned int LocalSize = TDim * TNumNodes;
289 if (rLeftHandSideMatrix.size1() != LocalSize)
290 rLeftHandSideMatrix.
resize(LocalSize, LocalSize,
false);
293 CalculateLaplacianMassMatrix(rLeftHandSideMatrix, rCurrentProcessInfo);
295 this->CalculateRightHandSide(rRightHandSideVector, rCurrentProcessInfo);
306 if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1) {
308 const unsigned int LocalSize = (TDim + 1) * TNumNodes;
310 if (rLeftHandSideMatrix.size1() != LocalSize) rLeftHandSideMatrix.
resize(LocalSize, LocalSize,
false);
316 const unsigned int LocalSize = TDim * TNumNodes;
318 if (rLeftHandSideMatrix.size1() != LocalSize) rLeftHandSideMatrix.
resize(LocalSize, LocalSize,
false);
321 CalculateLaplacianMassMatrix(rLeftHandSideMatrix, rCurrentProcessInfo);
346 this->EvaluateInPoint(Density, DENSITY,
N);
348 if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1) {
350 const unsigned int LocalSize = (TDim + 1) * TNumNodes;
353 if (rRightHandSideVector.size() != LocalSize)
354 rRightHandSideVector.
resize(LocalSize,
false);
359 this->AddMomentumRHS(rRightHandSideVector, Density,
N, Area);
361 const double& DeltaTime = rCurrentProcessInfo[DELTA_TIME];
362 static const double arr[] = {1.0,-1.0};
363 std::vector<double> SchemeWeights (arr, arr +
sizeof(arr) /
sizeof(arr[0]));
364 this->AddMassRHS(rRightHandSideVector, Density,
N, Area, SchemeWeights, DeltaTime);
368 const unsigned int LocalSize = TDim * TNumNodes;
371 if (rRightHandSideVector.size() != LocalSize)
372 rRightHandSideVector.
resize(LocalSize,
false);
376 this->AddRHSLaplacian(rRightHandSideVector, DN_DX, Area);
399 if (rCurrentProcessInfo[OSS_SWITCH] == 1)
402 this->GetAdvectiveVel(AdvVel,
N);
405 this->EvaluateInPoint(KinViscosity, VISCOSITY,
N);
408 this->GetEffectiveViscosity(Density,KinViscosity,
N, DN_DX, Viscosity, rCurrentProcessInfo);
411 double TauOne, TauTwo;
412 this->CalculateTau(TauOne, TauTwo, AdvVel, Area,Density, Viscosity, rCurrentProcessInfo);
414 this->AddProjectionToRHS(rRightHandSideVector, AdvVel, Density, TauOne, TauTwo,
N, DN_DX, Area,rCurrentProcessInfo[DELTA_TIME]);
470 const unsigned int LocalSize = (TDim + 1) * TNumNodes;
473 if (rMassMatrix.size1() != LocalSize)
474 rMassMatrix.
resize(LocalSize, LocalSize,
false);
476 rMassMatrix =
ZeroMatrix(LocalSize, LocalSize);
486 this->EvaluateInPoint(Density, DENSITY,
N);
489 double Coeff = Density * Area / TNumNodes;
491 this->CalculateLumpedMassMatrix(rMassMatrix, Coeff);
510 if (rCurrentProcessInfo[OSS_SWITCH] != 1)
513 this->EvaluateInPoint(KinViscosity, VISCOSITY,
N);
516 this->GetEffectiveViscosity(Density,KinViscosity,
N, DN_DX, Viscosity, rCurrentProcessInfo);
520 this->GetAdvectiveVel(AdvVel,
N);
523 double TauOne, TauTwo;
524 this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
536 this->AddMassStabTerms(rMassMatrix, Density, AdvVel, TauOne,
N, DN_DX, Area);
544 const unsigned int LocalSize = TDim * TNumNodes;
547 if (rMassMatrix.size1() != LocalSize)
548 rMassMatrix.
resize(LocalSize, LocalSize,
false);
550 rMassMatrix =
ZeroMatrix(LocalSize, LocalSize);
561 double Coeff = Area / TNumNodes;
563 this->CalculateLaplacianLumpedMassMatrix(rMassMatrix, Coeff);
580 const unsigned int LocalSize = (TDim + 1) * TNumNodes;
584 if (rDampingMatrix.size1() != LocalSize)
585 rDampingMatrix.
resize(LocalSize, LocalSize,
false);
596 double Density, KinViscosity;
597 this->EvaluateInPoint(Density, DENSITY,
N);
598 this->EvaluateInPoint(KinViscosity, VISCOSITY,
N);
601 this->GetEffectiveViscosity(Density,KinViscosity,
N, DN_DX, Viscosity, rCurrentProcessInfo);
605 this->GetAdvectiveVel(AdvVel,
N);
608 double TauOne, TauTwo;
609 this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
619 this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo,
N, DN_DX, Area);
647 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
650 for (
unsigned int d = 0;
d < TDim; ++
d)
652 U[LocalIndex] = rVel[
d];
655 U[LocalIndex] = this->GetGeometry()[iNode].FastGetSolutionStepValue(PRESSURE);
659 noalias(rRightHandSideVector) -=
prod(rDampingMatrix, U);
708 for(
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
710 this->GetGeometry()[iNode].SetLock();
711 this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION_OLD) = this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION);
712 this->GetGeometry()[iNode].UnSetLock();
733 if (rVariable == ERROR_RATIO)
742 double Density, KinViscosity;
743 this->EvaluateInPoint(Density, DENSITY,
N);
744 this->EvaluateInPoint(KinViscosity, VISCOSITY,
N);
747 this->GetEffectiveViscosity(Density,KinViscosity,
N, DN_DX, Viscosity, rCurrentProcessInfo);
751 this->GetAdvectiveVel(AdvVel,
N);
758 this->CalculateStaticTau(TauOne, AdvVel, Area,Density, Viscosity);
760 if ( rCurrentProcessInfo[OSS_SWITCH] != 1 )
763 this->ASGSMomResidual(AdvVel,Density,ElementalMomRes,
N,DN_DX,1.0);
782 ElementalMomRes *= TauOne;
786 this->OSSMomResidual(AdvVel,Density,ElementalMomRes,
N,DN_DX,1.0);;
787 ElementalMomRes *= TauOne;
791 double ErrorRatio(0.0);
795 for (
unsigned int i = 0;
i < TDim; ++
i)
797 ErrorRatio += ElementalMomRes[
i] * ElementalMomRes[
i];
800 ErrorRatio = sqrt(ErrorRatio);
801 ErrorRatio /= Density;
802 this->
SetValue(ERROR_RATIO, ErrorRatio);
803 rOutput = ErrorRatio;
805 else if (rVariable == NODAL_AREA)
814 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
816 this->GetGeometry()[
i].SetLock();
817 this->GetGeometry()[
i].FastGetSolutionStepValue(NODAL_AREA) += Area *
N[
i];
818 this->GetGeometry()[
i].UnSetLock();
839 if (rVariable == ADVPROJ)
849 this->EvaluateInPoint(Density, DENSITY,
N);
853 this->GetAdvectiveVel(AdvVel,
N);
857 double ElementalMassRes(0);
859 this->AddProjectionResidualContribution(AdvVel, Density, ElementalMomRes, ElementalMassRes,
N, DN_DX, Area);
861 if (rCurrentProcessInfo[OSS_SWITCH] == 1)
864 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
866 this->GetGeometry()[
i].SetLock();
868 for (
unsigned int d = 0;
d < TDim; ++
d)
869 rAdvProj[
d] +=
N[
i] * ElementalMomRes[
d];
871 this->GetGeometry()[
i].FastGetSolutionStepValue(DIVPROJ) +=
N[
i] * ElementalMassRes;
872 this->GetGeometry()[
i].FastGetSolutionStepValue(NODAL_AREA) += Area *
N[
i];
873 this->GetGeometry()[
i].UnSetLock();
878 rOutput = ElementalMomRes;
880 else if (rVariable == SUBSCALE_VELOCITY)
890 this->EvaluateInPoint(Density, DENSITY,
N);
894 this->GetAdvectiveVel(AdvVel,
N);
898 double ElementalMassRes(0.0);
900 this->AddProjectionResidualContribution(AdvVel, Density, ElementalMomRes, ElementalMassRes,
N, DN_DX, Area);
902 if (rCurrentProcessInfo[OSS_SWITCH] == 1)
907 const double Weight = ConsistentMassCoef(Area);
909 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
911 this->GetGeometry()[
i].SetLock();
915 double& rMassRHS = this->GetGeometry()[
i].GetValue(DIVPROJ);
916 for (
unsigned int d = 0;
d < TDim; ++
d)
917 rMomRHS[
d] +=
N[
i] * ElementalMomRes[
d];
919 rMassRHS +=
N[
i] * ElementalMassRes;
922 this->GetGeometry()[
i].FastGetSolutionStepValue(NODAL_AREA) += Area *
N[
i];
925 for(
unsigned int j = 0;
j < TNumNodes; ++
j)
927 for(
unsigned int d = 0;
d < TDim; ++
d)
928 rMomRHS[
d] -= Weight * this->GetGeometry()[
j].FastGetSolutionStepValue(ADVPROJ)[
d];
929 rMassRHS -= Weight * this->GetGeometry()[
j].FastGetSolutionStepValue(DIVPROJ);
931 for(
unsigned int d = 0;
d < TDim; ++
d)
932 rMomRHS[
d] -= Weight * this->GetGeometry()[
i].FastGetSolutionStepValue(ADVPROJ)[
d];
933 rMassRHS -= Weight * this->GetGeometry()[
i].FastGetSolutionStepValue(DIVPROJ);
935 this->GetGeometry()[
i].UnSetLock();
940 rOutput = ElementalMomRes;
1002 std::vector<double>& rValues,
1005 if (rVariable == TAUONE || rVariable == TAUTWO || rVariable == MU)
1007 double TauOne, TauTwo;
1014 this->GetAdvectiveVel(AdvVel,
N);
1016 double Density,KinViscosity;
1017 this->EvaluateInPoint(Density, DENSITY,
N);
1018 this->EvaluateInPoint(KinViscosity, VISCOSITY,
N);
1021 this->GetEffectiveViscosity(Density,KinViscosity,
N, DN_DX, Viscosity, rCurrentProcessInfo);
1023 this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
1025 rValues.resize(1,
false);
1026 if (rVariable == TAUONE)
1028 rValues[0] = TauOne;
1030 else if (rVariable == TAUTWO)
1032 rValues[0] = TauTwo;
1034 else if (rVariable == MU)
1036 rValues[0] = Density * Viscosity;
1039 else if(rVariable == SUBSCALE_PRESSURE)
1041 double TauOne, TauTwo;
1048 this->GetAdvectiveVel(AdvVel,
N);
1050 double Density,KinViscosity;
1051 this->EvaluateInPoint(Density, DENSITY,
N);
1052 this->EvaluateInPoint(KinViscosity, VISCOSITY,
N);
1055 this->GetEffectiveViscosity(Density,KinViscosity,
N, DN_DX, Viscosity, rCurrentProcessInfo);
1057 this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
1060 for(
unsigned int i=0;
i < TNumNodes;
i++)
1062 for(
unsigned int d = 0;
d < TDim;
d++)
1063 DivU -= DN_DX(
i,
d) * this->GetGeometry()[
i].FastGetSolutionStepValue(VELOCITY)[
d];
1066 rValues.
resize(1,
false);
1068 rValues[0] = TauTwo * DivU;
1070 if(rCurrentProcessInfo[OSS_SWITCH]==1)
1073 for(
unsigned int i=0;
i < TNumNodes;
i++)
1075 Proj +=
N[
i]*this->GetGeometry()[
i].FastGetSolutionStepValue(DIVPROJ);
1077 rValues[0] -= TauTwo*Proj;
1080 else if (rVariable == NODAL_AREA && TDim == 3)
1088 J(0,0) =
X1[0]-X0[0];
1089 J(0,1) =
X2[0]-X0[0];
1090 J(0,2) = X3[0]-X0[0];
1091 J(1,0) =
X1[1]-X0[1];
1092 J(1,1) =
X2[1]-X0[1];
1093 J(1,2) = X3[1]-X0[1];
1094 J(2,0) =
X1[2]-X0[2];
1095 J(2,1) =
X2[2]-X0[2];
1096 J(2,2) = X3[2]-X0[2];
1098 double DetJ =
J(0,0)*(
J(1,1)*
J(2,2) -
J(1,2)*
J(2,1) ) +
J(0,1)*(
J(1,2)*
J(2,0) -
J(1,0)*
J(2,2) ) +
J(0,2)*(
J(1,0)*
J(2,1) -
J(1,1)*
J(2,0) );
1099 rValues.resize(1,
false);
1104 rValues.resize(1,
false);
1112 rValues[0] = const_this->
GetValue(rVariable);
1122 std::vector<Vector>& rValues,
1127 std::vector<Matrix>& rValues,
1154 if(ErrorCode != 0)
return ErrorCode;
1162 for(
unsigned int i=0;
i<this->GetGeometry().size(); ++
i) {
1176 if (this->GetGeometry().WorkingSpaceDimension() == 2) {
1177 for (
unsigned int i=0;
i<this->GetGeometry().size(); ++
i) {
1178 KRATOS_ERROR_IF(this->GetGeometry()[
i].
Z() != 0.0) <<
"Node with non-zero Z coordinate found. Id: "<< this->GetGeometry()[
i].Id() << std::endl;
1197 virtual std::string
Info()
const override
1199 std::stringstream buffer;
1200 buffer <<
"MonolithicDEMCoupled #" << Id();
1201 return buffer.str();
1205 virtual void PrintInfo(std::ostream& rOStream)
const override
1207 rOStream <<
"MonolithicDEMCoupled" << TDim <<
"D";
1260 const double Density,
1261 const double KinViscosity,
1265 double AdvVelNorm = 0.0;
1266 for (
unsigned int d = 0;
d < TDim; ++
d)
1267 AdvVelNorm += rAdvVel[
d] * rAdvVel[
d];
1269 AdvVelNorm = sqrt(AdvVelNorm);
1271 const double Element_Size = this->ElementSize(Area);
1274 TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 5.6666666666 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
1276 TauTwo = Density * (KinViscosity + 0.5 * Element_Size * AdvVelNorm);
1297 const double Density,
1298 const double KinViscosity)
1301 double AdvVelNorm = 0.0;
1302 for (
unsigned int d = 0;
d < TDim; ++
d)
1303 AdvVelNorm += rAdvVel[
d] * rAdvVel[
d];
1305 AdvVelNorm = sqrt(AdvVelNorm);
1307 const double Element_Size = this->ElementSize(Area);
1309 TauOne = 1.0 / (Density*(4.0 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size));
1314 const double Density,
1316 const double Weight)
1318 double Coef = Density * Weight;
1321 this->EvaluateInPoint(BodyForce, BODY_FORCE, rShapeFunc);
1326 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1328 for (
unsigned int d = 0;
d < TDim; ++
d)
1330 F[LocalIndex++] += Coef * rShapeFunc[iNode] * BodyForce[
d];
1338 const double Weight)
1340 double Coef = Weight;
1344 int LocalNodalIndex = 0;
1346 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1348 Velocity = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY);
1349 for (
unsigned int d = 0;
d < TDim; ++
d)
1351 F[LocalIndex++] -= Coef * rShapeDeriv(LocalNodalIndex,
d) * Velocity[
d] * rShapeDeriv(iNode,
d);
1358 const double Density,
1360 const double Weight,
1361 const std::vector<double>& TimeSchemeWeights,
1362 const double& DeltaTime)
1364 double FluidFractionRate = 0.0;
1365 this->EvaluateTimeDerivativeInPoint(FluidFractionRate, FLUID_FRACTION_RATE, rShapeFunc, DeltaTime, TimeSchemeWeights);
1368 int LocalIndex = TDim;
1370 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
1371 F[LocalIndex] -= Weight * rShapeFunc[iNode] * FluidFractionRate;
1372 LocalIndex += TDim + 1;
1381 const double Density,
1382 const double TauOne,
1383 const double TauTwo,
1386 const double Weight,
1387 const double DeltaTime = 1.0)
1389 const unsigned int BlockSize = TDim + 1;
1392 this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv);
1394 double DivProj = 0.0;
1395 this->EvaluateInPoint(MomProj,ADVPROJ,rShapeFunc);
1396 this->EvaluateInPoint(DivProj,DIVPROJ,rShapeFunc);
1401 unsigned int FirstRow = 0;
1403 for (
unsigned int i = 0;
i < TNumNodes;
i++)
1406 double FluidFraction = this->GetGeometry()[
i].FastGetSolutionStepValue(FLUID_FRACTION);
1408 FluidFractionGradient[0] += rShapeDeriv(
i, 0) * FluidFraction;
1409 FluidFractionGradient[1] += rShapeDeriv(
i, 1) * FluidFraction;
1410 FluidFractionGradient[2] += rShapeDeriv(
i, 2) * FluidFraction;
1413 for (
unsigned int d = 0;
d < TDim;
d++)
1417 RHS[FirstRow+
d] -= Weight * (Density * AGradN[
i] * MomProj[
d] + (FluidFraction * rShapeDeriv(
i,
d) + FluidFractionGradient[
d] * rShapeFunc[
i]) * DivProj);
1419 RHS[FirstRow+TDim] -= Weight * rShapeDeriv(
i,
d) * MomProj[
d];
1421 FirstRow += BlockSize;
1435 unsigned int DofIndex = 0;
1436 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1438 for (
unsigned int d = 0;
d < TDim; ++
d)
1440 rLHSMatrix(DofIndex, DofIndex) += Mass;
1451 unsigned int DofIndex = 0;
1452 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1454 for (
unsigned int d = 0;
d < TDim; ++
d)
1456 rLHSMatrix(DofIndex, DofIndex) += Mass;
1465 const double Density,
1466 const double Weight)
1468 const unsigned int BlockSize = TDim + 1;
1470 double Coef = Density * Weight;
1471 unsigned int FirstRow(0), FirstCol(0);
1475 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
1478 for (
unsigned int j = 0;
j < TNumNodes; ++
j)
1481 K = Coef * rShapeFunc[
i] * rShapeFunc[
j];
1483 for (
unsigned int d = 0;
d < TDim; ++
d)
1485 rLHSMatrix(FirstRow +
d, FirstCol +
d) +=
K;
1488 FirstCol += BlockSize;
1491 FirstRow += BlockSize;
1511 const double Density,
1513 const double TauOne,
1516 const double Weight)
1518 const unsigned int BlockSize = TDim + 1;
1520 double Coef = Weight * TauOne;
1521 unsigned int FirstRow(0), FirstCol(0);
1526 this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv);
1528 double AdvVelDiv = 0.0;
1529 this->GetAdvectiveVelDivergence(AdvVelDiv, rShapeDeriv);
1531 double FluidFraction;
1532 this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
1536 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
1539 for (
unsigned int j = 0;
j < TNumNodes; ++
j)
1543 K = Coef * Density * AGradN[
i] * Density * rShapeFunc[
j];
1547 for (
unsigned int d = 0;
d < TDim; ++
d)
1549 rLHSMatrix(FirstRow +
d, FirstCol +
d) +=
K;
1553 rLHSMatrix(FirstRow + TDim, FirstCol +
d) += Coef * Density * FluidFraction * rShapeDeriv(
i,
d) * rShapeFunc[
j];
1557 FirstCol += BlockSize;
1560 FirstRow += BlockSize;
1568 const double Density,
1569 const double Viscosity,
1571 const double TauOne,
1572 const double TauTwo,
1575 const double Weight)
1577 const unsigned int BlockSize = TDim + 1;
1581 this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv);
1583 double AdvVelDiv = 0.0;
1584 this->GetAdvectiveVelDivergence(AdvVelDiv, rShapeDeriv);
1585 this->GetModifiedConvectionOperator(AGradNMod, rAdvVel, AdvVelDiv, rShapeFunc, rShapeDeriv);
1588 unsigned int FirstRow(0), FirstCol(0);
1589 double K, G, PDivV,
L, qF;
1592 this->EvaluateInPoint(BodyForce,BODY_FORCE,rShapeFunc);
1593 BodyForce *= Density;
1595 double PAlphaDivV, GAlpha, FluidFraction, FluidFractionRate;
1597 this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
1598 this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
1600 for (
unsigned int i = 0;
i < TNumNodes; ++
i) {
1601 this->GetGeometry()[
i].FastGetSolutionStepValue(FLUID_FRACTION_GRADIENT) = FluidFractionGradient;
1604 this->EvaluateInPoint(FluidFractionRate,FLUID_FRACTION_RATE,rShapeFunc);
1606 const double EpsilonInside =
false;
1610 rGradEpsOverEps = 1.0 / FluidFraction * FluidFractionGradient ;
1614 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
1616 for (
unsigned int j = 0;
j < TNumNodes; ++
j)
1621 K = Density * rShapeFunc[
i] * AGradN[
j];
1624 K += TauOne * Density * AGradN[
i] * Density * AGradN[
j];
1632 for (
unsigned int m = 0;
m < TDim; ++
m)
1639 G = TauOne * Density * AGradN[
i] * rShapeDeriv(
j,
m);
1642 GAlpha = TauOne * Density * AGradN[
i] * (FluidFraction * rShapeDeriv(
j,
m));
1644 PAlphaDivV = (FluidFraction * rShapeDeriv(
i,
m) + FluidFractionGradient[
m] * rShapeFunc[
i]) * rShapeFunc[
j];
1646 PDivV = rShapeDeriv(
i,
m) * rShapeFunc[
j];
1648 rDampingMatrix(FirstRow +
m, FirstCol + TDim) += Weight * (G - PDivV);
1652 rDampingMatrix(FirstCol + TDim, FirstRow +
m) += Weight * (GAlpha + PAlphaDivV);
1656 L += FluidFraction * rShapeDeriv(
i,
m) * rShapeDeriv(
j,
m);
1658 for (
unsigned int n = 0;
n < TDim; ++
n)
1663 rDampingMatrix(FirstRow +
m, FirstCol +
n) += Weight * TauTwo * rShapeDeriv(
i,
m) * (FluidFraction * rShapeDeriv(
j,
n) + FluidFractionGradient[
n] * rShapeFunc[
j]);
1670 for (
unsigned int d = 0;
d < TDim; ++
d)
1671 rDampingMatrix(FirstRow +
d, FirstCol +
d) +=
K;
1674 rDampingMatrix(FirstRow + TDim, FirstCol + TDim) += Weight * TauOne *
L;
1678 FirstCol += BlockSize;
1684 for (
unsigned int d = 0;
d < TDim; ++
d)
1689 rDampRHS[FirstRow +
d] += Weight * (TauOne * Density * AGradN[
i] * BodyForce[
d] - TauTwo * rShapeDeriv(
i,
d) * FluidFractionRate);
1692 qF += (FluidFraction * rShapeDeriv(
i,
d)) * BodyForce[
d];
1695 rDampRHS[FirstRow + TDim] += Weight * TauOne * qF;
1698 FirstRow += BlockSize;
1703 this->AddViscousTerm(rDampingMatrix,rShapeDeriv,Viscosity*Density*Weight);
1820 const double Density,
1822 double& rElementalMassRes,
1825 const double Weight)
1829 this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv);
1831 double FluidFraction, FluidFractionRate;
1833 this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
1834 this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
1835 this->EvaluateInPoint(FluidFractionRate,FLUID_FRACTION_RATE,rShapeFunc);
1838 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
1844 const double& rPressure = this->GetGeometry()[
i].FastGetSolutionStepValue(PRESSURE);
1847 for (
unsigned int d = 0;
d < TDim; ++
d)
1850 rElementalMomRes[
d] += Weight * (Density * (rShapeFunc[
i] * rBodyForce[
d] - AGradN[
i] * rVelocity[
d]) - rShapeDeriv(
i,
d) * rPressure);
1853 rElementalMassRes -= Weight * (FluidFraction * rShapeDeriv(
i,
d) * rVelocity[
d] + FluidFractionGradient[
d] * rShapeFunc[
i] * rVelocity[
d]);
1858 rElementalMassRes -= Weight * FluidFractionRate;
1873 const double Density,
1877 const double Weight)
1881 this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv);
1883 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
1888 const array_1d< double, 3 > & rAcceleration = this->GetGeometry()[
i].FastGetSolutionStepValue(ACCELERATION);
1890 const double& rPressure = this->GetGeometry()[
i].FastGetSolutionStepValue(PRESSURE);
1893 for (
unsigned int d = 0;
d < TDim; ++
d)
1895 rElementalMomRes[
d] += Weight * (Density * (rShapeFunc[
i] * (rBodyForce[
d] - rAcceleration[
d]) - AGradN[
i] * rVelocity[
d]) - rShapeDeriv(
i,
d) * rPressure);
1911 const double Density,
1915 const double Weight)
1919 this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv);
1921 for (
unsigned int i = 0;
i < TNumNodes; ++
i)
1928 const double& rPressure = this->GetGeometry()[
i].FastGetSolutionStepValue(PRESSURE);
1931 for (
unsigned int d = 0;
d < TDim; ++
d)
1933 rElementalMomRes[
d] += Weight * (Density * (rShapeFunc[
i] * rBodyForce[
d] - AGradN[
i] * rVelocity[
d]) - rShapeDeriv(
i,
d) * rPressure);
1934 rElementalMomRes[
d] -= Weight * rShapeFunc[
i] * rProjection[
d];
1951 const double MolecularViscosity,
1954 double& TotalViscosity,
1957 const double C = this->
GetValue(C_SMAGORINSKY);
1959 TotalViscosity = MolecularViscosity;
1963 const double FilterWidth = this->FilterWidth(rShapeDeriv);
1965 const double NormS = this->SymmetricGradientNorm(rShapeDeriv);
1968 TotalViscosity += 2.0 *
C *
C * FilterWidth * NormS;
1983 rAdvVel = rShapeFunc[0] * (this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY));
1984 for (
unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
1985 rAdvVel += rShapeFunc[iNode] * (this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY));
1998 const std::size_t Step)
2001 rAdvVel = rShapeFunc[0] * (this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY, Step));
2002 for (
unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
2003 rAdvVel += rShapeFunc[iNode] * (this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY, Step));
2020 for (
unsigned int iNode = 1; iNode < TNumNodes; ++iNode){
2021 const array_1d< double, 3 > vel_at_nodes = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY);
2023 for (
unsigned int d = 1;
d < TDim; ++
d){
2025 rAdvVelDiv += vel_at_nodes[
d] * rShapeDeriv(iNode,
d);
2034 const std::size_t Step)
2038 for (
unsigned int iNode = 1; iNode < TNumNodes; ++iNode){
2039 array_1d< double, 3 > vel_at_nodes = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY, Step);
2041 for (
unsigned int d = 1;
d < TDim; ++
d){
2043 rAdvVelDiv += vel_at_nodes[
d] * rShapeDeriv(iNode,
d);
2064 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
2067 rResult[iNode] = rVelocity[0] * rShapeDeriv(iNode, 0);
2068 for (
unsigned int d = 1;
d < TDim; ++
d)
2069 rResult[iNode] += rVelocity[
d] * rShapeDeriv(iNode,
d);
2096 const double & rVelocityDiv,
2101 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
2103 rResult[iNode] = rVelocityDiv * rShapeFunc[iNode];
2105 for (
unsigned int d = 0;
d < TDim; ++
d){
2106 rResult[iNode] += rVelocity[
d] * rShapeDeriv(iNode,
d);
2115 const double Weight = 1.0)
2118 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
2119 rResult += Weight * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2137 rResult = rShapeFunc[0] * this->GetGeometry()[0].FastGetSolutionStepValue(rVariable);
2139 for (
unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
2140 rResult += rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2150 for (
unsigned int i = 0;
i < TNumNodes; ++
i) {
2151 double& scalar = this->GetGeometry()[
i].FastGetSolutionStepValue(rVariable);
2153 for (
unsigned int d = 0;
d < TDim; ++
d){
2154 rResult[
d] += rShapeDeriv(
i,
d) * scalar;
2163 for (
unsigned int d1 = 0;
d1 < TDim; ++
d1) {
2165 for (
unsigned int i = 0;
i < TNumNodes; ++
i) {
2168 for (
unsigned int d2 = 0; d2 < TDim; ++d2) {
2169 rResult(
d1, d2) += rShapeDeriv(
i, d2) * vector[
d1];
2178 const double& DeltaTime,
2179 const std::vector<double>& rSchemeWeigths)
2183 if (rVariable == FLUID_FRACTION_RATE){
2185 double delta_time_inv = 1.0 / DeltaTime;
2199 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
2200 double rate = delta_time_inv * (this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION) - this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION_OLD));
2201 this->GetGeometry()[iNode].SetLock();
2202 this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION_RATE) = rate;
2203 this->GetGeometry()[iNode].UnSetLock();
2204 rResult += rShapeFunc[iNode] * rate;
2226 const double Weight = 1.0)
2229 for (
unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
2230 rResult += Weight * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2247 rResult = rShapeFunc[0] * this->GetGeometry()[0].FastGetSolutionStepValue(rVariable);
2248 for (
unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
2249 rResult += rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2275 const unsigned int GradientSize = (TDim*(TDim+1))/2;
2280 for (
unsigned int k = 0;
k < TNumNodes; ++
k)
2284 for (
unsigned int i = 0;
i < TDim; ++
i)
2286 for (
unsigned int j = 0;
j <
i; ++
j)
2287 GradientVector[
Index++] += 0.5 * (rShapeDeriv(
k,
j) * rNodeVel[
i] + rShapeDeriv(
k,
i) * rNodeVel[
j]);
2288 GradientVector[
Index++] += rShapeDeriv(
k,
i) * rNodeVel[
i];
2295 for (
unsigned int i = 0;
i < TDim; ++
i)
2297 for (
unsigned int j = 0;
j <
i; ++
j)
2299 NormS += 2.0 * GradientVector[
Index] * GradientVector[
Index];
2302 NormS += GradientVector[
Index] * GradientVector[
Index];
2306 NormS = sqrt( 2.0 * NormS );
2319 const double Weight);
2334 const double Weight)
2339 this->CalculateC(
C,Weight);
2341 const unsigned int BlockSize = TDim + 1;
2342 const unsigned int StrainSize = (TDim*TNumNodes)/2;
2345 for(
unsigned int i=0;
i<TNumNodes;
i++)
2347 int base_index = TDim*
i;
2348 int aux_index = BlockSize*
i;
2349 for(
unsigned int j=0;
j<TDim;
j++)
2351 aux[base_index+
j] = aux_index+
j;
2355 for(
unsigned int k=0;
k< StrainSize;
k++)
2357 for(
unsigned int l=0;
l< StrainSize;
l++)
2359 const double Ckl =
C(
k,
l);
2360 for (
unsigned int i = 0;
i < TDim*TNumNodes; ++
i)
2362 const double Bki=
B(
k,
i);
2363 for (
unsigned int j = 0;
j < TDim*TNumNodes; ++
j)
2365 rDampingMatrix(aux[
i],aux[
j]) += Bki*Ckl*
B(
l,
j);
2376 const double Weight)
2382 for (
unsigned int n = 0;
n < TNumNodes;
n++)
2384 const array_1d<double,3>& rVel = this->GetGeometry()[
n].FastGetSolutionStepValue(VELOCITY);
2385 for (
unsigned int i = 0;
i < TDim;
i++)
2386 for (
unsigned int j = 0;
j < TDim;
j++)
2387 GradU(
i,
j) += rDN_DX(
n,
j)*rVel[
i];
2392 Delta[0] = fabs(rGeom[TNumNodes-1].
X()-rGeom[0].
X());
2393 Delta[1] = fabs(rGeom[TNumNodes-1].
Y()-rGeom[0].
Y());
2394 Delta[2] = fabs(rGeom[TNumNodes-1].
Z()-rGeom[0].
Z());
2396 for (
unsigned int n = 1;
n < TNumNodes;
n++)
2398 double hx = fabs(rGeom[
n].
X()-rGeom[
n-1].
X());
2399 if (hx > Delta[0]) Delta[0] = hx;
2400 double hy = fabs(rGeom[
n].
Y()-rGeom[
n-1].
Y());
2401 if (hy > Delta[1]) Delta[1] = hy;
2402 double hz = fabs(rGeom[
n].
Z()-rGeom[
n-1].
Z());
2403 if (hz > Delta[2]) Delta[2] = hz;
2406 double AvgDeltaSq = Delta[0];
2407 for (
unsigned int d = 1;
d < TDim;
d++)
2408 AvgDeltaSq *= Delta[
d];
2409 AvgDeltaSq = std::pow(AvgDeltaSq,2./TDim);
2411 Delta[0] = Delta[0]*Delta[0]/12.0;
2412 Delta[1] = Delta[1]*Delta[1]/12.0;
2413 Delta[2] = Delta[2]*Delta[2]/12.0;
2417 for (
unsigned int i = 0;
i < TDim;
i++)
2418 for (
unsigned int j = 0;
j < TDim;
j++)
2419 for (
unsigned int d = 0;
d < TDim;
d++)
2420 G(
i,
j) += Delta[
d]*GradU(
i,
d)*GradU(
j,
d);
2423 double GijSij = 0.0;
2424 for (
unsigned int i = 0;
i < TDim;
i++)
2425 for (
unsigned int j = 0;
j < TDim;
j++)
2426 GijSij += 0.5*G(
i,
j)*( GradU(
i,
j) + GradU(
j,
i) );
2431 double Gkk = G(0,0);
2432 for (
unsigned int d = 1;
d < TDim;
d++)
2436 const double Ce = 1.0;
2439 double ksgs = -4*AvgDeltaSq*GijSij/(Ce*Ce*Gkk);
2442 unsigned int RowIndex = 0;
2443 unsigned int ColIndex = 0;
2445 for (
unsigned int i = 0;
i < TNumNodes;
i++)
2447 for (
unsigned int j = 0;
j < TNumNodes;
j++)
2449 for (
unsigned int d = 0;
d < TDim;
d++)
2451 double Aux = rDN_DX(
i,
d) * Delta[0] * G(
d,0)*rDN_DX(
j,0);
2452 for (
unsigned int k = 1;
k < TDim;
k++)
2453 Aux += rDN_DX(
i,
d) *Delta[
k] * G(
d,
k)*rDN_DX(
j,
k);
2454 rDampingMatrix(RowIndex+
d,ColIndex+
d) += Weight * 2.0*ksgs * Aux;
2484 const double Viscosity);
2520 virtual void save(
Serializer& rSerializer)
const override
2555 MonolithicDEMCoupled &
operator=(MonolithicDEMCoupled
const& rOther);
2558 MonolithicDEMCoupled(MonolithicDEMCoupled
const& rOther);
2576 template<
unsigned int TDim,
2577 unsigned int TNumNodes >
2585 template<
unsigned int TDim,
2586 unsigned int TNumNodes >
2591 rOStream << std::endl;
PeriodicInterfaceProcess & operator=(const PeriodicInterfaceProcess &)=delete
Base class for all Elements.
Definition: element.h:60
std::vector< DofType::Pointer > DofsVectorType
Definition: element.h:100
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: element.h:1135
virtual int Check(const ProcessInfo &rCurrentProcessInfo) const
Definition: element.h:904
std::vector< std::size_t > EquationIdVectorType
Definition: element.h:98
std::size_t IndexType
Definition: flags.h:74
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometrical_object.h:248
Geometry base class.
Definition: geometry.h:71
static void CalculateGeometryData(const GeometryType &rGeometry, BoundedMatrix< double, 4, 3 > &rDN_DX, array_1d< double, 4 > &rN, double &rVolume)
This function is designed to compute the shape function derivatives, shape functions and volume in 3D...
Definition: geometry_utilities.h:176
This object defines an indexed object.
Definition: indexed_object.h:54
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
A stabilized element for the incompressible Navier-Stokes equations.
Definition: monolithic_dem_coupled.h:137
void AddBTransCB(MatrixType &rDampingMatrix, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Adds the contribution of the viscous term to the momentum equation (alternate).
Definition: monolithic_dem_coupled.h:2332
virtual void AddPointContribution(double &rResult, const Variable< double > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight=1.0)
Definition: monolithic_dem_coupled.h:2112
virtual void EvaluateGradientOfVectorInPoint(BoundedMatrix< double, TDim, TDim > &rResult, const Variable< array_1d< double, 3 > > &rVariable, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Definition: monolithic_dem_coupled.h:2159
void GetConvectionOperator(array_1d< double, TNumNodes > &rResult, const array_1d< double, 3 > &rVelocity, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Write the convective operator evaluated at this point (for each nodal funciton) to an array.
Definition: monolithic_dem_coupled.h:2059
IndexedObject BaseType
base type: an IndexedObject that automatically has a unique number
Definition: monolithic_dem_coupled.h:146
virtual void CalculateLeftHandSide(MatrixType &rLeftHandSideMatrix, const ProcessInfo &rCurrentProcessInfo) override
Returns a zero matrix of appropiate size (provided for compatibility with scheme)
Definition: monolithic_dem_coupled.h:304
virtual void GetEffectiveViscosity(const double Density, const double MolecularViscosity, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, double &TotalViscosity, const ProcessInfo &rCurrentProcessInfo)
Add the contribution from Smagorinsky model to the (kinematic) viscosity.
Definition: monolithic_dem_coupled.h:1950
void CalculateB(BoundedMatrix< double,(TDim *TNumNodes)/2, TDim *TNumNodes > &rB, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Calculate the strain rate matrix.
virtual void CalculateMassMatrix(MatrixType &rMassMatrix, const ProcessInfo &rCurrentProcessInfo) override
Computes local contributions to the mass matrix.
Definition: monolithic_dem_coupled.h:467
virtual void FinalizeNonLinearIteration(const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:662
virtual void CalculateOnIntegrationPoints(const Variable< array_1d< double, 6 > > &rVariable, std::vector< array_1d< double, 6 > > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:1116
virtual void CalculateLocalVelocityContribution(MatrixType &rDampingMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Computes the local contribution associated to 'new' velocity and pressure values.
Definition: monolithic_dem_coupled.h:576
virtual ~MonolithicDEMCoupled()
Destructor.
Definition: monolithic_dem_coupled.h:230
void FinalizeSolutionStep(const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:707
Geometry< NodeType >::PointsArrayType NodesArrayType
definition of nodes container type, redefined from GeometryType
Definition: monolithic_dem_coupled.h:161
virtual void AddRHSLaplacian(VectorType &F, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Definition: monolithic_dem_coupled.h:1336
MonolithicDEMCoupled(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
Constuctor using geometry and properties.
Definition: monolithic_dem_coupled.h:225
virtual void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const override
Provides the global indices for each one of this element's local rows.
virtual void AddViscousTerm(MatrixType &rDampingMatrix, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Adds the contribution of the viscous term to the momentum equation.
virtual void AddMomentumRHS(VectorType &F, const double Density, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight)
Add the momentum equation contribution to the RHS (body forces)
Definition: monolithic_dem_coupled.h:1313
virtual void EvaluateInPoint(double &rResult, const Variable< double > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc)
Write the value of a variable at a point inside the element to a double.
Definition: monolithic_dem_coupled.h:2132
virtual void AddMassRHS(VectorType &F, const double Density, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight, const std::vector< double > &TimeSchemeWeights, const double &DeltaTime)
Definition: monolithic_dem_coupled.h:1357
MonolithicDEMCoupled(IndexType NewId, GeometryType::Pointer pGeometry)
Constructor using a geometry object.
Definition: monolithic_dem_coupled.h:215
virtual void EvaluateGradientOfScalarInPoint(array_1d< double, 3 > &rResult, const Variable< double > &rVariable, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Definition: monolithic_dem_coupled.h:2145
std::vector< Dof< double >::Pointer > DofsVectorType
Definition: monolithic_dem_coupled.h:173
void CalculateWeights(ShapeFunctionDerivativesArrayType &rDN_DX, Matrix &rNContainer, Vector &rGaussWeights)
Determine integration point weights and shape funcition derivatives from the element's geometry.
virtual void AddProjectionToRHS(VectorType &RHS, const array_1d< double, 3 > &rAdvVel, const double Density, const double TauOne, const double TauTwo, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight, const double DeltaTime=1.0)
Add OSS projection terms to the RHS.
Definition: monolithic_dem_coupled.h:1379
virtual void Calculate(const Variable< double > &rVariable, double &rOutput, const ProcessInfo &rCurrentProcessInfo) override
Implementation of Calculate to compute an error estimate.
Definition: monolithic_dem_coupled.h:729
void AddMassStabTerms(MatrixType &rLHSMatrix, const double Density, const array_1d< double, 3 > &rAdvVel, const double TauOne, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Add mass-like stabilization terms to LHS.
Definition: monolithic_dem_coupled.h:1510
void AddConsistentMassMatrixContribution(MatrixType &rLHSMatrix, const array_1d< double, TNumNodes > &rShapeFunc, const double Density, const double Weight)
Definition: monolithic_dem_coupled.h:1463
double SymmetricGradientNorm(const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Compute the norm of the symmetric gradient of velocity.
Definition: monolithic_dem_coupled.h:2273
virtual void CalculateOnIntegrationPoints(const Variable< array_1d< double, 3 > > &rVariable, std::vector< array_1d< double, 3 > > &rOutput, const ProcessInfo &rCurrentProcessInfo) override
Obtain an array_1d<double,3> elemental variable, evaluated on gauss points.
void GetModifiedConvectionOperator(array_1d< double, TNumNodes > &rResult, const array_1d< double, 3 > &rVelocity, const double &rVelocityDiv, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Add the weighted value of a variable at a point inside the element to a double.
Definition: monolithic_dem_coupled.h:2094
virtual void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: monolithic_dem_coupled.h:1205
virtual void GetSecondDerivativesVector(Vector &Values, int Step=0) const override
Returns ACCELERATION_X, ACCELERATION_Y, (ACCELERATION_Z,) 0 for each node.
virtual void GetAdvectiveVel(array_1d< double, 3 > &rAdvVel, const array_1d< double, TNumNodes > &rShapeFunc)
Write the advective velocity evaluated at this point to an array.
Definition: monolithic_dem_coupled.h:1979
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: monolithic_dem_coupled.h:158
virtual void AddPointContribution(array_1d< double, 3 > &rResult, const Variable< array_1d< double, 3 > > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight=1.0)
Add the weighted value of a variable at a point inside the element to a vector.
Definition: monolithic_dem_coupled.h:2223
void ModulatedGradientDiffusion(MatrixType &rDampingMatrix, const BoundedMatrix< double, TNumNodes, TDim > &rDN_DX, const double Weight)
Definition: monolithic_dem_coupled.h:2374
Vector VectorType
Definition: monolithic_dem_coupled.h:163
virtual std::string Info() const override
Turn back information as a string.
Definition: monolithic_dem_coupled.h:1197
std::vector< std::size_t > EquationIdVectorType
Definition: monolithic_dem_coupled.h:171
GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType
Type for an array of shape function gradient matrices.
Definition: monolithic_dem_coupled.h:185
void ASGSMomResidual(const array_1d< double, 3 > &rAdvVel, const double Density, array_1d< double, 3 > &rElementalMomRes, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Assemble the contribution from an integration point to the element's residual.
Definition: monolithic_dem_coupled.h:1872
void OSSMomResidual(const array_1d< double, 3 > &rAdvVel, const double Density, array_1d< double, 3 > &rElementalMomRes, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Assemble the contribution from an integration point to the element's residual.
Definition: monolithic_dem_coupled.h:1910
virtual void GetFirstDerivativesVector(Vector &Values, int Step=0) const override
Returns VELOCITY_X, VELOCITY_Y, (VELOCITY_Z,) PRESSURE for each node.
Element::Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const override
Create a new element of this type.
Definition: monolithic_dem_coupled.h:251
virtual void GetAdvectiveVelDivergence(double &rAdvVelDiv, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const std::size_t Step)
Definition: monolithic_dem_coupled.h:2032
double FilterWidth()
Return the filter width for the smagorinsky model (Delta squared)
virtual void GetDofList(DofsVectorType &rElementalDofList, const ProcessInfo &rCurrentProcessInfo) const override
Returns a list of the element's Dofs.
virtual void CalculateTau(double &TauOne, double &TauTwo, const array_1d< double, 3 > &rAdvVel, const double Area, const double Density, const double KinViscosity, const ProcessInfo &rCurrentProcessInfo)
Calculate Stabilization parameters.
Definition: monolithic_dem_coupled.h:1256
Matrix MatrixType
Definition: monolithic_dem_coupled.h:165
virtual void CalculateRightHandSide(VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Provides local contributions from body forces and projections to the RHS.
Definition: monolithic_dem_coupled.h:335
MonolithicDEMCoupled(IndexType NewId, const NodesArrayType &ThisNodes)
Constructor using an array of nodes.
Definition: monolithic_dem_coupled.h:206
virtual void Calculate(const Variable< array_1d< double, 3 > > &rVariable, array_1d< double, 3 > &rOutput, const ProcessInfo &rCurrentProcessInfo) override
Implementation of Calculate to compute the local OSS projections.
Definition: monolithic_dem_coupled.h:835
KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(MonolithicDEMCoupled)
Pointer definition of MonolithicDEMCoupled.
double ConsistentMassCoef(const double Area)
std::size_t SizeType
Definition: monolithic_dem_coupled.h:169
virtual int Check(const ProcessInfo &rCurrentProcessInfo) const override
Checks the input and that all required Kratos variables have been registered.
Definition: monolithic_dem_coupled.h:1148
PointerVectorSet< Dof< double >, IndexedObject > DofsArrayType
Definition: monolithic_dem_coupled.h:175
void AddProjectionResidualContribution(const array_1d< double, 3 > &rAdvVel, const double Density, array_1d< double, 3 > &rElementalMomRes, double &rElementalMassRes, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Assemble the contribution from an integration point to the element's residual.
Definition: monolithic_dem_coupled.h:1819
MonolithicDEMCoupled(IndexType NewId=0)
Default constuctor.
Definition: monolithic_dem_coupled.h:197
void CalculateLumpedMassMatrix(MatrixType &rLHSMatrix, const double Mass)
Add lumped mass matrix.
Definition: monolithic_dem_coupled.h:1432
void CalculateLaplacianLumpedMassMatrix(MatrixType &rLHSMatrix, const double Mass)
Definition: monolithic_dem_coupled.h:1448
virtual void EvaluateInPoint(array_1d< double, 3 > &rResult, const Variable< array_1d< double, 3 > > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc)
Write the value of a variable at a point inside the element to a double.
Definition: monolithic_dem_coupled.h:2242
Properties PropertiesType
Definition: monolithic_dem_coupled.h:155
Kratos::Matrix ShapeFunctionDerivativesType
Type for a matrix containing the shape function gradients.
Definition: monolithic_dem_coupled.h:182
virtual void CalculateOnIntegrationPoints(const Variable< Matrix > &rVariable, std::vector< Matrix > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:1126
virtual void CalculateOnIntegrationPoints(const Variable< double > &rVariable, std::vector< double > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Obtain a double elemental variable, evaluated on gauss points.
Definition: monolithic_dem_coupled.h:1001
Node NodeType
definition of node type (default is: Node)
Definition: monolithic_dem_coupled.h:149
virtual void CalculateOnIntegrationPoints(const Variable< Vector > &rVariable, std::vector< Vector > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:1121
virtual void GetAdvectiveVel(array_1d< double, 3 > &rAdvVel, const array_1d< double, TNumNodes > &rShapeFunc, const std::size_t Step)
Write the advective velocity evaluated at this point to an array.
Definition: monolithic_dem_coupled.h:1996
std::size_t IndexType
Definition: monolithic_dem_coupled.h:167
virtual void EvaluateTimeDerivativeInPoint(double &rResult, const Variable< double > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc, const double &DeltaTime, const std::vector< double > &rSchemeWeigths)
Definition: monolithic_dem_coupled.h:2175
double FilterWidth(const BoundedMatrix< double, TNumNodes, TDim > &DN_DX)
Return the filter width for the smagorinsky model (Delta squared)
Kratos::Vector ShapeFunctionsType
Type for shape function values container.
Definition: monolithic_dem_coupled.h:179
virtual void CalculateC(BoundedMatrix< double,(TDim *TNumNodes)/2,(TDim *TNumNodes)/2 > &rC, const double Viscosity)
Calculate a matrix that provides the stress given the strain rate.
virtual void CalculateStaticTau(double &TauOne, const array_1d< double, 3 > &rAdvVel, const double Area, const double Density, const double KinViscosity)
Calculate momentum stabilization parameter (without time term).
Definition: monolithic_dem_coupled.h:1294
virtual void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Provides local contributions from body forces and OSS projection terms.
Definition: monolithic_dem_coupled.h:269
double ElementSize(const double)
Return an estimate for the element size h, used to calculate the stabilization parameters.
virtual void CalculateLaplacianMassMatrix(MatrixType &rMassMatrix, const ProcessInfo &rCurrentProcessInfo)
Definition: monolithic_dem_coupled.h:542
void AddIntegrationPointVelocityContribution(MatrixType &rDampingMatrix, VectorType &rDampRHS, const double Density, const double Viscosity, const array_1d< double, 3 > &rAdvVel, const double TauOne, const double TauTwo, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Add a the contribution from a single integration point to the velocity contribution.
Definition: monolithic_dem_coupled.h:1566
virtual void GetAdvectiveVelDivergence(double &rAdvVelDiv, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Write the divergence of the advective velocity evaluated at this point to an array.
Definition: monolithic_dem_coupled.h:2015
This class defines the node.
Definition: node.h:65
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
#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_CHECK_DOF_IN_NODE(TheVariable, TheNode)
Definition: checks.h:176
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
#define KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(TheVariable, TheNode)
Definition: checks.h:171
Modeler::Pointer Create(const std::string &ModelerName, Model &rModel, const Parameters ModelParameters)
Checks if the modeler is registered.
Definition: modeler_factory.cpp:30
Parameters GetValue(Parameters &rParameters, const std::string &rEntry)
Definition: add_kratos_parameters_to_python.cpp:53
void CalculateB(const GeometricalObject &rElement, const TMatrixType1 &rDN_DX, TMatrixType2 &rB)
This method computes the deformation tensor B (for small deformation solid elements)
Definition: structural_mechanics_element_utilities.h:105
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
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
def SetValue(entity, variable, value)
Definition: coupling_interface_data.py:256
X2
Definition: generate_frictional_mortar_condition.py:120
X1
Definition: generate_frictional_mortar_condition.py:119
int C
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:27
F
Definition: hinsberg_optimization.py:144
def load(f)
Definition: ode_solve.py:307
int L
Definition: ode_solve.py:390
int d
Definition: ode_solve.py:397
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
int m
Definition: run_marine_rain_substepping.py:8
J
Definition: sensitivityMatrix.py:58
N
Definition: sensitivityMatrix.py:29
K
Definition: sensitivityMatrix.py:73
B
Definition: sensitivityMatrix.py:76
integer i
Definition: TensorModule.f:17
subroutine d1(DSTRAN, D, dtime, NDI, NSHR, NTENS)
Definition: TensorModule.f:594
integer l
Definition: TensorModule.f:17