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.
nodal_two_step_v_p_strategy.h
Go to the documentation of this file.
1 //
2 // Project Name: KratosPFEMFluidDynamicsApplication $
3 // Last modified by: $Author: AFranci $
4 // Date: $Date: June 2018 $
5 // Revision: $Revision: 0.0 $
6 //
7 //
8 
9 #ifndef KRATOS_NODAL_TWO_STEP_V_P_STRATEGY_H
10 #define KRATOS_NODAL_TWO_STEP_V_P_STRATEGY_H
11 
12 #include "includes/define.h"
13 #include "includes/model_part.h"
14 #include "includes/deprecated_variables.h"
15 #include "includes/cfd_variables.h"
16 #include "utilities/openmp_utils.h"
17 #include "processes/process.h"
18 #include "solving_strategies/schemes/scheme.h"
19 #include "solving_strategies/strategies/implicit_solving_strategy.h"
20 #include "custom_utilities/mesher_utilities.hpp"
21 #include "custom_utilities/boundary_normals_calculation_utilities.hpp"
22 #include "geometries/geometry.h"
23 #include "utilities/geometry_utilities.h"
24 
25 #include "solving_strategies/schemes/residualbased_incrementalupdate_static_scheme.h"
26 #include "custom_strategies/builders_and_solvers/nodal_residualbased_elimination_builder_and_solver.h"
27 #include "custom_strategies/builders_and_solvers/nodal_residualbased_elimination_builder_and_solver_continuity.h"
28 
29 #include "custom_utilities/solver_settings.h"
30 
31 #include "custom_strategies/strategies/gauss_seidel_linear_strategy.h"
32 
33 #include "pfem_fluid_dynamics_application_variables.h"
34 
35 #include <stdio.h>
36 #include <cmath>
37 #include <iostream>
38 #include <fstream>
39 
40 namespace Kratos
41 {
42 
45 
48 
52 
54 
57 
61 
65 
66  template <class TSparseSpace,
67  class TDenseSpace,
68  class TLinearSolver>
69  class NodalTwoStepVPStrategy : public ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>
70  {
71  public:
74  KRATOS_CLASS_POINTER_DEFINITION(NodalTwoStepVPStrategy);
75 
77  typedef ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver> BaseType;
78 
79  typedef typename BaseType::TDataType TDataType;
80 
82  typedef Node NodeType;
83 
85  typedef Geometry<NodeType> GeometryType;
86 
87  typedef std::size_t SizeType;
88 
89  typedef typename BaseType::DofsArrayType DofsArrayType;
90 
91  typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
92 
93  typedef typename BaseType::TSystemVectorType TSystemVectorType;
94 
95  typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
96 
97  typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
98 
99  typedef typename BaseType::ElementsArrayType ElementsArrayType;
100 
101  typedef typename ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>::Pointer StrategyPointerType;
102 
103  typedef TwoStepVPSolverSettings<TSparseSpace, TDenseSpace, TLinearSolver> SolverSettingsType;
104 
105  typedef GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType;
106 
107  typedef GlobalPointersVector<Node> NodeWeakPtrVectorType;
111 
112  NodalTwoStepVPStrategy(ModelPart &rModelPart,
113  SolverSettingsType &rSolverConfig) : BaseType(rModelPart)
114  {
115  InitializeStrategy(rSolverConfig);
116  }
117 
118  NodalTwoStepVPStrategy(ModelPart &rModelPart,
119  typename TLinearSolver::Pointer pVelocityLinearSolver,
120  typename TLinearSolver::Pointer pPressureLinearSolver,
121  bool ReformDofSet = true,
122  double VelTol = 0.0001,
123  double PresTol = 0.0001,
124  int MaxPressureIterations = 1, // Only for predictor-corrector
125  unsigned int TimeOrder = 2,
126  unsigned int DomainSize = 2) : BaseType(rModelPart), // Move Mesh flag, pass as input?
127  mVelocityTolerance(VelTol),
128  mPressureTolerance(PresTol),
129  mMaxPressureIter(MaxPressureIterations),
130  mDomainSize(DomainSize),
131  mTimeOrder(TimeOrder),
132  mReformDofSet(ReformDofSet)
133  {
134  KRATOS_TRY;
135 
136  BaseType::SetEchoLevel(1);
137 
138  // Check that input parameters are reasonable and sufficient.
139  this->Check();
140 
141  bool CalculateNormDxFlag = true;
142 
143  bool ReformDofAtEachIteration = false; // DofSet modifiaction is managed by the fractional step strategy, auxiliary strategies should not modify the DofSet directly.
144 
145  // Additional Typedefs
146  typedef typename BuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::Pointer BuilderSolverTypePointer;
147  typedef ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver> BaseType;
148 
149  // initializing fractional velocity solution step
150  typedef Scheme<TSparseSpace, TDenseSpace> SchemeType;
151  typename SchemeType::Pointer pScheme;
152 
153  typename SchemeType::Pointer Temp = typename SchemeType::Pointer(new ResidualBasedIncrementalUpdateStaticScheme<TSparseSpace, TDenseSpace>());
154  pScheme.swap(Temp);
155 
156  // CONSTRUCTION OF VELOCITY
157  BuilderSolverTypePointer vel_build = BuilderSolverTypePointer(new NodalResidualBasedEliminationBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>(pVelocityLinearSolver));
158 
159  this->mpMomentumStrategy = typename BaseType::Pointer(new GaussSeidelLinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, pScheme, pVelocityLinearSolver, vel_build, ReformDofAtEachIteration, CalculateNormDxFlag));
160 
161  this->mpMomentumStrategy->SetEchoLevel(BaseType::GetEchoLevel());
162 
163  vel_build->SetCalculateReactionsFlag(false);
164 
165  BuilderSolverTypePointer pressure_build = BuilderSolverTypePointer(new NodalResidualBasedEliminationBuilderAndSolverContinuity<TSparseSpace, TDenseSpace, TLinearSolver>(pPressureLinearSolver));
166 
167  this->mpPressureStrategy = typename BaseType::Pointer(new GaussSeidelLinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, pScheme, pPressureLinearSolver, pressure_build, ReformDofAtEachIteration, CalculateNormDxFlag));
168 
169  this->mpPressureStrategy->SetEchoLevel(BaseType::GetEchoLevel());
170 
171  pressure_build->SetCalculateReactionsFlag(false);
172 
173  KRATOS_CATCH("");
174  }
175 
177  virtual ~NodalTwoStepVPStrategy() {}
178 
179  int Check() override
180  {
181  KRATOS_TRY;
182 
183  // Check elements and conditions in the model part
184  int ierr = BaseType::Check();
185  if (ierr != 0)
186  return ierr;
187 
188  if (DELTA_TIME.Key() == 0)
189  KRATOS_THROW_ERROR(std::runtime_error, "DELTA_TIME Key is 0. Check that the application was correctly registered.", "");
190 
191  ModelPart &rModelPart = BaseType::GetModelPart();
192 
193  const auto &r_current_process_info = rModelPart.GetProcessInfo();
194  for (const auto &r_element : rModelPart.Elements())
195  {
196  ierr = r_element.Check(r_current_process_info);
197  if (ierr != 0)
198  {
199  break;
200  }
201  }
202 
203  return ierr;
204 
205  KRATOS_CATCH("");
206  }
207 
208  void SetTimeCoefficients(ProcessInfo &rCurrentProcessInfo)
209  {
210  KRATOS_TRY;
211 
212  if (mTimeOrder == 2)
213  {
214  // calculate the BDF coefficients
215  double Dt = rCurrentProcessInfo[DELTA_TIME];
216  double OldDt = rCurrentProcessInfo.GetPreviousTimeStepInfo(1)[DELTA_TIME];
217 
218  double Rho = OldDt / Dt;
219  double TimeCoeff = 1.0 / (Dt * Rho * Rho + Dt * Rho);
220 
221  Vector &BDFcoeffs = rCurrentProcessInfo[BDF_COEFFICIENTS];
222  BDFcoeffs.resize(3, false);
223 
224  BDFcoeffs[0] = TimeCoeff * (Rho * Rho + 2.0 * Rho); // coefficient for step n+1 (3/2Dt if Dt is constant)
225  BDFcoeffs[1] = -TimeCoeff * (Rho * Rho + 2.0 * Rho + 1.0); // coefficient for step n (-4/2Dt if Dt is constant)
226  BDFcoeffs[2] = TimeCoeff; // coefficient for step n-1 (1/2Dt if Dt is constant)
227  }
228  else if (mTimeOrder == 1)
229  {
230  double Dt = rCurrentProcessInfo[DELTA_TIME];
231  double TimeCoeff = 1.0 / Dt;
232 
233  Vector &BDFcoeffs = rCurrentProcessInfo[BDF_COEFFICIENTS];
234  BDFcoeffs.resize(2, false);
235 
236  BDFcoeffs[0] = TimeCoeff; // coefficient for step n+1 (1/Dt)
237  BDFcoeffs[1] = -TimeCoeff; // coefficient for step n (-1/Dt)
238  }
239 
240  KRATOS_CATCH("");
241  }
242 
243  bool SolveSolutionStep() override
244  {
245  // Initialize BDF2 coefficients
246  ModelPart &rModelPart = BaseType::GetModelPart();
247  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
248  double currentTime = rCurrentProcessInfo[TIME];
249  double timeInterval = rCurrentProcessInfo[DELTA_TIME];
250  bool timeIntervalChanged = rCurrentProcessInfo[TIME_INTERVAL_CHANGED];
251  bool converged = false;
252 
253  unsigned int maxNonLinearIterations = mMaxPressureIter;
254 
255  KRATOS_INFO("\n Solve with nodally_integrated_two_step_vp strategy at t=") << currentTime << "s" << std::endl;
256 
257  if (timeIntervalChanged == true && currentTime > 10 * timeInterval)
258  {
259  maxNonLinearIterations *= 2;
260  }
261  if (currentTime < 10 * timeInterval)
262  {
263  if (BaseType::GetEchoLevel() > 1)
264  std::cout << "within the first 10 time steps, I consider the given iteration number x3" << std::endl;
265  maxNonLinearIterations *= 3;
266  }
267  if (currentTime < 20 * timeInterval && currentTime >= 10 * timeInterval)
268  {
269  if (BaseType::GetEchoLevel() > 1)
270  std::cout << "within the second 10 time steps, I consider the given iteration number x2" << std::endl;
271  maxNonLinearIterations *= 2;
272  }
273 
274  bool momentumConverged = true;
275  bool continuityConverged = false;
276  bool fixedTimeStep = false;
277  double pressureNorm = 0;
278  double velocityNorm = 0;
279 
280  this->FillNodalSFDVector();
281  for (unsigned int it = 0; it < maxNonLinearIterations; ++it)
282  {
283  if (BaseType::GetEchoLevel() > 1 && rModelPart.GetCommunicator().MyPID() == 0)
284  std::cout << "----- > iteration: " << it << std::endl;
285 
286  if (it == 0)
287  {
288  this->ComputeNodalVolume();
289 
290  this->InitializeNonLinearIterations();
291  }
292 
293  this->CalcNodalStrainsAndStresses();
294 
295  momentumConverged = this->SolveMomentumIteration(it, maxNonLinearIterations, fixedTimeStep, velocityNorm);
296 
297  this->UpdateTopology(rModelPart, BaseType::GetEchoLevel());
298  this->ComputeNodalVolume();
299  this->InitializeNonLinearIterations();
300  this->CalcNodalStrains();
301 
302  if (fixedTimeStep == false)
303  {
304  continuityConverged = this->SolveContinuityIteration(it, maxNonLinearIterations, pressureNorm);
305  }
306 
307  // if((momentumConverged==true || it==maxNonLinearIterations-1) && momentumAlreadyConverged==false){
308  // std::ofstream myfile;
309  // myfile.open ("momentumConvergedIteration.txt",std::ios::app);
310  // myfile << currentTime << "\t" << it << "\n";
311  // myfile.close();
312  // momentumAlreadyConverged=true;
313  // }
314  // if((continuityConverged==true || it==maxNonLinearIterations-1) && continuityAlreadyConverged==false){
315  // std::ofstream myfile;
316  // myfile.open ("continuityConvergedIteration.txt",std::ios::app);
317  // myfile << currentTime << "\t" << it << "\n";
318  // myfile.close();
319  // continuityAlreadyConverged=true;
320  // }
321 
322  if (it == maxNonLinearIterations - 1 || ((continuityConverged && momentumConverged) && it > 1))
323  {
324  // this->ComputeErrorL2NormCaseImposedG();
325  // this->ComputeErrorL2NormCasePoiseuille();
326  this->CalculateAccelerations();
327  // std::ofstream myfile;
328  // myfile.open ("maxConvergedIteration.txt",std::ios::app);
329  // myfile << currentTime << "\t" << it << "\n";
330  // myfile.close();
331  }
332  if ((continuityConverged && momentumConverged) && it > 1)
333  {
334  rCurrentProcessInfo.SetValue(BAD_VELOCITY_CONVERGENCE, false);
335  rCurrentProcessInfo.SetValue(BAD_PRESSURE_CONVERGENCE, false);
336  converged = true;
337  std::cout << "nodal V-P strategy converged in " << it + 1 << " iterations." << std::endl;
338  break;
339  }
340  }
341 
342  if (!continuityConverged && !momentumConverged && BaseType::GetEchoLevel() > 0 && rModelPart.GetCommunicator().MyPID() == 0)
343  std::cout << "Convergence tolerance not reached." << std::endl;
344 
345  if (mReformDofSet)
346  this->Clear();
347 
348  return converged;
349  }
350 
351  void FinalizeSolutionStep() override
352  {
353  /* this->UpdateStressStrain(); */
354  }
355 
356  void Initialize() override
357  {
358 
359  ModelPart &rModelPart = BaseType::GetModelPart();
360  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
361  unsigned int sizeStrains = 3 * (dimension - 1);
362 
363  // #pragma omp parallel
364  // {
365  ModelPart::NodeIterator NodesBegin;
366  ModelPart::NodeIterator NodesEnd;
367  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
368 
369  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
370  {
371  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
372  unsigned int neighbourNodes = neighb_nodes.size();
373  unsigned int sizeSDFNeigh = neighbourNodes * dimension;
374 
375  if (itNode->SolutionStepsDataHas(NODAL_CAUCHY_STRESS))
376  {
377  Vector &rNodalStress = itNode->FastGetSolutionStepValue(NODAL_CAUCHY_STRESS);
378  if (rNodalStress.size() != sizeStrains)
379  {
380  rNodalStress.resize(sizeStrains, false);
381  }
382  noalias(rNodalStress) = ZeroVector(sizeStrains);
383  }
384  else
385  {
386  std::cout << "THIS node does not have NODAL_CAUCHY_STRESS... " << itNode->X() << " " << itNode->Y() << std::endl;
387  }
388 
389  if (itNode->SolutionStepsDataHas(NODAL_DEVIATORIC_CAUCHY_STRESS))
390  {
391  Vector &rNodalStress = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS);
392  if (rNodalStress.size() != sizeStrains)
393  {
394  rNodalStress.resize(sizeStrains, false);
395  }
396  noalias(rNodalStress) = ZeroVector(sizeStrains);
397  }
398  else
399  {
400  std::cout << "THIS node does not have NODAL_DEVIATORIC_CAUCHY_STRESS... " << itNode->X() << " " << itNode->Y() << std::endl;
401  }
402 
403  if (itNode->SolutionStepsDataHas(NODAL_VOLUME))
404  {
405  itNode->FastGetSolutionStepValue(NODAL_VOLUME) = 0;
406  }
407  else
408  {
409  std::cout << "THIS node does not have NODAL_VOLUME... " << itNode->X() << " " << itNode->Y() << std::endl;
410  }
411 
412  if (itNode->SolutionStepsDataHas(NODAL_MEAN_MESH_SIZE))
413  {
414  itNode->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE) = 0;
415  }
416  else
417  {
418  std::cout << "THIS node does not have NODAL_MEAN_MESH_SIZE... " << itNode->X() << " " << itNode->Y() << std::endl;
419  }
420 
421  if (itNode->SolutionStepsDataHas(NODAL_FREESURFACE_AREA))
422  {
423  itNode->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA) = 0;
424  }
425  else
426  {
427  std::cout << "THIS node does not have NODAL_FREESURFACE_AREA... " << itNode->X() << " " << itNode->Y() << std::endl;
428  }
429 
430  if (itNode->SolutionStepsDataHas(NODAL_SFD_NEIGHBOURS))
431  {
432  Vector &rNodalSFDneighbours = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS);
433  if (rNodalSFDneighbours.size() != sizeSDFNeigh)
434  {
435  rNodalSFDneighbours.resize(sizeSDFNeigh, false);
436  }
437  noalias(rNodalSFDneighbours) = ZeroVector(sizeSDFNeigh);
438  }
439  else
440  {
441  std::cout << "THIS node does not have NODAL_SFD_NEIGHBOURS... " << itNode->X() << " " << itNode->Y() << std::endl;
442  }
443 
444  if (itNode->SolutionStepsDataHas(NODAL_SPATIAL_DEF_RATE))
445  {
446  Vector &rSpatialDefRate = itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
447  if (rSpatialDefRate.size() != sizeStrains)
448  {
449  rSpatialDefRate.resize(sizeStrains, false);
450  }
451  noalias(rSpatialDefRate) = ZeroVector(sizeStrains);
452  }
453  else
454  {
455  std::cout << "THIS node does not have NODAL_SPATIAL_DEF_RATE... " << itNode->X() << " " << itNode->Y() << std::endl;
456  }
457 
458  if (itNode->SolutionStepsDataHas(NODAL_DEFORMATION_GRAD))
459  {
460  Matrix &rFgrad = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD);
461  if (rFgrad.size1() != dimension)
462  {
463  rFgrad.resize(dimension, dimension, false);
464  }
465  noalias(rFgrad) = ZeroMatrix(dimension, dimension);
466  }
467  else
468  {
469  std::cout << "THIS node does not have NODAL_DEFORMATION_GRAD... " << itNode->X() << " " << itNode->Y() << std::endl;
470  }
471 
472  if (itNode->SolutionStepsDataHas(NODAL_DEFORMATION_GRAD_VEL))
473  {
474  Matrix &rFgradVel = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD_VEL);
475  if (rFgradVel.size1() != dimension)
476  {
477  rFgradVel.resize(dimension, dimension, false);
478  }
479  noalias(rFgradVel) = ZeroMatrix(dimension, dimension);
480  }
481  else
482  {
483  std::cout << "THIS node does not have NODAL_DEFORMATION_GRAD_VEL... " << itNode->X() << " " << itNode->Y() << std::endl;
484  }
485 
486  this->InitialAssignFluidMaterialToEachNode(itNode);
487  }
488 
489  // }
490  }
491 
492  void InitialAssignFluidMaterialToEachNode(ModelPart::NodeIterator itNode)
493  {
494  ModelPart &rModelPart = BaseType::GetModelPart();
495  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
496  const double timeInterval = rCurrentProcessInfo[DELTA_TIME];
497 
498  double deviatoricCoeff = 0;
499  if (rModelPart.GetNodalSolutionStepVariablesList().Has(DYNAMIC_VISCOSITY))
500  {
501  deviatoricCoeff = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
502  }
503 
504  double volumetricCoeff = timeInterval * itNode->FastGetSolutionStepValue(BULK_MODULUS);
505  double currFirstLame = volumetricCoeff - 2.0 * deviatoricCoeff / 3.0;
506 
507  itNode->FastGetSolutionStepValue(VOLUMETRIC_COEFFICIENT) = currFirstLame;
508  itNode->FastGetSolutionStepValue(DEVIATORIC_COEFFICIENT) = deviatoricCoeff;
509  }
510 
511  void ComputeNodalVolume()
512  {
513  ModelPart &rModelPart = BaseType::GetModelPart();
514  ElementsArrayType &pElements = rModelPart.Elements();
515  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
516 
517 #ifdef _OPENMP
518  int number_of_threads = omp_get_max_threads();
519 #else
520  int number_of_threads = 1;
521 #endif
522 
523  vector<unsigned int> element_partition;
524  OpenMPUtils::CreatePartition(number_of_threads, pElements.size(), element_partition);
525 
526  // #pragma omp parallel
527  // {
528  int k = OpenMPUtils::ThisThread();
529  typename ElementsArrayType::iterator ElemBegin = pElements.begin() + element_partition[k];
530  typename ElementsArrayType::iterator ElemEnd = pElements.begin() + element_partition[k + 1];
531 
532  for (typename ElementsArrayType::iterator itElem = ElemBegin; itElem != ElemEnd; itElem++) // MSI: To be parallelized
533  {
534  Element::GeometryType &geometry = itElem->GetGeometry();
535  double elementalVolume = 0;
536 
537  if (dimension == 2)
538  {
539  elementalVolume = geometry.Area() / 3.0;
540  }
541  else if (dimension == 3)
542  {
543  elementalVolume = geometry.Volume() * 0.25;
544  }
545  // index = 0;
546  unsigned int numNodes = geometry.size();
547  for (unsigned int i = 0; i < numNodes; i++)
548  {
549  double &nodalVolume = geometry(i)->FastGetSolutionStepValue(NODAL_VOLUME);
550  nodalVolume += elementalVolume;
551  }
552  }
553  // }
554  }
555 
556  void SetNeighboursOrderToNode(ModelPart::NodeIterator itNode)
557  {
558  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
559  unsigned int neighbourNodes = neighb_nodes.size() + 1; // +1 becausealso the node itself must be considered as nieghbor node
560  Vector &rNodeOrderedNeighbours = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS_ORDER);
561 
562  if (rNodeOrderedNeighbours.size() != neighbourNodes)
563  rNodeOrderedNeighbours.resize(neighbourNodes, false);
564 
565  noalias(rNodeOrderedNeighbours) = ZeroVector(neighbourNodes);
566 
567  rNodeOrderedNeighbours[0] = itNode->Id();
568 
569  if (neighbourNodes > 1)
570  {
571  for (unsigned int k = 0; k < neighbourNodes - 1; k++)
572  {
573  rNodeOrderedNeighbours[k + 1] = neighb_nodes[k].Id();
574  }
575  }
576  }
577 
578  void InitializeNodalVariablesForRemeshedDomain(ModelPart::NodeIterator itNode)
579  {
580 
581  ModelPart &rModelPart = BaseType::GetModelPart();
582  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
583  unsigned int sizeStrains = 3 * (dimension - 1);
584  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
585  unsigned int neighbourNodes = neighb_nodes.size() + 1;
586  unsigned int sizeSDFNeigh = neighbourNodes * dimension;
587 
588  if (itNode->SolutionStepsDataHas(NODAL_CAUCHY_STRESS))
589  {
590  Vector &rNodalStress = itNode->FastGetSolutionStepValue(NODAL_CAUCHY_STRESS);
591  if (rNodalStress.size() != sizeStrains)
592  rNodalStress.resize(sizeStrains, false);
593  noalias(rNodalStress) = ZeroVector(sizeStrains);
594  }
595  if (itNode->SolutionStepsDataHas(NODAL_DEVIATORIC_CAUCHY_STRESS))
596  {
597  Vector &rNodalDevStress = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS);
598  if (rNodalDevStress.size() != sizeStrains)
599  rNodalDevStress.resize(sizeStrains, false);
600  noalias(rNodalDevStress) = ZeroVector(sizeStrains);
601  }
602  if (itNode->SolutionStepsDataHas(NODAL_SFD_NEIGHBOURS_ORDER))
603  {
604  Vector &rNodalSFDneighboursOrder = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS_ORDER);
605  if (rNodalSFDneighboursOrder.size() != neighbourNodes)
606  rNodalSFDneighboursOrder.resize(neighbourNodes, false);
607  noalias(rNodalSFDneighboursOrder) = ZeroVector(neighbourNodes);
608  }
609  if (itNode->SolutionStepsDataHas(NODAL_SFD_NEIGHBOURS))
610  {
611  Vector &rNodalSFDneighbours = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS);
612  if (rNodalSFDneighbours.size() != sizeSDFNeigh)
613  rNodalSFDneighbours.resize(sizeSDFNeigh, false);
614  noalias(rNodalSFDneighbours) = ZeroVector(sizeSDFNeigh);
615  }
616  if (itNode->SolutionStepsDataHas(NODAL_SPATIAL_DEF_RATE))
617  {
618  Vector &rSpatialDefRate = itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
619  if (rSpatialDefRate.size() != sizeStrains)
620  rSpatialDefRate.resize(sizeStrains, false);
621  noalias(rSpatialDefRate) = ZeroVector(sizeStrains);
622  }
623  if (itNode->SolutionStepsDataHas(NODAL_DEFORMATION_GRAD))
624  {
625  Matrix &rFgrad = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD);
626  if (rFgrad.size1() != dimension)
627  rFgrad.resize(dimension, dimension, false);
628  noalias(rFgrad) = ZeroMatrix(dimension, dimension);
629  }
630  if (itNode->SolutionStepsDataHas(NODAL_DEFORMATION_GRAD_VEL))
631  {
632  Matrix &rFgradVel = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD_VEL);
633  if (rFgradVel.size1() != dimension)
634  rFgradVel.resize(dimension, dimension, false);
635  noalias(rFgradVel) = ZeroMatrix(dimension, dimension);
636  }
637  if (itNode->SolutionStepsDataHas(NODAL_VOLUME))
638  {
639  itNode->FastGetSolutionStepValue(NODAL_VOLUME) = 0;
640  }
641  if (itNode->SolutionStepsDataHas(NODAL_MEAN_MESH_SIZE))
642  {
643  itNode->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE) = 0;
644  }
645  if (itNode->SolutionStepsDataHas(NODAL_FREESURFACE_AREA))
646  {
647  itNode->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA) = 0;
648  }
649  if (itNode->SolutionStepsDataHas(NODAL_VOLUMETRIC_DEF_RATE))
650  {
651  itNode->FastGetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = 0;
652  }
653  if (itNode->SolutionStepsDataHas(NODAL_EQUIVALENT_STRAIN_RATE))
654  {
655  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = 0;
656  }
657  }
658 
659  void InitializeNonLinearIterations()
660  {
661 
662  ModelPart &rModelPart = BaseType::GetModelPart();
663  ElementsArrayType &pElements = rModelPart.Elements();
664  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
665 
666 #ifdef _OPENMP
667  int number_of_threads = omp_get_max_threads();
668 #else
669  int number_of_threads = 1;
670 #endif
671 
672  vector<unsigned int> element_partition;
673  OpenMPUtils::CreatePartition(number_of_threads, pElements.size(), element_partition);
674 
675  // #pragma omp parallel
676  // {
677  int k = OpenMPUtils::ThisThread();
678  typename ElementsArrayType::iterator ElemBegin = pElements.begin() + element_partition[k];
679  typename ElementsArrayType::iterator ElemEnd = pElements.begin() + element_partition[k + 1];
680 
681  for (typename ElementsArrayType::iterator itElem = ElemBegin; itElem != ElemEnd; itElem++) // MSI: To be parallelized
682  {
683  itElem->InitializeNonLinearIteration(rCurrentProcessInfo);
684  }
685  // }
686  }
687 
688  void CalcNodalStrainsAndStresses()
689  {
690  ModelPart &rModelPart = BaseType::GetModelPart();
691 
692  // #pragma omp parallel
693  // {
694  ModelPart::NodeIterator NodesBegin;
695  ModelPart::NodeIterator NodesEnd;
696  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
697 
698  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
699  {
700 
701  double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
702 
703  double theta = 0.5;
704 
705  if (nodalVolume > 0)
706  {
707  this->ComputeAndStoreNodalDeformationGradient(itNode, theta);
708  this->CalcNodalStrainsAndStressesForNode(itNode);
709  }
710  else
711  { // if nodalVolume==0
712  InitializeNodalVariablesForRemeshedDomain(itNode);
713  }
714  }
715  // }
716  }
717 
718  void CalcNodalStrainsAndStressesForNode(ModelPart::NodeIterator itNode)
719  {
720 
721  ModelPart &rModelPart = BaseType::GetModelPart();
722  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
723  double currFirstLame = itNode->FastGetSolutionStepValue(VOLUMETRIC_COEFFICIENT);
724  double deviatoricCoeff = 0;
725 
726  Matrix Fgrad = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD);
727  Matrix FgradVel = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD_VEL);
728  double detFgrad = 1.0;
729  Matrix InvFgrad = ZeroMatrix(dimension, dimension);
730  Matrix SpatialVelocityGrad = ZeroMatrix(dimension, dimension);
731 
732  if (dimension == 2)
733  {
734  MathUtils<double>::InvertMatrix2(Fgrad, InvFgrad, detFgrad);
735  }
736  else if (dimension == 3)
737  {
738  MathUtils<double>::InvertMatrix3(Fgrad, InvFgrad, detFgrad);
739  }
740 
741  // it computes the spatial velocity gradient tensor --> [L_ij]=dF_ik*invF_kj
742  SpatialVelocityGrad = prod(FgradVel, InvFgrad);
743 
744  if (dimension == 2)
745  {
746  auto &r_stain_tensor2D = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
747  r_stain_tensor2D[0] = SpatialVelocityGrad(0, 0);
748  r_stain_tensor2D[1] = SpatialVelocityGrad(1, 1);
749  r_stain_tensor2D[2] = 0.5 * (SpatialVelocityGrad(1, 0) + SpatialVelocityGrad(0, 1));
750 
751  this->ComputeDeviatoricCoefficientForFluid(itNode, dimension, deviatoricCoeff);
752 
753  double DefVol = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] + itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1];
754 
755  itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = DefVol;
756 
757  double nodalSigmaTot_xx = currFirstLame * DefVol + 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0];
758  double nodalSigmaTot_yy = currFirstLame * DefVol + 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1];
759  double nodalSigmaTot_xy = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2];
760 
761  double nodalSigmaDev_xx = 2.0 * deviatoricCoeff * (itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] - DefVol / 3.0);
762  double nodalSigmaDev_yy = 2.0 * deviatoricCoeff * (itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] - DefVol / 3.0);
763  double nodalSigmaDev_xy = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2];
764 
765  auto &r_stress_tensor2D = itNode->GetSolutionStepValue(NODAL_CAUCHY_STRESS, 0);
766  r_stress_tensor2D[0] = nodalSigmaTot_xx;
767  r_stress_tensor2D[1] = nodalSigmaTot_yy;
768  r_stress_tensor2D[2] = nodalSigmaTot_xy;
769 
770  auto &r_dev_stress_tensor2D = itNode->GetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS, 0);
771  r_dev_stress_tensor2D[0] = nodalSigmaDev_xx;
772  r_dev_stress_tensor2D[1] = nodalSigmaDev_yy;
773  r_dev_stress_tensor2D[2] = nodalSigmaDev_xy;
774  }
775  else if (dimension == 3)
776  {
777  auto &r_stain_tensor3D = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
778  r_stain_tensor3D[0] = SpatialVelocityGrad(0, 0);
779  r_stain_tensor3D[1] = SpatialVelocityGrad(1, 1);
780  r_stain_tensor3D[2] = SpatialVelocityGrad(2, 2);
781  r_stain_tensor3D[3] = 0.5 * (SpatialVelocityGrad(1, 0) + SpatialVelocityGrad(0, 1));
782  r_stain_tensor3D[4] = 0.5 * (SpatialVelocityGrad(2, 0) + SpatialVelocityGrad(0, 2));
783  r_stain_tensor3D[5] = 0.5 * (SpatialVelocityGrad(2, 1) + SpatialVelocityGrad(1, 2));
784 
785  this->ComputeDeviatoricCoefficientForFluid(itNode, dimension, deviatoricCoeff);
786 
787  double DefVol = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] + itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] + itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2];
788  itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = DefVol;
789 
790  double nodalSigmaTot_xx = currFirstLame * DefVol + 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0];
791  double nodalSigmaTot_yy = currFirstLame * DefVol + 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1];
792  double nodalSigmaTot_zz = currFirstLame * DefVol + 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2];
793  double nodalSigmaTot_xy = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3];
794  double nodalSigmaTot_xz = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4];
795  double nodalSigmaTot_yz = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5];
796 
797  double nodalSigmaDev_xx = 2.0 * deviatoricCoeff * (itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] - DefVol / 3.0);
798  double nodalSigmaDev_yy = 2.0 * deviatoricCoeff * (itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] - DefVol / 3.0);
799  double nodalSigmaDev_zz = 2.0 * deviatoricCoeff * (itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] - DefVol / 3.0);
800  double nodalSigmaDev_xy = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3];
801  double nodalSigmaDev_xz = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4];
802  double nodalSigmaDev_yz = 2.0 * deviatoricCoeff * itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5];
803 
804  auto &r_stress_tensor3D = itNode->GetSolutionStepValue(NODAL_CAUCHY_STRESS, 0);
805  r_stress_tensor3D[0] = nodalSigmaTot_xx;
806  r_stress_tensor3D[1] = nodalSigmaTot_yy;
807  r_stress_tensor3D[2] = nodalSigmaTot_zz;
808  r_stress_tensor3D[3] = nodalSigmaTot_xy;
809  r_stress_tensor3D[4] = nodalSigmaTot_xz;
810  r_stress_tensor3D[5] = nodalSigmaTot_yz;
811 
812  auto &r_dev_stress_tensor3D = itNode->GetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS, 0);
813  r_dev_stress_tensor3D[0] = nodalSigmaDev_xx;
814  r_dev_stress_tensor3D[1] = nodalSigmaDev_yy;
815  r_dev_stress_tensor3D[2] = nodalSigmaDev_zz;
816  r_dev_stress_tensor3D[3] = nodalSigmaDev_xy;
817  r_dev_stress_tensor3D[4] = nodalSigmaDev_xz;
818  r_dev_stress_tensor3D[5] = nodalSigmaDev_yz;
819  }
820  }
821 
822  void ComputeDeviatoricCoefficientForFluid(ModelPart::NodeIterator itNode, const unsigned int dimension, double &deviatoricCoefficient)
823  {
824  ModelPart &rModelPart = BaseType::GetModelPart();
825  const double tolerance = 1e-12;
826 
827  if (rModelPart.GetNodalSolutionStepVariablesList().Has(STATIC_FRICTION)) // mu(I)-rheology
828  {
829  const double static_friction = itNode->FastGetSolutionStepValue(STATIC_FRICTION);
830  const double dynamic_friction = itNode->FastGetSolutionStepValue(DYNAMIC_FRICTION);
831  const double delta_friction = dynamic_friction - static_friction;
832  const double inertial_number_zero = itNode->FastGetSolutionStepValue(INERTIAL_NUMBER_ZERO);
833  const double grain_diameter = itNode->FastGetSolutionStepValue(GRAIN_DIAMETER);
834  const double grain_density = itNode->FastGetSolutionStepValue(GRAIN_DENSITY);
835  const double regularization_coeff = itNode->FastGetSolutionStepValue(REGULARIZATION_COEFFICIENT);
836 
837  const double theta = 0.5;
838  double mean_pressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) * theta + itNode->FastGetSolutionStepValue(PRESSURE, 1) * (1 - theta);
839 
840  double pressure_tolerance = -1.0e-07;
841  if (mean_pressure > pressure_tolerance)
842  {
843  mean_pressure = pressure_tolerance;
844  }
845 
846  if (dimension == 2)
847  {
848  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt((2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
849  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
850  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2]));
851  }
852  else if (dimension == 3)
853  {
854  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt(2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
855  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
856  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] +
857  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] +
858  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] +
859  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5]);
860  }
861  const double equivalent_strain_rate = itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE);
862  const double exponent = -equivalent_strain_rate / regularization_coeff;
863  const double second_viscous_term = delta_friction * grain_diameter / (inertial_number_zero * std::sqrt(std::fabs(mean_pressure) / grain_density) + equivalent_strain_rate * grain_diameter);
864 
865  if (std::fabs(equivalent_strain_rate) > tolerance)
866  {
867  const double first_viscous_term = static_friction * (1 - std::exp(exponent)) / equivalent_strain_rate;
868  deviatoricCoefficient = (first_viscous_term + second_viscous_term) * std::fabs(mean_pressure);
869  }
870  else
871  {
872  deviatoricCoefficient = 1.0; // this is for the first iteration and first time step
873  }
874  }
875  else if (rModelPart.GetNodalSolutionStepVariablesList().Has(INTERNAL_FRICTION_ANGLE)) // frictiional viscoplastic model
876  {
877  const double dynamic_viscosity = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
878  const double friction_angle = itNode->FastGetSolutionStepValue(INTERNAL_FRICTION_ANGLE);
879  const double cohesion = itNode->FastGetSolutionStepValue(COHESION);
880  const double adaptive_exponent = itNode->FastGetSolutionStepValue(ADAPTIVE_EXPONENT);
881 
882  const double theta = 0.5;
883  double mean_pressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) * theta + itNode->FastGetSolutionStepValue(PRESSURE, 1) * (1 - theta);
884 
885  double pressure_tolerance = -1.0e-07;
886  if (mean_pressure > pressure_tolerance)
887  {
888  mean_pressure = pressure_tolerance;
889  }
890 
891  if (dimension == 2)
892  {
893  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt((2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
894  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
895  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2]));
896  }
897  else if (dimension == 3)
898  {
899  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt(2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
900  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
901  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] +
902  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] +
903  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] +
904  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5]);
905  }
906 
907  const double equivalent_strain_rate = itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE);
908 
909  // Ensuring that the case of equivalent_strain_rate = 0 is not problematic
910  if (std::fabs(equivalent_strain_rate) > tolerance)
911  {
912  const double friction_angle_rad = friction_angle * Globals::Pi / 180.0;
913  const double tanFi = std::tan(friction_angle_rad);
914  double regularization = 1.0 - std::exp(-adaptive_exponent * equivalent_strain_rate);
915  deviatoricCoefficient = dynamic_viscosity + regularization * ((cohesion + tanFi * fabs(mean_pressure)) / equivalent_strain_rate);
916  }
917  else
918  {
919  deviatoricCoefficient = dynamic_viscosity;
920  }
921  }
922  else if (rModelPart.GetNodalSolutionStepVariablesList().Has(YIELD_SHEAR)) // bingham model
923  {
924  if (dimension == 2)
925  {
926  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt((2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
927  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
928  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2]));
929  }
930  else if (dimension == 3)
931  {
932  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt(2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
933  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
934  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] +
935  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] +
936  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] +
937  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5]);
938  }
939 
940  const double yieldShear = itNode->FastGetSolutionStepValue(YIELD_SHEAR);
941  const double equivalentStrainRate = itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE);
942  const double adaptiveExponent = itNode->FastGetSolutionStepValue(ADAPTIVE_EXPONENT);
943  const double exponent = -adaptiveExponent * equivalentStrainRate;
944  deviatoricCoefficient = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
945  if (equivalentStrainRate > tolerance)
946  {
947  deviatoricCoefficient += (yieldShear / equivalentStrainRate) * (1 - exp(exponent));
948  }
949  }
950  else if (rModelPart.GetNodalSolutionStepVariablesList().Has(DYNAMIC_VISCOSITY))
951  {
952  deviatoricCoefficient = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
953  }
954  itNode->FastGetSolutionStepValue(DEVIATORIC_COEFFICIENT) = deviatoricCoefficient;
955  }
956 
957  void CalcNodalStrainsForNode(ModelPart::NodeIterator itNode)
958  {
959 
960  ModelPart &rModelPart = BaseType::GetModelPart();
961  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
962 
963  double detFgrad = 1.0;
964  Matrix nodalFgrad = ZeroMatrix(dimension, dimension);
965  Matrix FgradVel = ZeroMatrix(dimension, dimension);
966  Matrix InvFgrad = ZeroMatrix(dimension, dimension);
967  Matrix SpatialVelocityGrad = ZeroMatrix(dimension, dimension);
968 
969  nodalFgrad = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD);
970  FgradVel = itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD_VEL);
971 
972  // Inverse
973 
974  if (dimension == 2)
975  {
976  MathUtils<double>::InvertMatrix2(nodalFgrad, InvFgrad, detFgrad);
977  }
978  else if (dimension == 3)
979  {
980  MathUtils<double>::InvertMatrix3(nodalFgrad, InvFgrad, detFgrad);
981  }
982 
983  // it computes the spatial velocity gradient tensor --> [L_ij]=dF_ik*invF_kj
984  SpatialVelocityGrad = prod(FgradVel, InvFgrad);
985 
986  if (dimension == 2)
987  {
988 
989  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] = SpatialVelocityGrad(0, 0);
990  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] = SpatialVelocityGrad(1, 1);
991  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] = 0.5 * (SpatialVelocityGrad(1, 0) + SpatialVelocityGrad(0, 1));
992 
993  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt((2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
994  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
995  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2]));
996 
997  double DefX = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0];
998  double DefY = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1];
999 
1000  double DefVol = DefX + DefY;
1001 
1002  itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = DefVol;
1003  }
1004  else if (dimension == 3)
1005  {
1006 
1007  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] = SpatialVelocityGrad(0, 0);
1008  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] = SpatialVelocityGrad(1, 1);
1009  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] = SpatialVelocityGrad(2, 2);
1010  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] = 0.5 * (SpatialVelocityGrad(1, 0) + SpatialVelocityGrad(0, 1));
1011  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] = 0.5 * (SpatialVelocityGrad(2, 0) + SpatialVelocityGrad(0, 2));
1012  itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5] = 0.5 * (SpatialVelocityGrad(2, 1) + SpatialVelocityGrad(1, 2));
1013 
1014  itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = sqrt(2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0] +
1015  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1] +
1016  2.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2] +
1017  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[3] +
1018  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[4] +
1019  4.0 * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5] * itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[5]);
1020 
1021  double DefX = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[0];
1022  double DefY = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[1];
1023  double DefZ = itNode->GetSolutionStepValue(NODAL_SPATIAL_DEF_RATE)[2];
1024 
1025  double DefVol = DefX + DefY + DefZ;
1026 
1027  itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = DefVol;
1028  }
1029  }
1030 
1031  void CalcNodalStrains()
1032  {
1033 
1034  ModelPart &rModelPart = BaseType::GetModelPart();
1035 
1036  // #pragma omp parallel
1037  // {
1038  ModelPart::NodeIterator NodesBegin;
1039  ModelPart::NodeIterator NodesEnd;
1040  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
1041 
1042  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
1043  {
1044 
1045  double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
1046 
1047  double theta = 1.0;
1048 
1049  if (nodalVolume > 0)
1050  {
1051  this->ComputeAndStoreNodalDeformationGradient(itNode, theta);
1052  this->CalcNodalStrainsForNode(itNode);
1053  }
1054  else
1055  {
1056  InitializeNodalVariablesForRemeshedDomain(itNode);
1057  }
1058  }
1059  // }
1060  }
1061 
1062  void ComputeAndStoreNodalDeformationGradient(ModelPart::NodeIterator itNode, double theta)
1063  {
1064 
1065  KRATOS_TRY;
1066 
1067  ModelPart &rModelPart = BaseType::GetModelPart();
1068  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
1069  Vector nodalSFDneighboursId = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS_ORDER);
1070  Vector rNodalSFDneigh = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS);
1071  const unsigned int neighSize = nodalSFDneighboursId.size();
1072  Matrix Fgrad = ZeroMatrix(dimension, dimension);
1073  Matrix FgradVel = ZeroMatrix(dimension, dimension);
1074  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
1075 
1076  if (dimension == 2)
1077  {
1078 
1079  double dNdXi = rNodalSFDneigh[0];
1080  double dNdYi = rNodalSFDneigh[1];
1081 
1082  Fgrad(0, 0) += dNdXi * itNode->X();
1083  Fgrad(0, 1) += dNdYi * itNode->X();
1084  Fgrad(1, 0) += dNdXi * itNode->Y();
1085  Fgrad(1, 1) += dNdYi * itNode->Y();
1086 
1087  double VelocityX = itNode->FastGetSolutionStepValue(VELOCITY_X, 0) * theta + itNode->FastGetSolutionStepValue(VELOCITY_X, 1) * (1 - theta);
1088  double VelocityY = itNode->FastGetSolutionStepValue(VELOCITY_Y, 0) * theta + itNode->FastGetSolutionStepValue(VELOCITY_Y, 1) * (1 - theta);
1089 
1090  FgradVel(0, 0) += dNdXi * VelocityX;
1091  FgradVel(0, 1) += dNdYi * VelocityX;
1092  FgradVel(1, 0) += dNdXi * VelocityY;
1093  FgradVel(1, 1) += dNdYi * VelocityY;
1094 
1095  unsigned int firstRow = 2;
1096 
1097  if (neighSize > 0)
1098  {
1099  for (unsigned int i = 0; i < neighSize - 1; i++) // neigh_nodes has one cell less than nodalSFDneighboursId becuase this has also the considered node ID at the beginning
1100  {
1101  dNdXi = rNodalSFDneigh[firstRow];
1102  dNdYi = rNodalSFDneigh[firstRow + 1];
1103  unsigned int neigh_nodes_id = neighb_nodes[i].Id();
1104  unsigned int other_neigh_nodes_id = nodalSFDneighboursId[i + 1];
1105  if (neigh_nodes_id != other_neigh_nodes_id)
1106  {
1107  std::cout << "node (x,y)=(" << itNode->X() << "," << itNode->Y() << ") with neigh_nodes_id " << neigh_nodes_id << " different than other_neigh_nodes_id " << other_neigh_nodes_id << std::endl;
1108  }
1109  Fgrad(0, 0) += dNdXi * neighb_nodes[i].X();
1110  Fgrad(0, 1) += dNdYi * neighb_nodes[i].X();
1111  Fgrad(1, 0) += dNdXi * neighb_nodes[i].Y();
1112  Fgrad(1, 1) += dNdYi * neighb_nodes[i].Y();
1113 
1114  VelocityX = neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_X, 0) * theta + neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_X, 1) * (1 - theta);
1115  VelocityY = neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_Y, 0) * theta + neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_Y, 1) * (1 - theta);
1116 
1117  FgradVel(0, 0) += dNdXi * VelocityX;
1118  FgradVel(0, 1) += dNdYi * VelocityX;
1119  FgradVel(1, 0) += dNdXi * VelocityY;
1120  FgradVel(1, 1) += dNdYi * VelocityY;
1121 
1122  firstRow += 2;
1123  }
1124  }
1125  }
1126  else
1127  {
1128 
1129  double dNdXi = rNodalSFDneigh[0];
1130  double dNdYi = rNodalSFDneigh[1];
1131  double dNdZi = rNodalSFDneigh[2];
1132 
1133  double VelocityX = itNode->FastGetSolutionStepValue(VELOCITY_X, 0) * theta + itNode->FastGetSolutionStepValue(VELOCITY_X, 1) * (1 - theta);
1134  double VelocityY = itNode->FastGetSolutionStepValue(VELOCITY_Y, 0) * theta + itNode->FastGetSolutionStepValue(VELOCITY_Y, 1) * (1 - theta);
1135  double VelocityZ = itNode->FastGetSolutionStepValue(VELOCITY_Z, 0) * theta + itNode->FastGetSolutionStepValue(VELOCITY_Z, 1) * (1 - theta);
1136 
1137  Fgrad(0, 0) += dNdXi * itNode->X();
1138  Fgrad(0, 1) += dNdYi * itNode->X();
1139  Fgrad(0, 2) += dNdZi * itNode->X();
1140 
1141  Fgrad(1, 0) += dNdXi * itNode->Y();
1142  Fgrad(1, 1) += dNdYi * itNode->Y();
1143  Fgrad(1, 2) += dNdZi * itNode->Y();
1144 
1145  Fgrad(2, 0) += dNdXi * itNode->Z();
1146  Fgrad(2, 1) += dNdYi * itNode->Z();
1147  Fgrad(2, 2) += dNdZi * itNode->Z();
1148 
1149  FgradVel(0, 0) += dNdXi * VelocityX;
1150  FgradVel(0, 1) += dNdYi * VelocityX;
1151  FgradVel(0, 2) += dNdZi * VelocityX;
1152 
1153  FgradVel(1, 0) += dNdXi * VelocityY;
1154  FgradVel(1, 1) += dNdYi * VelocityY;
1155  FgradVel(1, 2) += dNdZi * VelocityY;
1156 
1157  FgradVel(2, 0) += dNdXi * VelocityZ;
1158  FgradVel(2, 1) += dNdYi * VelocityZ;
1159  FgradVel(2, 2) += dNdZi * VelocityZ;
1160 
1161  unsigned int firstRow = 3;
1162 
1163  if (neighSize > 0)
1164  {
1165  for (unsigned int i = 0; i < neighSize - 1; i++)
1166  {
1167 
1168  dNdXi = rNodalSFDneigh[firstRow];
1169  dNdYi = rNodalSFDneigh[firstRow + 1];
1170  dNdZi = rNodalSFDneigh[firstRow + 2];
1171 
1172  VelocityX = neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_X, 0) * theta + neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_X, 1) * (1 - theta);
1173  VelocityY = neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_Y, 0) * theta + neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_Y, 1) * (1 - theta);
1174  VelocityZ = neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_Z, 0) * theta + neighb_nodes[i].FastGetSolutionStepValue(VELOCITY_Z, 1) * (1 - theta);
1175 
1176  Fgrad(0, 0) += dNdXi * neighb_nodes[i].X();
1177  Fgrad(0, 1) += dNdYi * neighb_nodes[i].X();
1178  Fgrad(0, 2) += dNdZi * neighb_nodes[i].X();
1179 
1180  Fgrad(1, 0) += dNdXi * neighb_nodes[i].Y();
1181  Fgrad(1, 1) += dNdYi * neighb_nodes[i].Y();
1182  Fgrad(1, 2) += dNdZi * neighb_nodes[i].Y();
1183 
1184  Fgrad(2, 0) += dNdXi * neighb_nodes[i].Z();
1185  Fgrad(2, 1) += dNdYi * neighb_nodes[i].Z();
1186  Fgrad(2, 2) += dNdZi * neighb_nodes[i].Z();
1187 
1188  FgradVel(0, 0) += dNdXi * VelocityX;
1189  FgradVel(0, 1) += dNdYi * VelocityX;
1190  FgradVel(0, 2) += dNdZi * VelocityX;
1191 
1192  FgradVel(1, 0) += dNdXi * VelocityY;
1193  FgradVel(1, 1) += dNdYi * VelocityY;
1194  FgradVel(1, 2) += dNdZi * VelocityY;
1195 
1196  FgradVel(2, 0) += dNdXi * VelocityZ;
1197  FgradVel(2, 1) += dNdYi * VelocityZ;
1198  FgradVel(2, 2) += dNdZi * VelocityZ;
1199 
1200  firstRow += 3;
1201  }
1202  }
1203  }
1204 
1205  itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD) = Fgrad;
1206  itNode->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD_VEL) = FgradVel;
1207  KRATOS_CATCH("");
1208  }
1209 
1210  void UpdateTopology(ModelPart &rModelPart, unsigned int echoLevel)
1211  {
1212  KRATOS_TRY;
1213  /* this->CalculateDisplacements(); */
1214  this->CalculateDisplacementsAndResetNodalVariables();
1215  BaseType::MoveMesh();
1216  BoundaryNormalsCalculationUtilities BoundaryComputation;
1217  BoundaryComputation.CalculateUnitBoundaryNormals(rModelPart, echoLevel);
1218 
1219  KRATOS_CATCH("");
1220  }
1221 
1222  void CalculatePressureVelocity()
1223  {
1224  ModelPart &rModelPart = BaseType::GetModelPart();
1225  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1226  const double timeInterval = rCurrentProcessInfo[DELTA_TIME];
1227  unsigned int timeStep = rCurrentProcessInfo[STEP];
1228 
1229  for (ModelPart::NodeIterator i = rModelPart.NodesBegin();
1230  i != rModelPart.NodesEnd(); ++i)
1231  {
1232  if (timeStep == 1)
1233  {
1234  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) = 0;
1235  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1) = 0;
1236  }
1237  else
1238  {
1239  double &CurrentPressure = (i)->FastGetSolutionStepValue(PRESSURE, 0);
1240  double &PreviousPressure = (i)->FastGetSolutionStepValue(PRESSURE, 1);
1241  double &CurrentPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0);
1242  CurrentPressureVelocity = (CurrentPressure - PreviousPressure) / timeInterval;
1243  }
1244  }
1245  }
1246 
1247  void CalculatePressureAcceleration()
1248  {
1249  ModelPart &rModelPart = BaseType::GetModelPart();
1250  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1251  const double timeInterval = rCurrentProcessInfo[DELTA_TIME];
1252  unsigned int timeStep = rCurrentProcessInfo[STEP];
1253 
1254  for (ModelPart::NodeIterator i = rModelPart.NodesBegin(); i != rModelPart.NodesEnd(); ++i)
1255  {
1256  if (timeStep == 1)
1257  {
1258  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0) = 0;
1259  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 1) = 0;
1260  }
1261  else
1262  {
1263  double &CurrentPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0);
1264  double &PreviousPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1);
1265  double &CurrentPressureAcceleration = (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0);
1266  CurrentPressureAcceleration = (CurrentPressureVelocity - PreviousPressureVelocity) / timeInterval;
1267  }
1268  }
1269  }
1270 
1271  void CalculateAccelerations()
1272  {
1273  ModelPart &rModelPart = BaseType::GetModelPart();
1274  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1275 
1276  for (ModelPart::NodeIterator i = rModelPart.NodesBegin(); i != rModelPart.NodesEnd(); ++i)
1277  {
1278 
1279  array_1d<double, 3> &CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 0);
1280  array_1d<double, 3> &PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 1);
1281 
1282  array_1d<double, 3> &CurrentAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION, 0);
1283  array_1d<double, 3> &PreviousAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION, 1);
1284 
1285  if ((i)->IsNot(ISOLATED) && ((i)->IsNot(RIGID) || (i)->Is(SOLID)))
1286  {
1287  UpdateAccelerations(CurrentAcceleration, CurrentVelocity, PreviousAcceleration, PreviousVelocity);
1288  }
1289  else if ((i)->Is(RIGID))
1290  {
1291  array_1d<double, 3> Zeros(3, 0.0);
1292  (i)->FastGetSolutionStepValue(ACCELERATION, 0) = Zeros;
1293  (i)->FastGetSolutionStepValue(ACCELERATION, 1) = Zeros;
1294  }
1295  else
1296  {
1297  (i)->FastGetSolutionStepValue(NODAL_VOLUME) = 0.0;
1298  (i)->FastGetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = 0.0;
1299  (i)->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = 0;
1300  (i)->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE) = 0.0;
1301  (i)->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA) = 0.0;
1302  (i)->FastGetSolutionStepValue(PRESSURE, 0) = 0.0;
1303  (i)->FastGetSolutionStepValue(PRESSURE, 1) = 0.0;
1304  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) = 0.0;
1305  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1) = 0.0;
1306  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0) = 0.0;
1307  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 1) = 0.0;
1308  if ((i)->SolutionStepsDataHas(VOLUME_ACCELERATION))
1309  {
1310  array_1d<double, 3> &VolumeAcceleration = (i)->FastGetSolutionStepValue(VOLUME_ACCELERATION);
1311  (i)->FastGetSolutionStepValue(ACCELERATION, 0) = VolumeAcceleration;
1312  (i)->FastGetSolutionStepValue(VELOCITY, 0) += VolumeAcceleration * rCurrentProcessInfo[DELTA_TIME];
1313  }
1314  }
1315  }
1316  }
1317 
1318  inline void UpdateAccelerations(array_1d<double, 3> &CurrentAcceleration,
1319  const array_1d<double, 3> &CurrentVelocity,
1320  array_1d<double, 3> &PreviousAcceleration,
1321  const array_1d<double, 3> &PreviousVelocity)
1322  {
1323  ModelPart &rModelPart = BaseType::GetModelPart();
1324  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1325  double Dt = rCurrentProcessInfo[DELTA_TIME];
1326  noalias(CurrentAcceleration) = 2.0 * (CurrentVelocity - PreviousVelocity) / Dt - PreviousAcceleration;
1327  }
1328 
1329  void CalculateDisplacements()
1330  {
1331  ModelPart &rModelPart = BaseType::GetModelPart();
1332  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1333  const double TimeStep = rCurrentProcessInfo[DELTA_TIME];
1334 
1335  for (ModelPart::NodeIterator i = rModelPart.NodesBegin(); i != rModelPart.NodesEnd(); ++i)
1336  {
1337 
1338  array_1d<double, 3> &CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 0);
1339  array_1d<double, 3> &PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 1);
1340 
1341  array_1d<double, 3> &CurrentDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT, 0);
1342  array_1d<double, 3> &PreviousDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT, 1);
1343 
1344  CurrentDisplacement[0] = 0.5 * TimeStep * (CurrentVelocity[0] + PreviousVelocity[0]) + PreviousDisplacement[0];
1345  CurrentDisplacement[1] = 0.5 * TimeStep * (CurrentVelocity[1] + PreviousVelocity[1]) + PreviousDisplacement[1];
1346  CurrentDisplacement[2] = 0.5 * TimeStep * (CurrentVelocity[2] + PreviousVelocity[2]) + PreviousDisplacement[2];
1347  }
1348  }
1349 
1350  void CalculateDisplacementsAndResetNodalVariables()
1351  {
1352  ModelPart &rModelPart = BaseType::GetModelPart();
1353  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1354  const double TimeStep = rCurrentProcessInfo[DELTA_TIME];
1355  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
1356  unsigned int sizeStrains = 3 * (dimension - 1);
1357 
1358  // #pragma omp parallel
1359  // {
1360  ModelPart::NodeIterator NodesBegin;
1361  ModelPart::NodeIterator NodesEnd;
1362  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
1363 
1364  for (ModelPart::NodeIterator i = NodesBegin; i != NodesEnd; ++i)
1365  {
1366  array_1d<double, 3> &CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 0);
1367  array_1d<double, 3> &PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 1);
1368 
1369  array_1d<double, 3> &CurrentDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT, 0);
1370  array_1d<double, 3> &PreviousDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT, 1);
1371 
1372  CurrentDisplacement[0] = 0.5 * TimeStep * (CurrentVelocity[0] + PreviousVelocity[0]) + PreviousDisplacement[0];
1373  CurrentDisplacement[1] = 0.5 * TimeStep * (CurrentVelocity[1] + PreviousVelocity[1]) + PreviousDisplacement[1];
1374  if (dimension == 3)
1375  {
1376  CurrentDisplacement[2] = 0.5 * TimeStep * (CurrentVelocity[2] + PreviousVelocity[2]) + PreviousDisplacement[2];
1377  }
1378 
1380  Vector &rNodalSFDneighbours = i->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS);
1381  unsigned int sizeSDFNeigh = rNodalSFDneighbours.size();
1382 
1383  i->FastGetSolutionStepValue(NODAL_VOLUME) = 0;
1384  i->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE) = 0;
1385  i->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA) = 0;
1386  i->FastGetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) = 0;
1387  i->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE) = 0;
1388 
1389  noalias(rNodalSFDneighbours) = ZeroVector(sizeSDFNeigh);
1390 
1391  Vector &rSpatialDefRate = i->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
1392  noalias(rSpatialDefRate) = ZeroVector(sizeStrains);
1393 
1394  Matrix &rFgrad = i->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD);
1395  noalias(rFgrad) = ZeroMatrix(dimension, dimension);
1396 
1397  Matrix &rFgradVel = i->FastGetSolutionStepValue(NODAL_DEFORMATION_GRAD_VEL);
1398  noalias(rFgradVel) = ZeroMatrix(dimension, dimension);
1399  }
1400  // }
1401  }
1402 
1403  void UpdatePressureAccelerations()
1404  {
1405  this->CalculateAccelerations();
1406  this->CalculatePressureVelocity();
1407  this->CalculatePressureAcceleration();
1408  }
1409 
1410  void Clear() override
1411  {
1412  mpMomentumStrategy->Clear();
1413  mpPressureStrategy->Clear();
1414  }
1415 
1419 
1420  void SetEchoLevel(int Level) override
1421  {
1422  BaseType::SetEchoLevel(Level);
1423  int StrategyLevel = Level > 0 ? Level - 1 : 0;
1424  mpMomentumStrategy->SetEchoLevel(StrategyLevel);
1425  mpPressureStrategy->SetEchoLevel(StrategyLevel);
1426  }
1427 
1431 
1435 
1437  std::string Info() const override
1438  {
1439  std::stringstream buffer;
1440  buffer << "NodalTwoStepVPStrategy";
1441  return buffer.str();
1442  }
1443 
1445  void PrintInfo(std::ostream &rOStream) const override
1446  {
1447  rOStream << "NodalTwoStepVPStrategy";
1448  }
1449 
1451  void PrintData(std::ostream &rOStream) const override
1452  {
1453  }
1454 
1458 
1460 
1461  protected:
1464 
1468 
1472 
1476 
1480 
1482 
1486  bool SolveMomentumIteration(unsigned int it, unsigned int maxIt, bool &fixedTimeStep, double &velocityNorm)
1487  {
1488  ModelPart &rModelPart = BaseType::GetModelPart();
1489  int Rank = rModelPart.GetCommunicator().MyPID();
1490  bool ConvergedMomentum = false;
1491  double NormDv = 0;
1492  fixedTimeStep = false;
1493  // build momentum system and solve for fractional step velocity increment
1494  rModelPart.GetProcessInfo().SetValue(FRACTIONAL_STEP, 1);
1495 
1496  if (it == 0)
1497  {
1498  mpMomentumStrategy->InitializeSolutionStep();
1499  }
1500 
1501  NormDv = mpMomentumStrategy->Solve();
1502 
1503  if (BaseType::GetEchoLevel() > 1 && Rank == 0)
1504  std::cout << "-------------- s o l v e d ! ------------------" << std::endl;
1505 
1506  if (it == 0)
1507  {
1508  velocityNorm = this->ComputeVelocityNorm();
1509  }
1510  double DvErrorNorm = NormDv / velocityNorm;
1511 
1512  unsigned int iterationForCheck = 3;
1513  KRATOS_INFO("TwoStepVPStrategy") << "iteration(" << it << ") Velocity error: " << DvErrorNorm << std::endl;
1514 
1515  // Check convergence
1516  if (it == maxIt - 1)
1517  {
1518  KRATOS_INFO("Iteration") << it << " Final Velocity error: " << DvErrorNorm << std::endl;
1519  fixedTimeStep = this->FixTimeStepMomentum(DvErrorNorm);
1520  }
1521  else if (it > iterationForCheck)
1522  {
1523  fixedTimeStep = this->CheckMomentumConvergence(DvErrorNorm);
1524  }
1525  // ProcessInfo& rCurrentProcessInfo = rModelPart.GetProcessInfo();
1526  // double currentTime = rCurrentProcessInfo[TIME];
1527  // double tolerance=0.0000000001;
1528  // if(currentTime>(0.25-tolerance) && currentTime<(0.25+tolerance)){
1529  // std::ofstream myfile;
1530  // myfile.open ("velocityConvergenceAt025s.txt",std::ios::app);
1531  // myfile << it << "\t" << DvErrorNorm << "\n";
1532  // myfile.close();
1533  // }
1534  // else if(currentTime>(0.5-tolerance) && currentTime<(0.5+tolerance)){
1535  // std::ofstream myfile;
1536  // myfile.open ("velocityConvergenceAt05s.txt",std::ios::app);
1537  // myfile << it << "\t" << DvErrorNorm << "\n";
1538  // myfile.close();
1539  // }
1540  // else if(currentTime>(0.75-tolerance) && currentTime<(0.75+tolerance)){
1541  // std::ofstream myfile;
1542  // myfile.open ("velocityConvergenceAt075s.txt",std::ios::app);
1543  // myfile << it << "\t" << DvErrorNorm << "\n";
1544  // myfile.close();
1545  // }
1546  // else if(currentTime>(1.0-tolerance) && currentTime<(1.0+tolerance)){
1547  // std::ofstream myfile;
1548  // myfile.open ("velocityConvergenceAt100s.txt",std::ios::app);
1549  // myfile << it << "\t" << DvErrorNorm << "\n";
1550  // myfile.close();
1551  // }
1552 
1553  if (!ConvergedMomentum && BaseType::GetEchoLevel() > 0 && Rank == 0)
1554  std::cout << "Momentum equations did not reach the convergence tolerance." << std::endl;
1555 
1556  return ConvergedMomentum;
1557  }
1558 
1559  bool SolveContinuityIteration(unsigned int it, unsigned int maxIt, double &NormP)
1560  {
1561  ModelPart &rModelPart = BaseType::GetModelPart();
1562  int Rank = rModelPart.GetCommunicator().MyPID();
1563  bool ConvergedContinuity = false;
1564  double NormDp = 0;
1565  // 2. Pressure solution
1566  rModelPart.GetProcessInfo().SetValue(FRACTIONAL_STEP, 5);
1567  if (it == 0)
1568  {
1569  mpPressureStrategy->InitializeSolutionStep();
1570  }
1571 
1572  NormDp = mpPressureStrategy->Solve();
1573 
1574  if (BaseType::GetEchoLevel() > 0 && Rank == 0)
1575  std::cout << "The norm of pressure is: " << NormDp << std::endl;
1576 
1577  if (it == 0)
1578  {
1579  NormP = this->ComputePressureNorm();
1580  }
1581 
1582  double DpErrorNorm = NormDp / (NormP);
1583 
1584  // double DpErrorNorm = 0;
1585  // ConvergedContinuity = this->CheckPressureConvergence(NormDp, DpErrorNorm);
1586  // Check convergence
1587  if (it == maxIt - 1)
1588  {
1589  KRATOS_INFO("Iteration") << it << " Final Pressure error: " << DpErrorNorm << std::endl;
1590  ConvergedContinuity = this->FixTimeStepContinuity(DpErrorNorm);
1591  }
1592  else
1593  {
1594  KRATOS_INFO("Iteration") << it << " Pressure error: " << DpErrorNorm << std::endl;
1595  }
1596 
1597  // ProcessInfo& rCurrentProcessInfo = rModelPart.GetProcessInfo();
1598  // double currentTime = rCurrentProcessInfo[TIME];
1599  // double tolerance=0.0000000001;
1600  // if(currentTime>(0.25-tolerance) && currentTime<(0.25+tolerance)){
1601  // std::ofstream myfile;
1602  // myfile.open ("pressureConvergenceAt025s.txt",std::ios::app);
1603  // myfile << it << "\t" << DpErrorNorm << "\n";
1604  // myfile.close();
1605  // }
1606  // else if(currentTime>(0.5-tolerance) && currentTime<(0.5+tolerance)){
1607  // std::ofstream myfile;
1608  // myfile.open ("pressureConvergenceAt05s.txt",std::ios::app);
1609  // myfile << it << "\t" << DpErrorNorm << "\n";
1610  // myfile.close();
1611  // }
1612  // else if(currentTime>(0.75-tolerance) && currentTime<(0.75+tolerance)){
1613  // std::ofstream myfile;
1614  // myfile.open ("pressureConvergenceAt075s.txt",std::ios::app);
1615  // myfile << it << "\t" << DpErrorNorm << "\n";
1616  // myfile.close();
1617  // }
1618  // else if(currentTime>(1.0-tolerance) && currentTime<(1.0+tolerance)){
1619  // std::ofstream myfile;
1620  // myfile.open ("pressureConvergenceAt100s.txt",std::ios::app);
1621  // myfile << it << "\t" << DpErrorNorm << "\n";
1622  // myfile.close();
1623  // }
1624 
1625  if (!ConvergedContinuity && BaseType::GetEchoLevel() > 0 && Rank == 0)
1626  std::cout << "Continuity equation did not reach the convergence tolerance." << std::endl;
1627 
1628  return ConvergedContinuity;
1629  }
1630 
1631  bool CheckVelocityConvergence(const double NormDv, double &errorNormDv)
1632  {
1633  ModelPart &rModelPart = BaseType::GetModelPart();
1634 
1635  double NormV = 0.00;
1636  errorNormDv = 0;
1637 
1638 #pragma omp parallel reduction(+ : NormV)
1639  {
1640  ModelPart::NodeIterator NodeBegin;
1641  ModelPart::NodeIterator NodeEnd;
1642  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1643  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1644  {
1645  const array_1d<double, 3> &Vel = itNode->FastGetSolutionStepValue(VELOCITY);
1646 
1647  double NormVelNode = 0;
1648 
1649  for (unsigned int d = 0; d < 3; ++d)
1650  {
1651  NormVelNode += Vel[d] * Vel[d];
1652  NormV += Vel[d] * Vel[d];
1653  }
1654  }
1655  }
1656  BaseType::GetModelPart().GetCommunicator().GetDataCommunicator().SumAll(NormV);
1657 
1658  NormV = sqrt(NormV);
1659 
1660  if (NormV == 0.0)
1661  NormV = 1.00;
1662 
1663  errorNormDv = NormDv / NormV;
1664 
1665  if (BaseType::GetEchoLevel() > 0 && rModelPart.GetCommunicator().MyPID() == 0)
1666  {
1667  std::cout << "The norm of velocity increment is: " << NormDv << std::endl;
1668  std::cout << "The norm of velocity is: " << NormV << std::endl;
1669  std::cout << "Velocity error: " << errorNormDv << "mVelocityTolerance: " << mVelocityTolerance << std::endl;
1670  }
1671 
1672  if (errorNormDv < mVelocityTolerance)
1673  {
1674  return true;
1675  }
1676  else
1677  {
1678  return false;
1679  }
1680  }
1681 
1682  double ComputeVelocityNorm()
1683  {
1684  ModelPart &rModelPart = BaseType::GetModelPart();
1685 
1686  double NormV = 0.00;
1687 
1688 #pragma omp parallel reduction(+ : NormV)
1689  {
1690  ModelPart::NodeIterator NodeBegin;
1691  ModelPart::NodeIterator NodeEnd;
1692  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1693  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1694  {
1695  const array_1d<double, 3> &Vel = itNode->FastGetSolutionStepValue(VELOCITY);
1696 
1697  double NormVelNode = 0;
1698 
1699  for (unsigned int d = 0; d < 3; ++d)
1700  {
1701  NormVelNode += Vel[d] * Vel[d];
1702  NormV += Vel[d] * Vel[d];
1703  }
1704  }
1705  }
1706 
1707  BaseType::GetModelPart().GetCommunicator().GetDataCommunicator().SumAll(NormV);
1708 
1709  NormV = sqrt(NormV);
1710 
1711  if (NormV == 0.0)
1712  NormV = 1.00;
1713 
1714  return NormV;
1715  }
1716 
1717  double ComputePressureNorm()
1718  {
1719  ModelPart &rModelPart = BaseType::GetModelPart();
1720 
1721  double NormP = 0.00;
1722 
1723 #pragma omp parallel reduction(+ : NormP)
1724  {
1725  ModelPart::NodeIterator NodeBegin;
1726  ModelPart::NodeIterator NodeEnd;
1727  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1728  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1729  {
1730  const double Pr = itNode->FastGetSolutionStepValue(PRESSURE);
1731  NormP += Pr * Pr;
1732  }
1733  }
1734 
1735  BaseType::GetModelPart().GetCommunicator().GetDataCommunicator().SumAll(NormP);
1736 
1737  NormP = sqrt(NormP);
1738 
1739  if (NormP == 0.0)
1740  NormP = 1.00;
1741 
1742  return NormP;
1743  }
1744 
1745  void ComputeErrorL2NormCaseImposedG()
1746  {
1747 
1748  ModelPart &rModelPart = BaseType::GetModelPart();
1749  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1750  const double currentTime = rCurrentProcessInfo[TIME];
1751 
1752  double sumErrorL2Velocity = 0;
1753  double sumErrorL2VelocityX = 0;
1754  double sumErrorL2VelocityY = 0;
1755  double sumErrorL2Pressure = 0;
1756  double sumErrorL2TauXX = 0;
1757  double sumErrorL2TauYY = 0;
1758  double sumErrorL2TauXY = 0;
1759 
1760 #pragma omp parallel
1761  {
1762  ModelPart::NodeIterator NodeBegin;
1763  ModelPart::NodeIterator NodeEnd;
1764  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1765  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1766  {
1767  const double posX = itNode->X();
1768  const double posY = itNode->Y();
1769  const double nodalArea = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
1770  const double velX = itNode->FastGetSolutionStepValue(VELOCITY_X);
1771  const double velY = itNode->FastGetSolutionStepValue(VELOCITY_Y);
1772  const double pressure = itNode->FastGetSolutionStepValue(PRESSURE);
1773  const double tauXX = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS)[0];
1774  const double tauYY = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS)[1];
1775  const double tauXY = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS)[2];
1776 
1777  double expectedVelocityX = pow(posX, 2) * (1.0 - posX) * (1.0 - posX) * (2.0 * posY - 6.0 * pow(posY, 2) + 4.0 * pow(posY, 3));
1778  double expectedVelocityY = -pow(posY, 2) * (1.0 - posY) * (1.0 - posY) * (2.0 * posX - 6.0 * pow(posX, 2) + 4.0 * pow(posX, 3));
1779  double expectedPressure = -posX * (1.0 - posX);
1780 
1781  double expectedTauXX = 2.0 * (-4.0 * (1 - posX) * posX * (-1.0 + 2.0 * posX) * posY * (1.0 - 3.0 * posY + 2.0 * pow(posY, 2)));
1782  double expectedTauYY = 2.0 * (4.0 * posX * (1.0 - 3.0 * posX + 2.0 * pow(posX, 2)) * (1 - posY) * posY * (-1.0 + 2.0 * posY));
1783  double expectedTauXY = (2.0 * (1.0 - 6.0 * posY + 6.0 * pow(posY, 2)) * (1 - posX) * (1 - posX) * pow(posX, 2) - 2.0 * (1.0 - 6.0 * posX + 6.0 * pow(posX, 2)) * (1 - posY) * (1 - posY) * pow(posY, 2));
1784 
1785  double nodalErrorVelocityX = velX - expectedVelocityX;
1786  double nodalErrorVelocityY = velY - expectedVelocityY;
1787  double nodalErrorPressure = pressure - expectedPressure;
1788  double nodalErrorTauXX = tauXX - expectedTauXX;
1789  double nodalErrorTauYY = tauYY - expectedTauYY;
1790  double nodalErrorTauXY = tauXY - expectedTauXY;
1791 
1792  sumErrorL2Velocity += (pow(nodalErrorVelocityX, 2) + pow(nodalErrorVelocityY, 2)) * nodalArea;
1793  sumErrorL2VelocityX += pow(nodalErrorVelocityX, 2) * nodalArea;
1794  sumErrorL2VelocityY += pow(nodalErrorVelocityY, 2) * nodalArea;
1795  sumErrorL2Pressure += pow(nodalErrorPressure, 2) * nodalArea;
1796  sumErrorL2TauXX += pow(nodalErrorTauXX, 2) * nodalArea;
1797  sumErrorL2TauYY += pow(nodalErrorTauYY, 2) * nodalArea;
1798  sumErrorL2TauXY += pow(nodalErrorTauXY, 2) * nodalArea;
1799 
1800  // itNode->FastGetSolutionStepValue(NODAL_ERROR_XX)=nodalErrorTauXX;
1801  }
1802  }
1803 
1804  double errorL2Velocity = sqrt(sumErrorL2Velocity);
1805  double errorL2VelocityX = sqrt(sumErrorL2VelocityX);
1806  double errorL2VelocityY = sqrt(sumErrorL2VelocityY);
1807  double errorL2Pressure = sqrt(sumErrorL2Pressure);
1808  double errorL2TauXX = sqrt(sumErrorL2TauXX);
1809  double errorL2TauYY = sqrt(sumErrorL2TauYY);
1810  double errorL2TauXY = sqrt(sumErrorL2TauXY);
1811 
1812  std::ofstream myfileVelocity;
1813  myfileVelocity.open("errorL2VelocityFile.txt", std::ios::app);
1814  myfileVelocity << currentTime << "\t" << errorL2Velocity << "\n";
1815  myfileVelocity.close();
1816 
1817  std::ofstream myfileVelocityX;
1818  myfileVelocityX.open("errorL2VelocityXFile.txt", std::ios::app);
1819  myfileVelocityX << currentTime << "\t" << errorL2VelocityX << "\n";
1820  myfileVelocityX.close();
1821 
1822  std::ofstream myfileVelocityY;
1823  myfileVelocityY.open("errorL2VelocityYFile.txt", std::ios::app);
1824  myfileVelocityY << currentTime << "\t" << errorL2VelocityY << "\n";
1825  myfileVelocityY.close();
1826 
1827  std::ofstream myfilePressure;
1828  myfilePressure.open("errorL2PressureFile.txt", std::ios::app);
1829  myfilePressure << currentTime << "\t" << errorL2Pressure << "\n";
1830  myfilePressure.close();
1831 
1832  std::ofstream myfileTauXX;
1833  myfileTauXX.open("errorL2TauXXFile.txt", std::ios::app);
1834  myfileTauXX << currentTime << "\t" << errorL2TauXX << "\n";
1835  myfileTauXX.close();
1836 
1837  std::ofstream myfileTauYY;
1838  myfileTauYY.open("errorL2TauYYFile.txt", std::ios::app);
1839  myfileTauYY << currentTime << "\t" << errorL2TauYY << "\n";
1840  myfileTauYY.close();
1841 
1842  std::ofstream myfileTauXY;
1843  myfileTauXY.open("errorL2TauXYFile.txt", std::ios::app);
1844  myfileTauXY << currentTime << "\t" << errorL2TauXY << "\n";
1845  myfileTauXY.close();
1846  }
1847 
1848  void ComputeErrorL2NormCasePoiseuille()
1849  {
1850 
1851  ModelPart &rModelPart = BaseType::GetModelPart();
1852  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1853  const double currentTime = rCurrentProcessInfo[TIME];
1854 
1855  double sumErrorL2VelocityTheta = 0;
1856  double sumErrorL2TauTheta = 0;
1857 
1858  double r_in = 0.2;
1859  double R_out = 0.5;
1860  double kappa = r_in / R_out;
1861  double omega = 0.5;
1862  double viscosity = 100.0;
1863 
1864 #pragma omp parallel
1865  {
1866  ModelPart::NodeIterator NodeBegin;
1867  ModelPart::NodeIterator NodeEnd;
1868  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1869  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1870  {
1871  const double posX = itNode->X();
1872  const double posY = itNode->Y();
1873  const double rPos = sqrt(pow(posX, 2) + pow(posY, 2));
1874  const double cosalfa = posX / rPos;
1875  const double sinalfa = posY / rPos;
1876  const double sin2alfa = 2.0 * cosalfa * sinalfa;
1877  const double cos2alfa = 1.0 - 2.0 * pow(sinalfa, 2);
1878  const double nodalArea = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
1879  const double velX = itNode->FastGetSolutionStepValue(VELOCITY_X);
1880  const double velY = itNode->FastGetSolutionStepValue(VELOCITY_Y);
1881  const double tauXX = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS)[0];
1882  const double tauYY = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS)[1];
1883  const double tauXY = itNode->FastGetSolutionStepValue(NODAL_DEVIATORIC_CAUCHY_STRESS)[2];
1884 
1885  double expectedVelocityTheta = pow(kappa, 2) * omega * R_out / (1.0 - pow(kappa, 2)) * (R_out / rPos - rPos / R_out);
1886  double computedVelocityTheta = sqrt(pow(velX, 2) + pow(velY, 2));
1887  double nodalErrorVelocityTheta = computedVelocityTheta - expectedVelocityTheta;
1888 
1889  double expectedTauTheta = (2.0 * viscosity * pow(kappa, 2) * omega * pow(R_out, 2)) / (1.0 - pow(kappa, 2)) / pow(rPos, 2);
1890  double computedTauTheta = +(tauXX - tauYY) * sin2alfa / 2.0 - tauXY * cos2alfa;
1891  double nodalErrorTauTheta = computedTauTheta - expectedTauTheta;
1892  itNode->FastGetSolutionStepValue(NODAL_ERROR_XX) = computedVelocityTheta;
1893 
1894  sumErrorL2VelocityTheta += pow(nodalErrorVelocityTheta, 2) * nodalArea;
1895  sumErrorL2TauTheta += pow(nodalErrorTauTheta, 2) * nodalArea;
1896  }
1897  }
1898 
1899  double errorL2VelocityTheta = sqrt(sumErrorL2VelocityTheta);
1900  double errorL2TauTheta = sqrt(sumErrorL2TauTheta);
1901 
1902  std::ofstream myfileVelocity;
1903  myfileVelocity.open("errorL2Poiseuille.txt", std::ios::app);
1904  myfileVelocity << currentTime << "\t" << errorL2VelocityTheta << "\t" << errorL2TauTheta << "\n";
1905  myfileVelocity.close();
1906  }
1907 
1908  bool CheckPressureConvergence(const double NormDp, double &errorNormDp)
1909  {
1910  ModelPart &rModelPart = BaseType::GetModelPart();
1911 
1912  double NormP = 0.00;
1913  errorNormDp = 0;
1914 
1915  // #pragma omp parallel reduction(+:NormP)
1916  // {
1917  ModelPart::NodeIterator NodeBegin;
1918  ModelPart::NodeIterator NodeEnd;
1919  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1920  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1921  {
1922  const double Pr = itNode->FastGetSolutionStepValue(PRESSURE);
1923  NormP += Pr * Pr;
1924  }
1925  // }
1926 
1927  BaseType::GetModelPart().GetCommunicator().GetDataCommunicator().SumAll(NormP);
1928 
1929  NormP = sqrt(NormP);
1930 
1931  if (NormP == 0.0)
1932  NormP = 1.00;
1933 
1934  errorNormDp = NormDp / NormP;
1935 
1936  if (BaseType::GetEchoLevel() > 0 && rModelPart.GetCommunicator().MyPID() == 0)
1937  {
1938  std::cout << " The norm of pressure increment is: " << NormDp << std::endl;
1939  std::cout << " The norm of pressure is: " << NormP << std::endl;
1940  std::cout << " Pressure error: " << errorNormDp << std::endl;
1941  }
1942 
1943  if (errorNormDp < mPressureTolerance)
1944  {
1945  return true;
1946  }
1947  else
1948  return false;
1949  }
1950 
1951  bool FixTimeStepMomentum(const double DvErrorNorm)
1952  {
1953  ModelPart &rModelPart = BaseType::GetModelPart();
1954  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
1955  double currentTime = rCurrentProcessInfo[TIME];
1956  double timeInterval = rCurrentProcessInfo[DELTA_TIME];
1957  double minTolerance = 0.005;
1958  bool fixedTimeStep = false;
1959  if (currentTime < 3 * timeInterval)
1960  {
1961  minTolerance = 10;
1962  }
1963 
1964  if ((DvErrorNorm > minTolerance || (DvErrorNorm < 0 && DvErrorNorm > 0) || (DvErrorNorm != DvErrorNorm)) &&
1965  DvErrorNorm != 0 &&
1966  (DvErrorNorm != 1 || currentTime > timeInterval))
1967  {
1968  rCurrentProcessInfo.SetValue(BAD_VELOCITY_CONVERGENCE, true);
1969  std::cout << "NOT GOOD CONVERGENCE!!! I'll reduce the next time interval" << DvErrorNorm << std::endl;
1970  minTolerance = 0.05;
1971  if (DvErrorNorm > minTolerance)
1972  {
1973  std::cout << "BAD CONVERGENCE!!! I GO AHEAD WITH THE PREVIOUS VELOCITY AND PRESSURE FIELDS" << DvErrorNorm << std::endl;
1974  fixedTimeStep = true;
1975  // #pragma omp parallel
1976  // {
1977  ModelPart::NodeIterator NodeBegin;
1978  ModelPart::NodeIterator NodeEnd;
1979  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
1980  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
1981  {
1982  itNode->FastGetSolutionStepValue(VELOCITY, 0) = itNode->FastGetSolutionStepValue(VELOCITY, 1);
1983  itNode->FastGetSolutionStepValue(PRESSURE, 0) = itNode->FastGetSolutionStepValue(PRESSURE, 1);
1984  itNode->FastGetSolutionStepValue(ACCELERATION, 0) = itNode->FastGetSolutionStepValue(ACCELERATION, 1);
1985  }
1986  // }
1987  }
1988  }
1989  else
1990  {
1991  rCurrentProcessInfo.SetValue(BAD_VELOCITY_CONVERGENCE, false);
1992  }
1993  return fixedTimeStep;
1994  }
1995 
1996  bool CheckMomentumConvergence(const double DvErrorNorm)
1997  {
1998  ModelPart &rModelPart = BaseType::GetModelPart();
1999  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
2000  double currentTime = rCurrentProcessInfo[TIME];
2001  double timeInterval = rCurrentProcessInfo[DELTA_TIME];
2002  double minTolerance = 0.99999;
2003  bool fixedTimeStep = false;
2004 
2005  if ((DvErrorNorm > minTolerance || (DvErrorNorm < 0 && DvErrorNorm > 0) || (DvErrorNorm != DvErrorNorm)) &&
2006  DvErrorNorm != 0 &&
2007  (DvErrorNorm != 1 || currentTime > timeInterval))
2008  {
2009  rCurrentProcessInfo.SetValue(BAD_VELOCITY_CONVERGENCE, true);
2010  std::cout << " BAD CONVERGENCE DETECTED DURING THE ITERATIVE LOOP!!! error: " << DvErrorNorm << " higher than 0.9999" << std::endl;
2011  std::cout << " I GO AHEAD WITH THE PREVIOUS VELOCITY AND PRESSURE FIELDS" << std::endl;
2012  fixedTimeStep = true;
2013 #pragma omp parallel
2014  {
2015  ModelPart::NodeIterator NodeBegin;
2016  ModelPart::NodeIterator NodeEnd;
2017  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodeBegin, NodeEnd);
2018  for (ModelPart::NodeIterator itNode = NodeBegin; itNode != NodeEnd; ++itNode)
2019  {
2020  itNode->FastGetSolutionStepValue(VELOCITY, 0) = itNode->FastGetSolutionStepValue(VELOCITY, 1);
2021  itNode->FastGetSolutionStepValue(PRESSURE, 0) = itNode->FastGetSolutionStepValue(PRESSURE, 1);
2022  itNode->FastGetSolutionStepValue(ACCELERATION, 0) = itNode->FastGetSolutionStepValue(ACCELERATION, 1);
2023  }
2024  }
2025  }
2026  else
2027  {
2028  rCurrentProcessInfo.SetValue(BAD_VELOCITY_CONVERGENCE, false);
2029  }
2030  return fixedTimeStep;
2031  }
2032 
2033  bool FixTimeStepContinuity(const double DvErrorNorm)
2034  {
2035  ModelPart &rModelPart = BaseType::GetModelPart();
2036  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
2037  double currentTime = rCurrentProcessInfo[TIME];
2038  double timeInterval = rCurrentProcessInfo[DELTA_TIME];
2039  double minTolerance = 0.01;
2040  bool fixedTimeStep = false;
2041  if (currentTime < 3 * timeInterval)
2042  {
2043  minTolerance = 10;
2044  }
2045 
2046  if ((DvErrorNorm > minTolerance || (DvErrorNorm < 0 && DvErrorNorm > 0) || (DvErrorNorm != DvErrorNorm)) &&
2047  DvErrorNorm != 0 &&
2048  (DvErrorNorm != 1 || currentTime > timeInterval))
2049  {
2050  fixedTimeStep = true;
2051  rCurrentProcessInfo.SetValue(BAD_PRESSURE_CONVERGENCE, true);
2052  }
2053  else
2054  {
2055  rCurrentProcessInfo.SetValue(BAD_PRESSURE_CONVERGENCE, false);
2056  }
2057  return fixedTimeStep;
2058  }
2059 
2063 
2067 
2071 
2073 
2074  // private:
2077 
2081 
2082  double mVelocityTolerance;
2083 
2084  double mPressureTolerance;
2085 
2086  unsigned int mMaxPressureIter;
2087 
2088  unsigned int mDomainSize;
2089 
2090  unsigned int mTimeOrder;
2091 
2092  bool mReformDofSet;
2093 
2094  // Fractional step index.
2095  /* 1 : Momentum step (calculate fractional step velocity)
2096  * 2-3 : Unused (reserved for componentwise calculation of frac step velocity)
2097  * 4 : Pressure step
2098  * 5 : Computation of projections
2099  * 6 : End of step velocity
2100  */
2101  // unsigned int mStepId;
2102 
2104  StrategyPointerType mpMomentumStrategy;
2105 
2107  StrategyPointerType mpPressureStrategy;
2108 
2112 
2116 
2117  void InitializeStrategy(SolverSettingsType &rSolverConfig)
2118  {
2119  KRATOS_TRY;
2120 
2121  // Check that input parameters are reasonable and sufficient.
2122  this->Check();
2123 
2124  // ModelPart& rModelPart = this->GetModelPart();
2125 
2126  mDomainSize = rSolverConfig.GetDomainSize();
2127 
2128  mReformDofSet = rSolverConfig.GetReformDofSet();
2129 
2130  BaseType::SetEchoLevel(rSolverConfig.GetEchoLevel());
2131 
2132  // Initialize strategies for each step
2133  bool HaveVelStrategy = rSolverConfig.FindStrategy(SolverSettingsType::Velocity, mpMomentumStrategy);
2134 
2135  if (HaveVelStrategy)
2136  {
2137  rSolverConfig.FindTolerance(SolverSettingsType::Velocity, mVelocityTolerance);
2138  }
2139  else
2140  {
2141  KRATOS_THROW_ERROR(std::runtime_error, "NodalTwoStepVPStrategy error: No Velocity strategy defined in FractionalStepSettings", "");
2142  }
2143 
2144  bool HavePressStrategy = rSolverConfig.FindStrategy(SolverSettingsType::Pressure, mpPressureStrategy);
2145 
2146  if (HavePressStrategy)
2147  {
2148  rSolverConfig.FindTolerance(SolverSettingsType::Pressure, mPressureTolerance);
2149  rSolverConfig.FindMaxIter(SolverSettingsType::Pressure, mMaxPressureIter);
2150  }
2151  else
2152  {
2153  KRATOS_THROW_ERROR(std::runtime_error, "NodalTwoStepVPStrategy error: No Pressure strategy defined in FractionalStepSettings", "");
2154  }
2155 
2156  // Check input parameters
2157  this->Check();
2158 
2159  KRATOS_CATCH("");
2160  }
2161 
2162  private:
2163  void FillNodalSFDVector()
2164  {
2165 
2166  ModelPart &rModelPart = BaseType::GetModelPart();
2167 
2168  // #pragma omp parallel
2169  // {
2170  // ModelPart::NodeIterator NodesBegin;
2171  // ModelPart::NodeIterator NodesEnd;
2172  // OpenMPUtils::PartitionedIterators(rModelPart.Nodes(),NodesBegin,NodesEnd);
2173 
2174  // for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
2175  // {
2176 
2177  for (ModelPart::NodeIterator itNode = rModelPart.NodesBegin(); itNode != rModelPart.NodesEnd(); itNode++)
2178  {
2179  InitializeNodalVariablesForRemeshedDomain(itNode);
2180 
2181  SetNeighboursOrderToNode(itNode); // it assigns neighbours to inner nodes, filling NODAL_SFD_NEIGHBOURS_ORDER
2182  }
2183  }
2184 
2188 
2192 
2196 
2198  NodalTwoStepVPStrategy &operator=(NodalTwoStepVPStrategy const &rOther) {}
2199 
2201  NodalTwoStepVPStrategy(NodalTwoStepVPStrategy const &rOther) {}
2202 
2204 
2205  };
2206 
2210 
2212 
2214 
2215 } // namespace Kratos.
2216 
2217 #endif // KRATOS_NODAL_TWO_STEP_V_P_STRATEGY_H
boundary_normals_calculation_utilities.hpp
cfd_variables.h
Kratos::BoundaryNormalsCalculationUtilities
Definition: boundary_normals_calculation_utilities.hpp:48
Kratos::BoundaryNormalsCalculationUtilities::CalculateUnitBoundaryNormals
void CalculateUnitBoundaryNormals(ModelPart &rModelPart, int EchoLevel=0)
Calculates the area normal (vector oriented as the normal with a dimension proportional to the area).
Definition: boundary_normals_calculation_utilities.hpp:83
Kratos::BuilderAndSolver
Current class provides an implementation for the base builder and solving operations.
Definition: builder_and_solver.h:64
Kratos::Communicator::GetDataCommunicator
virtual const DataCommunicator & GetDataCommunicator() const
Definition: communicator.cpp:340
Kratos::Communicator::MyPID
virtual int MyPID() const
Definition: communicator.cpp:91
Kratos::DataValueContainer::SetValue
void SetValue(const Variable< TDataType > &rThisVariable, TDataType const &rValue)
Sets the value for a given variable.
Definition: data_value_container.h:320
Kratos::GaussSeidelLinearStrategy
Short class definition.
Definition: gauss_seidel_linear_strategy.h:83
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::size
SizeType size() const
Definition: geometry.h:518
Kratos::Geometry::Volume
virtual double Volume() const
This method calculate and return volume of this geometry.
Definition: geometry.h:1358
Kratos::Geometry::Area
virtual double Area() const
This method calculate and return area or surface area of this geometry depending to it's dimension.
Definition: geometry.h:1345
Kratos::GlobalPointersVector
This class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
Kratos::GlobalPointersVector::size
size_type size() const
Definition: global_pointers_vector.h:307
Kratos::ImplicitSolvingStrategy
Implicit solving strategy base class This is the base class from which we will derive all the implici...
Definition: implicit_solving_strategy.h:61
Kratos::ImplicitSolvingStrategy::ElementsArrayType
BaseType::ElementsArrayType ElementsArrayType
Definition: implicit_solving_strategy.h:94
Kratos::Internals::Matrix
Definition: amatrix_interface.h:41
Kratos::Internals::Matrix::resize
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
Kratos::MathUtils::InvertMatrix3
static void InvertMatrix3(const TMatrix1 &rInputMatrix, TMatrix2 &rInvertedMatrix, double &rInputMatrixDet)
It inverts matrices of order 3.
Definition: math_utils.h:449
Kratos::MathUtils::InvertMatrix2
static void InvertMatrix2(const TMatrix1 &rInputMatrix, TMatrix2 &rInvertedMatrix, double &rInputMatrixDet)
It inverts matrices of order 2.
Definition: math_utils.h:416
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::GetCommunicator
Communicator & GetCommunicator()
Definition: model_part.h:1821
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::NodeIterator
MeshType::NodeIterator NodeIterator
Definition: model_part.h:134
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::ModelPart::GetNodalSolutionStepVariablesList
VariablesList & GetNodalSolutionStepVariablesList()
Definition: model_part.h:549
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity
Current class provides an implementation for standard builder and solving operations.
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:80
Kratos::NodalResidualBasedEliminationBuilderAndSolver
Current class provides an implementation for standard builder and solving operations.
Definition: nodal_residualbased_elimination_builder_and_solver.h:80
Kratos::NodalTwoStepVPStrategy
Definition: nodal_two_step_v_p_strategy.h:70
Kratos::NodalTwoStepVPStrategy::InitializeNonLinearIterations
void InitializeNonLinearIterations()
Definition: nodal_two_step_v_p_strategy.h:659
Kratos::NodalTwoStepVPStrategy::SetTimeCoefficients
void SetTimeCoefficients(ProcessInfo &rCurrentProcessInfo)
Definition: nodal_two_step_v_p_strategy.h:208
Kratos::NodalTwoStepVPStrategy::ComputeErrorL2NormCasePoiseuille
void ComputeErrorL2NormCasePoiseuille()
Definition: nodal_two_step_v_p_strategy.h:1848
Kratos::NodalTwoStepVPStrategy::ComputeNodalVolume
void ComputeNodalVolume()
Definition: nodal_two_step_v_p_strategy.h:511
Kratos::NodalTwoStepVPStrategy::mPressureTolerance
double mPressureTolerance
Definition: nodal_two_step_v_p_strategy.h:2084
Kratos::NodalTwoStepVPStrategy::ComputeAndStoreNodalDeformationGradient
void ComputeAndStoreNodalDeformationGradient(ModelPart::NodeIterator itNode, double theta)
Definition: nodal_two_step_v_p_strategy.h:1062
Kratos::NodalTwoStepVPStrategy::CalculatePressureVelocity
void CalculatePressureVelocity()
Definition: nodal_two_step_v_p_strategy.h:1222
Kratos::NodalTwoStepVPStrategy::CalculatePressureAcceleration
void CalculatePressureAcceleration()
Definition: nodal_two_step_v_p_strategy.h:1247
Kratos::NodalTwoStepVPStrategy::CalculateAccelerations
void CalculateAccelerations()
Definition: nodal_two_step_v_p_strategy.h:1271
Kratos::NodalTwoStepVPStrategy::CalcNodalStrainsForNode
void CalcNodalStrainsForNode(ModelPart::NodeIterator itNode)
Definition: nodal_two_step_v_p_strategy.h:957
Kratos::NodalTwoStepVPStrategy::TSystemVectorType
BaseType::TSystemVectorType TSystemVectorType
Definition: nodal_two_step_v_p_strategy.h:93
Kratos::NodalTwoStepVPStrategy::NodalTwoStepVPStrategy
NodalTwoStepVPStrategy(ModelPart &rModelPart, typename TLinearSolver::Pointer pVelocityLinearSolver, typename TLinearSolver::Pointer pPressureLinearSolver, bool ReformDofSet=true, double VelTol=0.0001, double PresTol=0.0001, int MaxPressureIterations=1, unsigned int TimeOrder=2, unsigned int DomainSize=2)
Definition: nodal_two_step_v_p_strategy.h:118
Kratos::NodalTwoStepVPStrategy::CheckMomentumConvergence
bool CheckMomentumConvergence(const double DvErrorNorm)
Definition: nodal_two_step_v_p_strategy.h:1996
Kratos::NodalTwoStepVPStrategy::LocalSystemMatrixType
BaseType::LocalSystemMatrixType LocalSystemMatrixType
Definition: nodal_two_step_v_p_strategy.h:97
Kratos::NodalTwoStepVPStrategy::SolverSettingsType
TwoStepVPSolverSettings< TSparseSpace, TDenseSpace, TLinearSolver > SolverSettingsType
Definition: nodal_two_step_v_p_strategy.h:103
Kratos::NodalTwoStepVPStrategy::CheckVelocityConvergence
bool CheckVelocityConvergence(const double NormDv, double &errorNormDv)
Definition: nodal_two_step_v_p_strategy.h:1631
Kratos::NodalTwoStepVPStrategy::Clear
void Clear() override
Clears the internal storage.
Definition: nodal_two_step_v_p_strategy.h:1410
Kratos::NodalTwoStepVPStrategy::Initialize
void Initialize() override
Initialization of member variables and prior operations.
Definition: nodal_two_step_v_p_strategy.h:356
Kratos::NodalTwoStepVPStrategy::FinalizeSolutionStep
void FinalizeSolutionStep() override
Performs all the required operations that should be done (for each step) after solving the solution s...
Definition: nodal_two_step_v_p_strategy.h:351
Kratos::NodalTwoStepVPStrategy::ComputeVelocityNorm
double ComputeVelocityNorm()
Definition: nodal_two_step_v_p_strategy.h:1682
Kratos::NodalTwoStepVPStrategy::mpMomentumStrategy
StrategyPointerType mpMomentumStrategy
Scheme for the solution of the momentum equation.
Definition: nodal_two_step_v_p_strategy.h:2104
Kratos::NodalTwoStepVPStrategy::PrintInfo
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: nodal_two_step_v_p_strategy.h:1445
Kratos::NodalTwoStepVPStrategy::NodalTwoStepVPStrategy
NodalTwoStepVPStrategy(ModelPart &rModelPart, SolverSettingsType &rSolverConfig)
Definition: nodal_two_step_v_p_strategy.h:112
Kratos::NodalTwoStepVPStrategy::mMaxPressureIter
unsigned int mMaxPressureIter
Definition: nodal_two_step_v_p_strategy.h:2086
Kratos::NodalTwoStepVPStrategy::CalcNodalStrainsAndStresses
void CalcNodalStrainsAndStresses()
Definition: nodal_two_step_v_p_strategy.h:688
Kratos::NodalTwoStepVPStrategy::LocalSystemVectorType
BaseType::LocalSystemVectorType LocalSystemVectorType
Definition: nodal_two_step_v_p_strategy.h:95
Kratos::NodalTwoStepVPStrategy::BaseType
ImplicitSolvingStrategy< TSparseSpace, TDenseSpace, TLinearSolver > BaseType
Counted pointer of NodalTwoStepVPStrategy.
Definition: nodal_two_step_v_p_strategy.h:77
Kratos::NodalTwoStepVPStrategy::UpdateTopology
void UpdateTopology(ModelPart &rModelPart, unsigned int echoLevel)
Definition: nodal_two_step_v_p_strategy.h:1210
Kratos::NodalTwoStepVPStrategy::SolveSolutionStep
bool SolveSolutionStep() override
Solves the current step. This function returns true if a solution has been found, false otherwise.
Definition: nodal_two_step_v_p_strategy.h:243
Kratos::NodalTwoStepVPStrategy::SetNeighboursOrderToNode
void SetNeighboursOrderToNode(ModelPart::NodeIterator itNode)
Definition: nodal_two_step_v_p_strategy.h:556
Kratos::NodalTwoStepVPStrategy::SolveMomentumIteration
bool SolveMomentumIteration(unsigned int it, unsigned int maxIt, bool &fixedTimeStep, double &velocityNorm)
Calculate the coefficients for time iteration.
Definition: nodal_two_step_v_p_strategy.h:1486
Kratos::NodalTwoStepVPStrategy::InitializeStrategy
void InitializeStrategy(SolverSettingsType &rSolverConfig)
Definition: nodal_two_step_v_p_strategy.h:2117
Kratos::NodalTwoStepVPStrategy::UpdatePressureAccelerations
void UpdatePressureAccelerations()
Definition: nodal_two_step_v_p_strategy.h:1403
Kratos::NodalTwoStepVPStrategy::FixTimeStepContinuity
bool FixTimeStepContinuity(const double DvErrorNorm)
Definition: nodal_two_step_v_p_strategy.h:2033
Kratos::NodalTwoStepVPStrategy::FixTimeStepMomentum
bool FixTimeStepMomentum(const double DvErrorNorm)
Definition: nodal_two_step_v_p_strategy.h:1951
Kratos::NodalTwoStepVPStrategy::UpdateAccelerations
void UpdateAccelerations(array_1d< double, 3 > &CurrentAcceleration, const array_1d< double, 3 > &CurrentVelocity, array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &PreviousVelocity)
Definition: nodal_two_step_v_p_strategy.h:1318
Kratos::NodalTwoStepVPStrategy::StrategyPointerType
ImplicitSolvingStrategy< TSparseSpace, TDenseSpace, TLinearSolver >::Pointer StrategyPointerType
Definition: nodal_two_step_v_p_strategy.h:101
Kratos::NodalTwoStepVPStrategy::ComputeErrorL2NormCaseImposedG
void ComputeErrorL2NormCaseImposedG()
Definition: nodal_two_step_v_p_strategy.h:1745
Kratos::NodalTwoStepVPStrategy::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(NodalTwoStepVPStrategy)
Kratos::NodalTwoStepVPStrategy::GeometryType
Geometry< NodeType > GeometryType
Geometry type (using with given NodeType)
Definition: nodal_two_step_v_p_strategy.h:85
Kratos::NodalTwoStepVPStrategy::Info
std::string Info() const override
Turn back information as a string.
Definition: nodal_two_step_v_p_strategy.h:1437
Kratos::NodalTwoStepVPStrategy::CalcNodalStrainsAndStressesForNode
void CalcNodalStrainsAndStressesForNode(ModelPart::NodeIterator itNode)
Definition: nodal_two_step_v_p_strategy.h:718
Kratos::NodalTwoStepVPStrategy::mTimeOrder
unsigned int mTimeOrder
Definition: nodal_two_step_v_p_strategy.h:2090
Kratos::NodalTwoStepVPStrategy::InitializeNodalVariablesForRemeshedDomain
void InitializeNodalVariablesForRemeshedDomain(ModelPart::NodeIterator itNode)
Definition: nodal_two_step_v_p_strategy.h:578
Kratos::NodalTwoStepVPStrategy::InitialAssignFluidMaterialToEachNode
void InitialAssignFluidMaterialToEachNode(ModelPart::NodeIterator itNode)
Definition: nodal_two_step_v_p_strategy.h:492
Kratos::NodalTwoStepVPStrategy::SizeType
std::size_t SizeType
Definition: nodal_two_step_v_p_strategy.h:87
Kratos::NodalTwoStepVPStrategy::Check
int Check() override
Function to perform expensive checks.
Definition: nodal_two_step_v_p_strategy.h:179
Kratos::NodalTwoStepVPStrategy::ShapeFunctionDerivativesArrayType
GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType
Definition: nodal_two_step_v_p_strategy.h:105
Kratos::NodalTwoStepVPStrategy::ComputeDeviatoricCoefficientForFluid
void ComputeDeviatoricCoefficientForFluid(ModelPart::NodeIterator itNode, const unsigned int dimension, double &deviatoricCoefficient)
Definition: nodal_two_step_v_p_strategy.h:822
Kratos::NodalTwoStepVPStrategy::CheckPressureConvergence
bool CheckPressureConvergence(const double NormDp, double &errorNormDp)
Definition: nodal_two_step_v_p_strategy.h:1908
Kratos::NodalTwoStepVPStrategy::TDataType
BaseType::TDataType TDataType
Definition: nodal_two_step_v_p_strategy.h:79
Kratos::NodalTwoStepVPStrategy::TSystemMatrixType
BaseType::TSystemMatrixType TSystemMatrixType
Definition: nodal_two_step_v_p_strategy.h:91
Kratos::NodalTwoStepVPStrategy::mpPressureStrategy
StrategyPointerType mpPressureStrategy
Scheme for the solution of the mass equation.
Definition: nodal_two_step_v_p_strategy.h:2107
Kratos::NodalTwoStepVPStrategy::DofsArrayType
BaseType::DofsArrayType DofsArrayType
Definition: nodal_two_step_v_p_strategy.h:89
Kratos::NodalTwoStepVPStrategy::CalcNodalStrains
void CalcNodalStrains()
Definition: nodal_two_step_v_p_strategy.h:1031
Kratos::NodalTwoStepVPStrategy::mDomainSize
unsigned int mDomainSize
Definition: nodal_two_step_v_p_strategy.h:2088
Kratos::NodalTwoStepVPStrategy::ElementsArrayType
BaseType::ElementsArrayType ElementsArrayType
Definition: nodal_two_step_v_p_strategy.h:99
Kratos::NodalTwoStepVPStrategy::CalculateDisplacements
void CalculateDisplacements()
Definition: nodal_two_step_v_p_strategy.h:1329
Kratos::NodalTwoStepVPStrategy::~NodalTwoStepVPStrategy
virtual ~NodalTwoStepVPStrategy()
Destructor.
Definition: nodal_two_step_v_p_strategy.h:177
Kratos::NodalTwoStepVPStrategy::SolveContinuityIteration
bool SolveContinuityIteration(unsigned int it, unsigned int maxIt, double &NormP)
Definition: nodal_two_step_v_p_strategy.h:1559
Kratos::NodalTwoStepVPStrategy::ComputePressureNorm
double ComputePressureNorm()
Definition: nodal_two_step_v_p_strategy.h:1717
Kratos::NodalTwoStepVPStrategy::mVelocityTolerance
double mVelocityTolerance
Definition: nodal_two_step_v_p_strategy.h:2082
Kratos::NodalTwoStepVPStrategy::mReformDofSet
bool mReformDofSet
Definition: nodal_two_step_v_p_strategy.h:2092
Kratos::NodalTwoStepVPStrategy::NodeType
Node NodeType
Node type (default is: Node)
Definition: nodal_two_step_v_p_strategy.h:82
Kratos::NodalTwoStepVPStrategy::CalculateDisplacementsAndResetNodalVariables
void CalculateDisplacementsAndResetNodalVariables()
Definition: nodal_two_step_v_p_strategy.h:1350
Kratos::NodalTwoStepVPStrategy::PrintData
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: nodal_two_step_v_p_strategy.h:1451
Kratos::NodalTwoStepVPStrategy::NodeWeakPtrVectorType
GlobalPointersVector< Node > NodeWeakPtrVectorType
Definition: nodal_two_step_v_p_strategy.h:107
Kratos::NodalTwoStepVPStrategy::SetEchoLevel
void SetEchoLevel(int Level) override
This sets the level of echo for the solving strategy.
Definition: nodal_two_step_v_p_strategy.h:1420
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::OpenMPUtils::PartitionedIterators
static void PartitionedIterators(TVector &rVector, typename TVector::iterator &rBegin, typename TVector::iterator &rEnd)
Generate a partition for an std::vector-like array, providing iterators to the begin and end position...
Definition: openmp_utils.h:179
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::ProcessInfo::GetPreviousTimeStepInfo
ProcessInfo & GetPreviousTimeStepInfo(IndexType StepsBefore=1)
Definition: process_info.h:187
Kratos::ResidualBasedIncrementalUpdateStaticScheme
This class provides the implementation of a static scheme.
Definition: residualbased_incrementalupdate_static_scheme.h:57
Kratos::Scheme
This class provides the implementation of the basic tasks that are needed by the solution strategy.
Definition: scheme.h:56
Kratos::SolvingStrategy
Solving strategy base class This is the base class from which we will derive all the strategies (impl...
Definition: solving_strategy.h:64
Kratos::SolvingStrategy::LocalSystemMatrixType
TDenseSpace::MatrixType LocalSystemMatrixType
Definition: solving_strategy.h:79
Kratos::SolvingStrategy::TDataType
TSparseSpace::DataType TDataType
Definition: solving_strategy.h:69
Kratos::SolvingStrategy::SetEchoLevel
virtual void SetEchoLevel(const int Level)
This sets the level of echo for the solving strategy.
Definition: solving_strategy.h:255
Kratos::SolvingStrategy::GetModelPart
ModelPart & GetModelPart()
Operations to get the pointer to the model.
Definition: solving_strategy.h:350
Kratos::SolvingStrategy::ElementsArrayType
ModelPart::ElementsContainerType ElementsArrayType
Definition: solving_strategy.h:91
Kratos::SolvingStrategy::GetEchoLevel
int GetEchoLevel()
This returns the level of echo for the solving strategy.
Definition: solving_strategy.h:271
Kratos::SolvingStrategy::TSystemMatrixType
TSparseSpace::MatrixType TSystemMatrixType
Definition: solving_strategy.h:71
Kratos::SolvingStrategy::Check
virtual int Check()
Function to perform expensive checks.
Definition: solving_strategy.h:377
Kratos::SolvingStrategy::TSystemVectorType
TSparseSpace::VectorType TSystemVectorType
Definition: solving_strategy.h:73
Kratos::SolvingStrategy::MoveMesh
virtual void MoveMesh()
This function is designed to move the mesh.
Definition: solving_strategy.h:330
Kratos::SolvingStrategy::LocalSystemVectorType
TDenseSpace::VectorType LocalSystemVectorType
Definition: solving_strategy.h:81
Kratos::TwoStepVPSolverSettings
Helper class to define solution strategies for TwoStepVPStrategy.
Definition: solver_settings.h:57
Kratos::TwoStepVPSolverSettings::FindTolerance
virtual bool FindTolerance(StrategyLabel const &rStrategyLabel, double &rTolerance)
Definition: solver_settings.h:143
Kratos::TwoStepVPSolverSettings::Velocity
@ Velocity
Definition: solver_settings.h:71
Kratos::TwoStepVPSolverSettings::Pressure
@ Pressure
Definition: solver_settings.h:72
Kratos::TwoStepVPSolverSettings::FindMaxIter
virtual bool FindMaxIter(StrategyLabel const &rStrategyLabel, unsigned int &rMaxIter)
Definition: solver_settings.h:156
Kratos::TwoStepVPSolverSettings::GetDomainSize
virtual unsigned int GetDomainSize() const
Definition: solver_settings.h:124
Kratos::TwoStepVPSolverSettings::GetEchoLevel
virtual unsigned int GetEchoLevel()
Definition: solver_settings.h:176
Kratos::TwoStepVPSolverSettings::FindStrategy
virtual bool FindStrategy(StrategyLabel const &rStrategyLabel, StrategyPointerType &pThisStrategy)
Definition: solver_settings.h:129
Kratos::TwoStepVPSolverSettings::GetReformDofSet
bool GetReformDofSet()
Definition: solver_settings.h:181
Kratos::VariablesList::Has
bool Has(const VariableData &rThisVariable) const
Definition: variables_list.h:372
Kratos::array_1d< double, 3 >
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_TRY
#define KRATOS_TRY
Definition: define.h:109
deprecated_variables.h
gauss_seidel_linear_strategy.h
KRATOS_INFO
#define KRATOS_INFO(label)
Definition: logger.h:250
implicit_solving_strategy.h
mesher_utilities.hpp
model_part.h
Kratos::Globals::Pi
constexpr double Pi
Definition: global_variables.h:25
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
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
PecletTest.Temp
Temp
Definition: PecletTest.py:105
edgebased_var.viscosity
float viscosity
Definition: edgebased_var.py:8
face_heat.Dt
Dt
Definition: face_heat.py:78
isotropic_damage_automatic_differentiation.dimension
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
ode_solve.kappa
kappa
Definition: ode_solve.py:416
ode_solve.d
int d
Definition: ode_solve.py:397
quadrature.k
int k
Definition: quadrature.py:595
tensors::i
integer i
Definition: TensorModule.f:17
test_pureconvectionsolver_benchmarking.pressure
float pressure
Definition: test_pureconvectionsolver_benchmarking.py:101
test_pureconvectionsolver_benchmarking.ReformDofAtEachIteration
ReformDofAtEachIteration
Definition: test_pureconvectionsolver_benchmarking.py:131
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
nodal_residualbased_elimination_builder_and_solver.h
nodal_residualbased_elimination_builder_and_solver_continuity.h
openmp_utils.h
pfem_fluid_dynamics_application_variables.h
residualbased_incrementalupdate_static_scheme.h
Generated on Sat Feb 3 2024 02:09:36 for KratosMultiphysics by doxygen 1.9.1