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.
dam_UP_scheme.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosDamApplication $
3 // Last Modified by: $Author: Lorenzo Gracia $
4 // Date: $Date: January 2016 $
5 // Revision: $Revision: 1.0 $
6 //
7 
8 #if !defined(KRATOS_DAM_UP_SCHEME )
9 #define KRATOS_DAM_UP_SCHEME
10 
11 // Project includes
12 #include "includes/define.h"
13 #include "includes/model_part.h"
14 #include "solving_strategies/schemes/scheme.h"
15 #include "includes/variables.h"
16 #include "containers/array_1d.h"
17 #include "includes/element.h"
18 
19 // Application includes
20 #include "dam_application_variables.h"
21 
22 namespace Kratos
23 {
24 
25 template<class TSparseSpace, class TDenseSpace>
26 
27 class DamUPScheme : public Scheme<TSparseSpace,TDenseSpace>
28 {
29 
30 public:
31 
32  KRATOS_CLASS_POINTER_DEFINITION( DamUPScheme );
33 
34  typedef Scheme<TSparseSpace,TDenseSpace> BaseType;
35  typedef typename BaseType::DofsArrayType DofsArrayType;
36  typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
37  typedef typename BaseType::TSystemVectorType TSystemVectorType;
38  typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
39  typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
40 
41 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
42 
44  DamUPScheme(double beta, double gamma, double rayleigh_m ,double rayleigh_k ): Scheme<TSparseSpace,TDenseSpace>()
45  {
46  mBeta = beta;
47  mGamma = gamma;
48  mrayleigh_m = rayleigh_m;
49  mrayleigh_k = rayleigh_k;
50 
51 
52  //Allocate auxiliary memory
53  int NumThreads = ParallelUtilities::GetNumThreads();
54  mMassMatrix.resize(NumThreads);
55  mAccelerationVector.resize(NumThreads);
56  mDampingMatrix.resize(NumThreads);
57  mVelocityVector.resize(NumThreads);
58  }
59 
60  //------------------------------------------------------------------------------------
61 
63  virtual ~DamUPScheme() {}
64 
65 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
66 
67  int Check(ModelPart& r_model_part) override
68  {
69  KRATOS_TRY
70 
71  //check for variables keys (verify that the variables are correctly initialized)
72  if(DISPLACEMENT.Key() == 0)
73  KRATOS_THROW_ERROR( std::invalid_argument,"DISPLACEMENT has Key zero! (check if the application is correctly registered", "" )
74  if(VELOCITY.Key() == 0)
75  KRATOS_THROW_ERROR( std::invalid_argument,"VELOCITY has Key zero! (check if the application is correctly registered", "" )
76  if(ACCELERATION.Key() == 0)
77  KRATOS_THROW_ERROR( std::invalid_argument,"ACCELERATION has Key zero! (check if the application is correctly registered", "" )
78  if(PRESSURE.Key() == 0)
79  KRATOS_THROW_ERROR( std::invalid_argument, "PRESSURE has Key zero! (check if the application is correctly registered", "" )
80  if(Dt_PRESSURE.Key() == 0)
81  KRATOS_THROW_ERROR( std::invalid_argument, "Dt_PRESSURE has Key zero! (check if the application is correctly registered", "" )
82  if(Dt2_PRESSURE.Key() == 0)
83  KRATOS_THROW_ERROR( std::invalid_argument, "Dt2_PRESSURE has Key zero! (check if the application is correctly registered", "" )
84  if ( VELOCITY_PRESSURE_COEFFICIENT.Key() == 0 )
85  KRATOS_THROW_ERROR( std::invalid_argument, "VELOCITY_PRESSURE_COEFFICIENT has Key zero! (check if the application is correctly registered", "" )
86  if ( ACCELERATION_PRESSURE_COEFFICIENT.Key() == 0 )
87  KRATOS_THROW_ERROR( std::invalid_argument, "ACCELERATION_PRESSURE_COEFFICIENT has Key zero! (check if the application is correctly registered", "" )
88 
89  //check that variables are correctly allocated
90  for(ModelPart::NodesContainerType::iterator it=r_model_part.NodesBegin(); it!=r_model_part.NodesEnd(); it++)
91  {
92  if(it->SolutionStepsDataHas(DISPLACEMENT) == false)
93  KRATOS_THROW_ERROR( std::logic_error, "DISPLACEMENT variable is not allocated for node ", it->Id() )
94  if(it->SolutionStepsDataHas(VELOCITY) == false)
95  KRATOS_THROW_ERROR( std::logic_error, "VELOCITY variable is not allocated for node ", it->Id() )
96  if(it->SolutionStepsDataHas(ACCELERATION) == false)
97  KRATOS_THROW_ERROR( std::logic_error, "ACCELERATION variable is not allocated for node ", it->Id() )
98  if(it->SolutionStepsDataHas(PRESSURE) == false)
99  KRATOS_THROW_ERROR( std::logic_error, "PRESSURE variable is not allocated for node ", it->Id() )
100  if(it->SolutionStepsDataHas(Dt_PRESSURE) == false)
101  KRATOS_THROW_ERROR( std::logic_error, "Dt_PRESSURE variable is not allocated for node ", it->Id() )
102  if(it->SolutionStepsDataHas(Dt2_PRESSURE) == false)
103  KRATOS_THROW_ERROR( std::logic_error, "Dt2_PRESSURE variable is not allocated for node ", it->Id() )
104 
105 
106  if(it->HasDofFor(DISPLACEMENT_X) == false)
107  KRATOS_THROW_ERROR( std::invalid_argument,"missing DISPLACEMENT_X dof on node ",it->Id() )
108  if(it->HasDofFor(DISPLACEMENT_Y) == false)
109  KRATOS_THROW_ERROR( std::invalid_argument,"missing DISPLACEMENT_Y dof on node ",it->Id() )
110  if(it->HasDofFor(DISPLACEMENT_Z) == false)
111  KRATOS_THROW_ERROR( std::invalid_argument,"missing DISPLACEMENT_Z dof on node ",it->Id() )
112  if(it->HasDofFor(PRESSURE) == false)
113  KRATOS_THROW_ERROR( std::invalid_argument,"missing PRESSURE dof on node ",it->Id() )
114  }
115 
116  //check for minimum value of the buffer index.
117  if (r_model_part.GetBufferSize() < 2)
118  KRATOS_THROW_ERROR( std::logic_error, "insufficient buffer size. Buffer size should be greater than 2. Current size is", r_model_part.GetBufferSize() )
119 
120  // Check beta, gamma and theta
121  if(mBeta <= 0.0 || mGamma<= 0.0)
122  KRATOS_THROW_ERROR( std::invalid_argument,"Some of the scheme variables: beta or gamma has an invalid value ", "" )
123 
124  return 0;
125 
126  KRATOS_CATCH( "" )
127  }
128 
129 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
130 
131  void Initialize(ModelPart& r_model_part) override
132  {
133  KRATOS_TRY
134 
135  mDeltaTime = r_model_part.GetProcessInfo()[DELTA_TIME];
136  r_model_part.GetProcessInfo()[VELOCITY_PRESSURE_COEFFICIENT] = mGamma/(mBeta*mDeltaTime);
137  r_model_part.GetProcessInfo()[ACCELERATION_PRESSURE_COEFFICIENT] = 1.0/(mBeta*mDeltaTime*mDeltaTime);
138 
139  r_model_part.GetProcessInfo()[RAYLEIGH_ALPHA] = mrayleigh_m;
140  r_model_part.GetProcessInfo()[RAYLEIGH_BETA] = mrayleigh_k;
141 
142  KRATOS_WATCH(mrayleigh_m)
143  KRATOS_WATCH(mrayleigh_k)
144 
145  BaseType::mSchemeIsInitialized = true;
146 
147  KRATOS_CATCH("")
148  }
149 
150 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
151 
152  void InitializeSolutionStep(
153  ModelPart& r_model_part,
154  TSystemMatrixType& A,
155  TSystemVectorType& Dx,
156  TSystemVectorType& b) override
157  {
158  KRATOS_TRY
159 
160  mDeltaTime = r_model_part.GetProcessInfo()[DELTA_TIME];
161  r_model_part.GetProcessInfo()[VELOCITY_PRESSURE_COEFFICIENT] = mGamma/(mBeta*mDeltaTime);
162  r_model_part.GetProcessInfo()[ACCELERATION_PRESSURE_COEFFICIENT] = 1.0/(mBeta*mDeltaTime*mDeltaTime);
163 
164  const ProcessInfo& CurrentProcessInfo = r_model_part.GetProcessInfo();
165 
166  int NElems = static_cast<int>(r_model_part.Elements().size());
167  ModelPart::ElementsContainerType::iterator el_begin = r_model_part.ElementsBegin();
168 
169  #pragma omp parallel for
170  for(int i = 0; i < NElems; i++)
171  {
172  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
173  itElem -> InitializeSolutionStep(CurrentProcessInfo);
174  }
175 
176  int NCons = static_cast<int>(r_model_part.Conditions().size());
177  ModelPart::ConditionsContainerType::iterator con_begin = r_model_part.ConditionsBegin();
178 
179  #pragma omp parallel for
180  for(int i = 0; i < NCons; i++)
181  {
182  ModelPart::ConditionsContainerType::iterator itCond = con_begin + i;
183  itCond -> InitializeSolutionStep(CurrentProcessInfo);
184  }
185 
186  KRATOS_CATCH("")
187  }
188 
189 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
190 
191  void Predict(
192  ModelPart& r_model_part,
193  DofsArrayType& rDofSet,
194  TSystemMatrixType& A,
195  TSystemVectorType& Dx,
196  TSystemVectorType& b) override
197  {
198 
199  KRATOS_TRY
200 
201  // Updating Velocities, Acclerations, DtPressure and Dt2Pressure
202  array_1d<double,3> DeltaDisplacement;
203  double DeltaPressure;
204 
205  const int NNodes = static_cast<int>(r_model_part.Nodes().size());
206  ModelPart::NodesContainerType::iterator node_begin = r_model_part.NodesBegin();
207 
208  #pragma omp parallel for private(DeltaDisplacement,DeltaPressure)
209  for(int i = 0; i < NNodes; i++)
210  {
211  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
212 
213  // Termes related to displacement field
214 
215  array_1d<double,3>& CurrentDisplacement = itNode->FastGetSolutionStepValue(DISPLACEMENT);
216  array_1d<double,3>& CurrentAcceleration = itNode->FastGetSolutionStepValue(ACCELERATION);
217  array_1d<double,3>& CurrentVelocity = itNode->FastGetSolutionStepValue(VELOCITY);
218 
219  const array_1d<double,3>& PreviousDisplacement = itNode->FastGetSolutionStepValue(DISPLACEMENT, 1);
220  const array_1d<double,3>& PreviousAcceleration = itNode->FastGetSolutionStepValue(ACCELERATION, 1);
221  const array_1d<double,3>& PreviousVelocity = itNode->FastGetSolutionStepValue(VELOCITY, 1);
222 
223  if (itNode -> IsFixed(ACCELERATION_X))
224  {
225  CurrentDisplacement[0] = PreviousDisplacement[0] + mDeltaTime * PreviousVelocity[0] + std::pow(mDeltaTime, 2) * ( ( 0.5 - mBeta) * PreviousAcceleration[0] + mBeta * CurrentAcceleration[0] );
226  }
227  else if (itNode -> IsFixed(VELOCITY_X))
228  {
229  CurrentDisplacement[0] = PreviousDisplacement[0] + mDeltaTime*(mBeta/mGamma*(CurrentVelocity[0]-PreviousVelocity[0])+PreviousVelocity[0]);
230  }
231  else if (itNode -> IsFixed(DISPLACEMENT_X) == false)
232  {
233  CurrentDisplacement[0] = PreviousDisplacement[0] + mDeltaTime * PreviousVelocity[0] + 0.5 * std::pow(mDeltaTime, 2) * PreviousAcceleration[0];
234  }
235 
236  if (itNode -> IsFixed(ACCELERATION_Y))
237  {
238  CurrentDisplacement[1] = PreviousDisplacement[1] + mDeltaTime * PreviousVelocity[1] + std::pow(mDeltaTime, 2) * ( ( 0.5 - mBeta) * PreviousAcceleration[1] + mBeta * CurrentAcceleration[1] );
239  }
240  else if (itNode -> IsFixed(VELOCITY_Y))
241  {
242  CurrentDisplacement[1] = PreviousDisplacement[1] + mDeltaTime*(mBeta/mGamma*(CurrentVelocity[1]-PreviousVelocity[1])+PreviousVelocity[1]);
243  }
244  else if (itNode -> IsFixed(DISPLACEMENT_Y) == false)
245  {
246  CurrentDisplacement[1] = PreviousDisplacement[1] + mDeltaTime * PreviousVelocity[1] + 0.5 * std::pow(mDeltaTime, 2) * PreviousAcceleration[1];
247  }
248 
249  // For 3D cases
250  if (itNode -> HasDofFor(DISPLACEMENT_Z))
251  {
252  if (itNode -> IsFixed(ACCELERATION_Z))
253  {
254  CurrentDisplacement[2] = PreviousDisplacement[2] + mDeltaTime * PreviousVelocity[2] + std::pow(mDeltaTime, 2) * ( ( 0.5 - mBeta) * PreviousAcceleration[2] + mBeta * CurrentAcceleration[2] );
255  }
256  else if (itNode -> IsFixed(VELOCITY_Z))
257  {
258  CurrentDisplacement[2] = PreviousDisplacement[2] + mDeltaTime*(mBeta/mGamma*(CurrentVelocity[2]-PreviousVelocity[2])+PreviousVelocity[2]);
259  }
260  else if (itNode -> IsFixed(DISPLACEMENT_Z) == false)
261  {
262  CurrentDisplacement[2] = PreviousDisplacement[2] + mDeltaTime * PreviousVelocity[2] + 0.5 * std::pow(mDeltaTime, 2) * PreviousAcceleration[2];
263  }
264  }
265 
266  noalias(DeltaDisplacement) = CurrentDisplacement - PreviousDisplacement;
267 
268  noalias(CurrentAcceleration) = 1.0/(mBeta*mDeltaTime*mDeltaTime)*(DeltaDisplacement - mDeltaTime*PreviousVelocity - (0.5-mBeta)*mDeltaTime*mDeltaTime*PreviousAcceleration);
269  noalias(CurrentVelocity) = PreviousVelocity + (1.0-mGamma)*mDeltaTime*PreviousAcceleration + mGamma*mDeltaTime*CurrentAcceleration;
270 
271  // Terms related to Pressure field
272 
273  double& CurrentDt2Pressure = itNode->FastGetSolutionStepValue(Dt2_PRESSURE);
274  double& CurrentDtPressure = itNode->FastGetSolutionStepValue(Dt_PRESSURE);
275  DeltaPressure = itNode->FastGetSolutionStepValue(PRESSURE) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
276  const double& PreviousDt2Pressure = itNode->FastGetSolutionStepValue(Dt2_PRESSURE, 1);
277  const double& PreviousDtPressure = itNode->FastGetSolutionStepValue(Dt_PRESSURE, 1);
278 
279  CurrentDt2Pressure = 1.0/(mBeta*mDeltaTime*mDeltaTime)*(DeltaPressure - mDeltaTime*PreviousDtPressure - (0.5-mBeta)*mDeltaTime*mDeltaTime*PreviousDt2Pressure);
280  CurrentDtPressure = PreviousDtPressure + (1.0-mGamma)*mDeltaTime*PreviousDt2Pressure + mGamma*mDeltaTime*CurrentDt2Pressure;
281 
282  }
283 
284  KRATOS_CATCH( "" )
285  }
286 
287 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
288 
289  void InitializeNonLinIteration(
290  ModelPart& r_model_part,
291  TSystemMatrixType& A,
292  TSystemVectorType& Dx,
293  TSystemVectorType& b) override
294  {
295  KRATOS_TRY
296 
297  const ProcessInfo& CurrentProcessInfo = r_model_part.GetProcessInfo();
298 
299  int NElems = static_cast<int>(r_model_part.Elements().size());
300  ModelPart::ElementsContainerType::iterator el_begin = r_model_part.ElementsBegin();
301 
302  #pragma omp parallel for
303  for(int i = 0; i < NElems; i++)
304  {
305  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
306  itElem -> InitializeNonLinearIteration(CurrentProcessInfo);
307  }
308 
309  int NCons = static_cast<int>(r_model_part.Conditions().size());
310  ModelPart::ConditionsContainerType::iterator con_begin = r_model_part.ConditionsBegin();
311 
312  #pragma omp parallel for
313  for(int i = 0; i < NCons; i++)
314  {
315  ModelPart::ConditionsContainerType::iterator itCond = con_begin + i;
316  itCond -> InitializeNonLinearIteration(CurrentProcessInfo);
317  }
318 
319  KRATOS_CATCH("")
320  }
321 
322 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
323 
324  void FinalizeNonLinIteration(
325  ModelPart& r_model_part,
326  TSystemMatrixType& A,
327  TSystemVectorType& Dx,
328  TSystemVectorType& b) override
329  {
330  KRATOS_TRY
331 
332  const ProcessInfo& CurrentProcessInfo = r_model_part.GetProcessInfo();
333 
334  int NElems = static_cast<int>(r_model_part.Elements().size());
335  ModelPart::ElementsContainerType::iterator el_begin = r_model_part.ElementsBegin();
336 
337  #pragma omp parallel for
338  for(int i = 0; i < NElems; i++)
339  {
340  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
341  itElem -> FinalizeNonLinearIteration(CurrentProcessInfo);
342  }
343 
344  int NCons = static_cast<int>(r_model_part.Conditions().size());
345  ModelPart::ConditionsContainerType::iterator con_begin = r_model_part.ConditionsBegin();
346 
347  #pragma omp parallel for
348  for(int i = 0; i < NCons; i++)
349  {
350  ModelPart::ConditionsContainerType::iterator itCond = con_begin + i;
351  itCond -> FinalizeNonLinearIteration(CurrentProcessInfo);
352  }
353 
354  KRATOS_CATCH("")
355  }
356 
357 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
358 
359 // Note: this is in a parallel loop
360 
361  void CalculateSystemContributions(
362  Element& rCurrentElement,
363  LocalSystemMatrixType& LHS_Contribution,
364  LocalSystemVectorType& RHS_Contribution,
365  Element::EquationIdVectorType& EquationId,
366  const ProcessInfo& CurrentProcessInfo) override
367  {
368  KRATOS_TRY
369 
370  int thread = OpenMPUtils::ThisThread();
371 
372  rCurrentElement.CalculateLocalSystem(LHS_Contribution,RHS_Contribution, CurrentProcessInfo);
373 
374  rCurrentElement.EquationIdVector(EquationId, CurrentProcessInfo);
375 
376  rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
377 
378  rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
379 
380  AddDynamicsToLHS (LHS_Contribution, mDampingMatrix[thread], mMassMatrix[thread], CurrentProcessInfo);
381 
382  AddDynamicsToRHS (rCurrentElement, RHS_Contribution, mDampingMatrix[thread], mMassMatrix[thread], CurrentProcessInfo);
383 
384  KRATOS_CATCH( "" )
385  }
386 
387 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
388 
389 // Note: this is in a parallel loop
390 
391  void CalculateSystemContributions(
392  Condition& rCurrentCondition,
393  LocalSystemMatrixType& LHS_Contribution,
394  LocalSystemVectorType& RHS_Contribution,
395  Element::EquationIdVectorType& EquationId,
396  const ProcessInfo& CurrentProcessInfo) override
397  {
398  KRATOS_TRY
399 
400  rCurrentCondition.CalculateLocalSystem(LHS_Contribution, RHS_Contribution, CurrentProcessInfo);
401 
402  rCurrentCondition.EquationIdVector(EquationId, CurrentProcessInfo);
403 
404  KRATOS_CATCH( "" )
405  }
406 
407 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
408 
409 // Note: this is in a parallel loop
410 
411  void Calculate_RHS_Contribution(
412  Element& rCurrentElement,
413  LocalSystemVectorType& RHS_Contribution,
414  Element::EquationIdVectorType& EquationId,
415  const ProcessInfo& CurrentProcessInfo)
416  {
417  KRATOS_TRY
418 
419  int thread = OpenMPUtils::ThisThread();
420 
421  rCurrentElement.CalculateRightHandSide(RHS_Contribution, CurrentProcessInfo);
422 
423  rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
424 
425  rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
426 
427  rCurrentElement.EquationIdVector(EquationId, CurrentProcessInfo);
428 
429  AddDynamicsToRHS (rCurrentElement, RHS_Contribution, mDampingMatrix[thread], mMassMatrix[thread], CurrentProcessInfo);
430 
431  KRATOS_CATCH( "" )
432  }
433 
434 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
435 
436 // Note: this is in a parallel loop
437 
438  void Calculate_RHS_Contribution(
439  Condition& rCurrentCondition,
440  LocalSystemVectorType& RHS_Contribution,
441  Element::EquationIdVectorType& EquationId,
442  const ProcessInfo& CurrentProcessInfo)
443  {
444  KRATOS_TRY
445 
446  rCurrentCondition.CalculateRightHandSide(RHS_Contribution, CurrentProcessInfo);
447 
448  rCurrentCondition.EquationIdVector(EquationId, CurrentProcessInfo);
449 
450 
451  KRATOS_CATCH( "" )
452  }
453 
454 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
455 
456 // Note: this is in a parallel loop
457 
458  void Calculate_LHS_Contribution(
459  Element& rCurrentElement,
460  LocalSystemMatrixType& LHS_Contribution,
461  Element::EquationIdVectorType& EquationId,
462  const ProcessInfo& CurrentProcessInfo)
463  {
464  KRATOS_TRY
465 
466  int thread = OpenMPUtils::ThisThread();
467 
468  rCurrentElement.CalculateLeftHandSide(LHS_Contribution, CurrentProcessInfo);
469 
470  rCurrentElement.EquationIdVector(EquationId,CurrentProcessInfo);
471 
472  rCurrentElement.CalculateMassMatrix(mMassMatrix[thread],CurrentProcessInfo);
473 
474  rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread],CurrentProcessInfo);
475 
476  AddDynamicsToLHS (LHS_Contribution, mDampingMatrix[thread], mMassMatrix[thread], CurrentProcessInfo);
477 
478  KRATOS_CATCH( "" )
479  }
480 
481 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
482 
483 // Note: this is in a parallel loop
484 
485  void Calculate_LHS_Contribution(
486  Condition& rCurrentCondition,
487  LocalSystemMatrixType& LHS_Contribution,
488  Element::EquationIdVectorType& EquationId,
489  const ProcessInfo& CurrentProcessInfo)
490  {
491  KRATOS_TRY
492 
493  rCurrentCondition.CalculateLeftHandSide(LHS_Contribution, CurrentProcessInfo);
494 
495  rCurrentCondition.EquationIdVector(EquationId, CurrentProcessInfo);
496 
497  KRATOS_CATCH( "" )
498  }
499 
500 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
501 
502  void Update(
503  ModelPart& r_model_part,
504  DofsArrayType& rDofSet,
505  TSystemMatrixType& A,
506  TSystemVectorType& Dx,
507  TSystemVectorType& b) override
508  {
509  KRATOS_TRY
510 
511  // Initialize
512  block_for_each(rDofSet, [&Dx](auto& dof)
513  {
514  if (dof.IsFree())
515  {
516  dof.GetSolutionStepValue() += TSparseSpace::GetValue(Dx, dof.EquationId());
517  }
518 
519  }
520  );
521 
522  this->UpdateVariablesDerivatives(r_model_part);
523 
524  KRATOS_CATCH( "" )
525  }
526 
527 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
528 
529 protected:
530 
532 
533  double mBeta;
534  double mGamma;
535  double mDeltaTime;
536  double mrayleigh_m;
537  double mrayleigh_k;
538 
539  std::vector< Matrix > mMassMatrix;
540  std::vector< Vector > mAccelerationVector;
541  std::vector< Matrix > mDampingMatrix;
542  std::vector< Vector > mVelocityVector;
543 
544 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
545 
546  inline void UpdateVariablesDerivatives(ModelPart& r_model_part)
547  {
548  KRATOS_TRY
549 
550  // Updating Velocities, Acclerations, DtPressure and Dt2Pressure
551  array_1d<double,3> DeltaDisplacement;
552  double DeltaPressure;
553 
554  const int NNodes = static_cast<int>(r_model_part.Nodes().size());
555  ModelPart::NodesContainerType::iterator node_begin = r_model_part.NodesBegin();
556 
557  #pragma omp parallel for private(DeltaDisplacement,DeltaPressure)
558  for(int i = 0; i < NNodes; i++)
559  {
560  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
561 
562  // Termes related to displacement field
563 
564  array_1d<double,3>& CurrentDisplacement = itNode->FastGetSolutionStepValue(DISPLACEMENT);
565  array_1d<double,3>& CurrentAcceleration = itNode->FastGetSolutionStepValue(ACCELERATION);
566  array_1d<double,3>& CurrentVelocity = itNode->FastGetSolutionStepValue(VELOCITY);
567 
568  const array_1d<double,3>& PreviousDisplacement = itNode->FastGetSolutionStepValue(DISPLACEMENT, 1);
569  const array_1d<double,3>& PreviousAcceleration = itNode->FastGetSolutionStepValue(ACCELERATION, 1);
570  const array_1d<double,3>& PreviousVelocity = itNode->FastGetSolutionStepValue(VELOCITY, 1);
571 
572  noalias(DeltaDisplacement) = CurrentDisplacement - PreviousDisplacement;
573 
574  noalias(CurrentAcceleration) = 1.0/(mBeta*mDeltaTime*mDeltaTime)*(DeltaDisplacement - mDeltaTime*PreviousVelocity - (0.5-mBeta)*mDeltaTime*mDeltaTime*PreviousAcceleration);
575  noalias(CurrentVelocity) = PreviousVelocity + (1.0-mGamma)*mDeltaTime*PreviousAcceleration + mGamma*mDeltaTime*CurrentAcceleration;
576 
577  // Terms related to Pressure field
578 
579  double& CurrentDt2Pressure = itNode->FastGetSolutionStepValue(Dt2_PRESSURE);
580  double& CurrentDtPressure = itNode->FastGetSolutionStepValue(Dt_PRESSURE);
581  DeltaPressure = itNode->FastGetSolutionStepValue(PRESSURE) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
582  const double& PreviousDt2Pressure = itNode->FastGetSolutionStepValue(Dt2_PRESSURE, 1);
583  const double& PreviousDtPressure = itNode->FastGetSolutionStepValue(Dt_PRESSURE, 1);
584 
585  CurrentDt2Pressure = 1.0/(mBeta*mDeltaTime*mDeltaTime)*(DeltaPressure - mDeltaTime*PreviousDtPressure - (0.5-mBeta)*mDeltaTime*mDeltaTime*PreviousDt2Pressure);
586  CurrentDtPressure = PreviousDtPressure + (1.0-mGamma)*mDeltaTime*PreviousDt2Pressure + mGamma*mDeltaTime*CurrentDt2Pressure;
587 
588  }
589 
590  KRATOS_CATCH( "" )
591  }
592 
593 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
594 
595  void AddDynamicsToLHS(LocalSystemMatrixType& LHS_Contribution, LocalSystemMatrixType& M, LocalSystemMatrixType& C, const ProcessInfo& CurrentProcessInfo)
596  {
597  // adding mass contribution
598  if (M.size1() != 0)
599  noalias(LHS_Contribution) += 1.0/(mBeta*mDeltaTime*mDeltaTime)*M;
600 
601  // adding damping contribution
602  if (C.size1() != 0)
603  noalias(LHS_Contribution) += mGamma/(mBeta*mDeltaTime)*C;
604  }
605 
606 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
607 
608  void AddDynamicsToRHS(Element& rCurrentElement,
609  LocalSystemVectorType& RHS_Contribution,
610  LocalSystemMatrixType& M,
611  LocalSystemMatrixType& C,
612  const ProcessInfo& CurrentProcessInfo)
613  {
614  int thread = OpenMPUtils::ThisThread();
615 
616  //adding inertia contribution
617  if (M.size1() != 0)
618  {
619  rCurrentElement.GetSecondDerivativesVector(mAccelerationVector[thread], 0);
620 
621  noalias(RHS_Contribution) -= prod(M, mAccelerationVector[thread]);
622  }
623 
624  //adding damping contribution
625  if (C.size1() != 0)
626  {
627  rCurrentElement.GetFirstDerivativesVector(mVelocityVector[thread], 0);
628 
629  noalias(RHS_Contribution) -= prod(C, mVelocityVector[thread]);
630  }
631  }
632 
633 
634 }; // Class DamUPScheme
635 } // namespace Kratos
636 
637 #endif // KRATOS_DAM_UP_SCHEME defined
638 
array_1d.h
Kratos::Condition
Base class for all Conditions.
Definition: condition.h:59
Kratos::Condition::CalculateLeftHandSide
virtual void CalculateLeftHandSide(MatrixType &rLeftHandSideMatrix, const ProcessInfo &rCurrentProcessInfo)
Definition: condition.h:426
Kratos::Condition::EquationIdVector
virtual void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const
Definition: condition.h:260
Kratos::Condition::CalculateRightHandSide
virtual void CalculateRightHandSide(VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo)
Definition: condition.h:440
Kratos::Condition::CalculateLocalSystem
virtual void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo)
Definition: condition.h:408
Kratos::DamUPScheme
Definition: dam_UP_scheme.hpp:28
Kratos::DamUPScheme::DofsArrayType
BaseType::DofsArrayType DofsArrayType
Definition: dam_UP_scheme.hpp:35
Kratos::DamUPScheme::Initialize
void Initialize(ModelPart &r_model_part) override
This is the place to initialize the Scheme.
Definition: dam_UP_scheme.hpp:131
Kratos::DamUPScheme::InitializeSolutionStep
void InitializeSolutionStep(ModelPart &r_model_part, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
Function called once at the beginning of each solution step.
Definition: dam_UP_scheme.hpp:152
Kratos::DamUPScheme::Calculate_RHS_Contribution
void Calculate_RHS_Contribution(Condition &rCurrentCondition, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, const ProcessInfo &CurrentProcessInfo)
Definition: dam_UP_scheme.hpp:438
Kratos::DamUPScheme::TSystemVectorType
BaseType::TSystemVectorType TSystemVectorType
Definition: dam_UP_scheme.hpp:37
Kratos::DamUPScheme::Calculate_LHS_Contribution
void Calculate_LHS_Contribution(Element &rCurrentElement, LocalSystemMatrixType &LHS_Contribution, Element::EquationIdVectorType &EquationId, const ProcessInfo &CurrentProcessInfo)
Definition: dam_UP_scheme.hpp:458
Kratos::DamUPScheme::mDampingMatrix
std::vector< Matrix > mDampingMatrix
Definition: dam_UP_scheme.hpp:541
Kratos::DamUPScheme::FinalizeNonLinIteration
void FinalizeNonLinIteration(ModelPart &r_model_part, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
Function to be called when it is needed to finalize an iteration. It is designed to be called at the ...
Definition: dam_UP_scheme.hpp:324
Kratos::DamUPScheme::LocalSystemVectorType
BaseType::LocalSystemVectorType LocalSystemVectorType
Definition: dam_UP_scheme.hpp:38
Kratos::DamUPScheme::DamUPScheme
DamUPScheme(double beta, double gamma, double rayleigh_m, double rayleigh_k)
Constructor.
Definition: dam_UP_scheme.hpp:44
Kratos::DamUPScheme::mAccelerationVector
std::vector< Vector > mAccelerationVector
Definition: dam_UP_scheme.hpp:540
Kratos::DamUPScheme::Calculate_RHS_Contribution
void Calculate_RHS_Contribution(Element &rCurrentElement, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, const ProcessInfo &CurrentProcessInfo)
Definition: dam_UP_scheme.hpp:411
Kratos::DamUPScheme::~DamUPScheme
virtual ~DamUPScheme()
Destructor.
Definition: dam_UP_scheme.hpp:63
Kratos::DamUPScheme::BaseType
Scheme< TSparseSpace, TDenseSpace > BaseType
Definition: dam_UP_scheme.hpp:34
Kratos::DamUPScheme::InitializeNonLinIteration
void InitializeNonLinIteration(ModelPart &r_model_part, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
unction to be called when it is needed to initialize an iteration. It is designed to be called at the...
Definition: dam_UP_scheme.hpp:289
Kratos::DamUPScheme::Calculate_LHS_Contribution
void Calculate_LHS_Contribution(Condition &rCurrentCondition, LocalSystemMatrixType &LHS_Contribution, Element::EquationIdVectorType &EquationId, const ProcessInfo &CurrentProcessInfo)
Definition: dam_UP_scheme.hpp:485
Kratos::DamUPScheme::LocalSystemMatrixType
BaseType::LocalSystemMatrixType LocalSystemMatrixType
Definition: dam_UP_scheme.hpp:39
Kratos::DamUPScheme::mVelocityVector
std::vector< Vector > mVelocityVector
Definition: dam_UP_scheme.hpp:542
Kratos::DamUPScheme::mGamma
double mGamma
Definition: dam_UP_scheme.hpp:534
Kratos::DamUPScheme::mDeltaTime
double mDeltaTime
Definition: dam_UP_scheme.hpp:535
Kratos::DamUPScheme::CalculateSystemContributions
void CalculateSystemContributions(Element &rCurrentElement, LocalSystemMatrixType &LHS_Contribution, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, const ProcessInfo &CurrentProcessInfo) override
This function is designed to be called in the builder and solver to introduce the selected time integ...
Definition: dam_UP_scheme.hpp:361
Kratos::DamUPScheme::AddDynamicsToRHS
void AddDynamicsToRHS(Element &rCurrentElement, LocalSystemVectorType &RHS_Contribution, LocalSystemMatrixType &M, LocalSystemMatrixType &C, const ProcessInfo &CurrentProcessInfo)
Definition: dam_UP_scheme.hpp:608
Kratos::DamUPScheme::mrayleigh_m
double mrayleigh_m
Definition: dam_UP_scheme.hpp:536
Kratos::DamUPScheme::Check
int Check(ModelPart &r_model_part) override
Definition: dam_UP_scheme.hpp:67
Kratos::DamUPScheme::UpdateVariablesDerivatives
void UpdateVariablesDerivatives(ModelPart &r_model_part)
Definition: dam_UP_scheme.hpp:546
Kratos::DamUPScheme::mBeta
double mBeta
Member Variables.
Definition: dam_UP_scheme.hpp:533
Kratos::DamUPScheme::Predict
void Predict(ModelPart &r_model_part, DofsArrayType &rDofSet, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
Performing the prediction of the solution.
Definition: dam_UP_scheme.hpp:191
Kratos::DamUPScheme::TSystemMatrixType
BaseType::TSystemMatrixType TSystemMatrixType
Definition: dam_UP_scheme.hpp:36
Kratos::DamUPScheme::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(DamUPScheme)
Kratos::DamUPScheme::mMassMatrix
std::vector< Matrix > mMassMatrix
Definition: dam_UP_scheme.hpp:539
Kratos::DamUPScheme::mrayleigh_k
double mrayleigh_k
Definition: dam_UP_scheme.hpp:537
Kratos::DamUPScheme::CalculateSystemContributions
void CalculateSystemContributions(Condition &rCurrentCondition, LocalSystemMatrixType &LHS_Contribution, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, const ProcessInfo &CurrentProcessInfo) override
Functions totally analogous to the precedent but applied to the "condition" objects.
Definition: dam_UP_scheme.hpp:391
Kratos::DamUPScheme::AddDynamicsToLHS
void AddDynamicsToLHS(LocalSystemMatrixType &LHS_Contribution, LocalSystemMatrixType &M, LocalSystemMatrixType &C, const ProcessInfo &CurrentProcessInfo)
Definition: dam_UP_scheme.hpp:595
Kratos::DamUPScheme::Update
void Update(ModelPart &r_model_part, DofsArrayType &rDofSet, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
Performing the update of the solution.
Definition: dam_UP_scheme.hpp:502
Kratos::Element
Base class for all Elements.
Definition: element.h:60
Kratos::Element::GetFirstDerivativesVector
virtual void GetFirstDerivativesVector(Vector &values, int Step=0) const
Definition: element.h:310
Kratos::Element::CalculateLeftHandSide
virtual void CalculateLeftHandSide(MatrixType &rLeftHandSideMatrix, const ProcessInfo &rCurrentProcessInfo)
Definition: element.h:423
Kratos::Element::CalculateMassMatrix
virtual void CalculateMassMatrix(MatrixType &rMassMatrix, const ProcessInfo &rCurrentProcessInfo)
Definition: element.h:570
Kratos::Element::CalculateRightHandSide
virtual void CalculateRightHandSide(VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo)
Definition: element.h:437
Kratos::Element::CalculateDampingMatrix
virtual void CalculateDampingMatrix(MatrixType &rDampingMatrix, const ProcessInfo &rCurrentProcessInfo)
Definition: element.h:583
Kratos::Element::EquationIdVector
virtual void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const
Definition: element.h:258
Kratos::Element::GetSecondDerivativesVector
virtual void GetSecondDerivativesVector(Vector &values, int Step=0) const
Definition: element.h:320
Kratos::Element::CalculateLocalSystem
virtual void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo)
Definition: element.h:405
Kratos::Element::EquationIdVectorType
std::vector< std::size_t > EquationIdVectorType
Definition: element.h:98
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::ElementsBegin
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
Kratos::ModelPart::ConditionsBegin
ConditionIterator ConditionsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1361
Kratos::ModelPart::GetBufferSize
IndexType GetBufferSize() const
This method gets the suffer size of the model part database.
Definition: model_part.h:1876
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::Conditions
ConditionsContainerType & Conditions(IndexType ThisIndex=0)
Definition: model_part.h:1381
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::Elements
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
Kratos::ModelPart::NodesEnd
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
Kratos::ModelPart::Nodes
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
Kratos::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::ParallelUtilities::GetNumThreads
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
Kratos::PointerVectorSet
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
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::Scheme
This class provides the implementation of the basic tasks that are needed by the solution strategy.
Definition: scheme.h:56
Kratos::Scheme::TSystemMatrixType
typename TSparseSpace::MatrixType TSystemMatrixType
Matrix type definition.
Definition: scheme.h:71
Kratos::Scheme::EquationId
virtual void EquationId(const Element &rElement, Element::EquationIdVectorType &rEquationId, const ProcessInfo &rCurrentProcessInfo)
This method gets the eqaution id corresponding to the current element.
Definition: scheme.h:636
Kratos::Scheme::TSystemVectorType
typename TSparseSpace::VectorType TSystemVectorType
Vector type definition.
Definition: scheme.h:74
Kratos::Scheme::LocalSystemVectorType
typename TDenseSpace::VectorType LocalSystemVectorType
Local system vector type definition.
Definition: scheme.h:80
Kratos::Scheme::LocalSystemMatrixType
typename TDenseSpace::MatrixType LocalSystemMatrixType
Local system matrix type definition.
Definition: scheme.h:77
Kratos::Scheme::mSchemeIsInitialized
bool mSchemeIsInitialized
Definition: scheme.h:755
Kratos::array_1d< double, 3 >
dam_application_variables.h
KRATOS_THROW_ERROR
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_WATCH
#define KRATOS_WATCH(variable)
Definition: define.h:806
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
model_part.h
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::block_for_each
void block_for_each(TIterator itBegin, TIterator itEnd, TFunction &&rFunction)
Execute a functor on all items of a range in parallel.
Definition: parallel_utilities.h:299
Kratos::prod
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
generate_hyper_elastic_simo_taylor_neo_hookean.C
int C
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:27
generate_total_lagrangian_mixed_volumetric_strain_element.b
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
generate_two_fluid_navier_stokes.gamma
float gamma
Definition: generate_two_fluid_navier_stokes.py:131
ode_solve.dof
int dof
Definition: ode_solve.py:393
sensitivityMatrix.A
A
Definition: sensitivityMatrix.py:70
tensors::i
integer i
Definition: TensorModule.f:17
Generated on Sat Feb 3 2024 01:59:53 for KratosMultiphysics by doxygen 1.9.1