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.
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: January 2016 $
5 // Revision: $Revision: 0.0 $
6 //
7 //
8 
9 #ifndef KRATOS_V_P_STRATEGY_H
10 #define KRATOS_V_P_STRATEGY_H
11 
12 #include "includes/define.h"
13 #include "includes/model_part.h"
14 #include "includes/deprecated_variables.h"
15 #include "utilities/openmp_utils.h"
16 #include "processes/process.h"
17 #include "solving_strategies/schemes/scheme.h"
18 #include "solving_strategies/strategies/solving_strategy.h"
19 #include "custom_utilities/mesher_utilities.hpp"
20 #include "custom_utilities/boundary_normals_calculation_utilities.hpp"
21 
22 #include "custom_utilities/solver_settings.h"
23 
24 #include "pfem_fluid_dynamics_application_variables.h"
25 
26 #include <stdio.h>
27 #include <cmath>
28 
29 namespace Kratos
30 {
31 
34 
37 
41 
43 
46 
50 
54 
55  template <class TSparseSpace,
56  class TDenseSpace,
57  class TLinearSolver>
58  class VPStrategy : public SolvingStrategy<TSparseSpace, TDenseSpace>
59  {
60  public:
63  KRATOS_CLASS_POINTER_DEFINITION(VPStrategy);
64 
65  typedef SolvingStrategy<TSparseSpace, TDenseSpace> BaseType;
66 
67  typedef TwoStepVPSolverSettings<TSparseSpace, TDenseSpace, TLinearSolver> SolverSettingsType;
68 
72 
73  VPStrategy(ModelPart &rModelPart,
74  SolverSettingsType &rSolverConfig) : BaseType(rModelPart)
75  {
76  std::cout << "VPStrategy INITIALIZE STRATEGY" << std::endl;
77  InitializeStrategy(rSolverConfig);
78  }
79 
80  VPStrategy(ModelPart &rModelPart,
81  typename TLinearSolver::Pointer pVelocityLinearSolver,
82  typename TLinearSolver::Pointer pPressureLinearSolver,
83  bool ReformDofSet = true,
84  unsigned int DomainSize = 2) : BaseType(rModelPart)
85  {
86  KRATOS_TRY;
87  KRATOS_CATCH("");
88  }
89 
91  virtual ~VPStrategy() {}
92 
93  virtual int Check() override
94  {
95  return false;
96  }
97 
98  virtual bool SolveSolutionStep() override
99  {
100  return false;
101  }
102 
103  virtual void FinalizeSolutionStep() override {}
104 
105  virtual void InitializeSolutionStep() override {}
106 
107  void UpdateTopology(ModelPart &rModelPart, unsigned int echoLevel)
108  {
109  KRATOS_TRY;
110 
111  this->CalculateDisplacementsAndPorosity();
112  BaseType::MoveMesh();
113  KRATOS_CATCH("");
114  }
115 
116  void SetBlockedAndIsolatedFlags()
117  {
118  KRATOS_TRY;
119 
120  ModelPart &rModelPart = BaseType::GetModelPart();
121  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
122 
123 #pragma omp parallel
124  {
125  ModelPart::ElementIterator ElemBegin;
126  ModelPart::ElementIterator ElemEnd;
127  OpenMPUtils::PartitionedIterators(rModelPart.Elements(), ElemBegin, ElemEnd);
128  for (ModelPart::ElementIterator itElem = ElemBegin; itElem != ElemEnd; ++itElem)
129  {
130  unsigned int numNodes = itElem->GetGeometry().size();
131  std::vector<array_1d<double, 3>> nodesCoordinates;
132  nodesCoordinates.resize(numNodes);
133 
134  (itElem)->Set(BLOCKED, false);
135  (itElem)->Set(ISOLATED, false);
136 
137  unsigned int freeSurfaceNodes = 0;
138  unsigned int freeSurfaceRigidNodes = 0;
139  unsigned int rigidNodes = 0;
140  unsigned int isolatedNodes = 0;
141  for (unsigned int i = 0; i < numNodes; i++)
142  {
143  if (itElem->GetGeometry()[i].Is(FREE_SURFACE))
144  {
145  freeSurfaceNodes++;
146  if (itElem->GetGeometry()[i].Is(RIGID))
147  {
148  freeSurfaceRigidNodes++;
149  }
150  }
151  else if (itElem->GetGeometry()[i].Is(RIGID))
152  {
153  rigidNodes++;
154  }
155  nodesCoordinates[i] = itElem->GetGeometry()[i].Coordinates();
156  ElementWeakPtrVectorType &neighb_elems = itElem->GetGeometry()[i].GetValue(NEIGHBOUR_ELEMENTS);
157  if (neighb_elems.size() == 1)
158  {
159  isolatedNodes++;
160  }
161  }
162 
163  if (dimension == 3)
164  {
165  double a1 = 0; // slope x for plane on the first triangular face of the tetrahedra (nodes A,B,C)
166  double b1 = 0; // slope y for plane on the first triangular face of the tetrahedra (nodes A,B,C)
167  double c1 = 0; // slope z for plane on the first triangular face of the tetrahedra (nodes A,B,C)
168  a1 = (nodesCoordinates[1][1] - nodesCoordinates[0][1]) * (nodesCoordinates[2][2] - nodesCoordinates[0][2]) - (nodesCoordinates[2][1] - nodesCoordinates[0][1]) * (nodesCoordinates[1][2] - nodesCoordinates[0][2]);
169  b1 = (nodesCoordinates[1][2] - nodesCoordinates[0][2]) * (nodesCoordinates[2][0] - nodesCoordinates[0][0]) - (nodesCoordinates[2][2] - nodesCoordinates[0][2]) * (nodesCoordinates[1][0] - nodesCoordinates[0][0]);
170  c1 = (nodesCoordinates[1][0] - nodesCoordinates[0][0]) * (nodesCoordinates[2][1] - nodesCoordinates[0][1]) - (nodesCoordinates[2][0] - nodesCoordinates[0][0]) * (nodesCoordinates[1][1] - nodesCoordinates[0][1]);
171  double a2 = 0; // slope x for plane on the second triangular face of the tetrahedra (nodes A,B,D)
172  double b2 = 0; // slope y for plane on the second triangular face of the tetrahedra (nodes A,B,D)
173  double c2 = 0; // slope z for plane on the second triangular face of the tetrahedra (nodes A,B,D)
174  a2 = (nodesCoordinates[1][1] - nodesCoordinates[0][1]) * (nodesCoordinates[3][2] - nodesCoordinates[0][2]) - (nodesCoordinates[3][1] - nodesCoordinates[0][1]) * (nodesCoordinates[1][2] - nodesCoordinates[0][2]);
175  b2 = (nodesCoordinates[1][2] - nodesCoordinates[0][2]) * (nodesCoordinates[3][0] - nodesCoordinates[0][0]) - (nodesCoordinates[3][2] - nodesCoordinates[0][2]) * (nodesCoordinates[1][0] - nodesCoordinates[0][0]);
176  c2 = (nodesCoordinates[1][0] - nodesCoordinates[0][0]) * (nodesCoordinates[3][1] - nodesCoordinates[0][1]) - (nodesCoordinates[3][0] - nodesCoordinates[0][0]) * (nodesCoordinates[1][1] - nodesCoordinates[0][1]);
177  double a3 = 0; // slope x for plane on the third triangular face of the tetrahedra (nodes B,C,D)
178  double b3 = 0; // slope y for plane on the third triangular face of the tetrahedra (nodes B,C,D)
179  double c3 = 0; // slope z for plane on the third triangular face of the tetrahedra (nodes B,C,D)
180  a3 = (nodesCoordinates[1][1] - nodesCoordinates[2][1]) * (nodesCoordinates[3][2] - nodesCoordinates[2][2]) - (nodesCoordinates[3][1] - nodesCoordinates[2][1]) * (nodesCoordinates[1][2] - nodesCoordinates[2][2]);
181  b3 = (nodesCoordinates[1][2] - nodesCoordinates[2][2]) * (nodesCoordinates[3][0] - nodesCoordinates[2][0]) - (nodesCoordinates[3][2] - nodesCoordinates[2][2]) * (nodesCoordinates[1][0] - nodesCoordinates[2][0]);
182  c3 = (nodesCoordinates[1][0] - nodesCoordinates[2][0]) * (nodesCoordinates[3][1] - nodesCoordinates[2][1]) - (nodesCoordinates[3][0] - nodesCoordinates[2][0]) * (nodesCoordinates[1][1] - nodesCoordinates[2][1]);
183  double a4 = 0; // slope x for plane on the fourth triangular face of the tetrahedra (nodes A,C,D)
184  double b4 = 0; // slope y for plane on the fourth triangular face of the tetrahedra (nodes A,C,D)
185  double c4 = 0; // slope z for plane on the fourth triangular face of the tetrahedra (nodes A,C,D)
186  a4 = (nodesCoordinates[0][1] - nodesCoordinates[2][1]) * (nodesCoordinates[3][2] - nodesCoordinates[2][2]) - (nodesCoordinates[3][1] - nodesCoordinates[2][1]) * (nodesCoordinates[0][2] - nodesCoordinates[2][2]);
187  b4 = (nodesCoordinates[0][2] - nodesCoordinates[2][2]) * (nodesCoordinates[3][0] - nodesCoordinates[2][0]) - (nodesCoordinates[3][2] - nodesCoordinates[2][2]) * (nodesCoordinates[0][0] - nodesCoordinates[2][0]);
188  c4 = (nodesCoordinates[0][0] - nodesCoordinates[2][0]) * (nodesCoordinates[3][1] - nodesCoordinates[2][1]) - (nodesCoordinates[3][0] - nodesCoordinates[2][0]) * (nodesCoordinates[0][1] - nodesCoordinates[2][1]);
189 
190  double cosAngle12 = (a1 * a2 + b1 * b2 + c1 * c2) / (sqrt(pow(a1, 2) + pow(b1, 2) + pow(c1, 2)) * sqrt(pow(a2, 2) + pow(b2, 2) + pow(c2, 2)));
191  double cosAngle13 = (a1 * a3 + b1 * b3 + c1 * c3) / (sqrt(pow(a1, 2) + pow(b1, 2) + pow(c1, 2)) * sqrt(pow(a3, 2) + pow(b3, 2) + pow(c3, 2)));
192  double cosAngle14 = (a1 * a4 + b1 * b4 + c1 * c4) / (sqrt(pow(a1, 2) + pow(b1, 2) + pow(c1, 2)) * sqrt(pow(a4, 2) + pow(b4, 2) + pow(c4, 2)));
193  double cosAngle23 = (a3 * a2 + b3 * b2 + c3 * c2) / (sqrt(pow(a3, 2) + pow(b3, 2) + pow(c3, 2)) * sqrt(pow(a2, 2) + pow(b2, 2) + pow(c2, 2)));
194  double cosAngle24 = (a4 * a2 + b4 * b2 + c4 * c2) / (sqrt(pow(a4, 2) + pow(b4, 2) + pow(c4, 2)) * sqrt(pow(a2, 2) + pow(b2, 2) + pow(c2, 2)));
195  double cosAngle34 = (a4 * a3 + b4 * b3 + c4 * c3) / (sqrt(pow(a4, 2) + pow(b4, 2) + pow(c4, 2)) * sqrt(pow(a3, 2) + pow(b3, 2) + pow(c3, 2)));
196 
197  if ((fabs(cosAngle12) > 0.99 || fabs(cosAngle13) > 0.99 || fabs(cosAngle14) > 0.99 || fabs(cosAngle23) > 0.99 || fabs(cosAngle24) > 0.99 || fabs(cosAngle34) > 0.99) && (freeSurfaceNodes == numNodes) && isolatedNodes > 1)
198  {
199  (itElem)->Set(BLOCKED, true);
200  // std::cout << "in the strategy BLOCKED ELEMENT: " << (itElem)->Id() << std::endl;
201  }
202  else if ((fabs(cosAngle12) > 0.995 || fabs(cosAngle13) > 0.995 || fabs(cosAngle14) > 0.995 || fabs(cosAngle23) > 0.995 || fabs(cosAngle24) > 0.995 || fabs(cosAngle34) > 0.995) && (freeSurfaceNodes == numNodes) && isolatedNodes == 1)
203  {
204  (itElem)->Set(BLOCKED, true);
205  // std::cout << "in the strategy BLOCKED ELEMENT: " << (itElem)->Id() << std::endl;
206  }
207  else if ((fabs(cosAngle12) > 0.999 || fabs(cosAngle13) > 0.999 || fabs(cosAngle14) > 0.999 || fabs(cosAngle23) > 0.999 || fabs(cosAngle24) > 0.999 || fabs(cosAngle34) > 0.999) && (freeSurfaceNodes == numNodes))
208  {
209  (itElem)->Set(BLOCKED, true);
210  // std::cout << "in the strategy BLOCKED ELEMENT: " << (itElem)->Id() << std::endl;
211  }
212  }
213 
214  if (freeSurfaceNodes == numNodes && rigidNodes == 0 && isolatedNodes >= (numNodes - 1))
215  {
216  (itElem)->Set(ISOLATED, true);
217  (itElem)->Set(BLOCKED, false);
218  }
219  }
220  }
221  KRATOS_CATCH("");
222  }
223 
224  void CalculatePressureVelocity()
225  {
226  ModelPart &rModelPart = BaseType::GetModelPart();
227  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
228  const double timeInterval = rCurrentProcessInfo[DELTA_TIME];
229  unsigned int timeStep = rCurrentProcessInfo[STEP];
230 
231  for (ModelPart::NodeIterator i = rModelPart.NodesBegin();
232  i != rModelPart.NodesEnd(); ++i)
233  {
234  if (timeStep == 1)
235  {
236  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) = 0;
237  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1) = 0;
238  }
239  else
240  {
241  double &CurrentPressure = (i)->FastGetSolutionStepValue(PRESSURE, 0);
242  double &PreviousPressure = (i)->FastGetSolutionStepValue(PRESSURE, 1);
243  double &CurrentPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0);
244  CurrentPressureVelocity = (CurrentPressure - PreviousPressure) / timeInterval;
245  }
246  }
247  }
248 
249  void CalculatePressureAcceleration()
250  {
251  ModelPart &rModelPart = BaseType::GetModelPart();
252  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
253  const double timeInterval = rCurrentProcessInfo[DELTA_TIME];
254  unsigned int timeStep = rCurrentProcessInfo[STEP];
255 
256  for (ModelPart::NodeIterator i = rModelPart.NodesBegin();
257  i != rModelPart.NodesEnd(); ++i)
258  {
259  if (timeStep == 1)
260  {
261  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0) = 0;
262  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 1) = 0;
263  }
264  else
265  {
266  double &CurrentPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0);
267  double &PreviousPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1);
268  double &CurrentPressureAcceleration = (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0);
269  CurrentPressureAcceleration = (CurrentPressureVelocity - PreviousPressureVelocity) / timeInterval;
270  }
271  }
272  }
273 
274  virtual void CalculateTemporalVariables()
275  {
276  ModelPart &rModelPart = BaseType::GetModelPart();
277  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
278 
279  for (ModelPart::NodeIterator i = rModelPart.NodesBegin();
280  i != rModelPart.NodesEnd(); ++i)
281  {
282 
283  array_1d<double, 3> &CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 0);
284  array_1d<double, 3> &PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 1);
285 
286  array_1d<double, 3> &CurrentAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION, 0);
287  array_1d<double, 3> &PreviousAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION, 1);
288 
289  /* if((i)->IsNot(ISOLATED) || (i)->Is(SOLID)){ */
290  if ((i)->IsNot(ISOLATED) && ((i)->IsNot(RIGID) || (i)->Is(SOLID)))
291  {
292  UpdateAccelerations(CurrentAcceleration, CurrentVelocity, PreviousAcceleration, PreviousVelocity);
293  }
294  else if ((i)->Is(RIGID))
295  {
296  array_1d<double, 3> Zeros(3, 0.0);
297  (i)->FastGetSolutionStepValue(ACCELERATION, 0) = Zeros;
298  (i)->FastGetSolutionStepValue(ACCELERATION, 1) = Zeros;
299  }
300  else
301  {
302  (i)->FastGetSolutionStepValue(PRESSURE, 0) = 0.0;
303  (i)->FastGetSolutionStepValue(PRESSURE, 1) = 0.0;
304  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) = 0.0;
305  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1) = 0.0;
306  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0) = 0.0;
307  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 1) = 0.0;
308  if ((i)->SolutionStepsDataHas(VOLUME_ACCELERATION))
309  {
310  array_1d<double, 3> &VolumeAcceleration = (i)->FastGetSolutionStepValue(VOLUME_ACCELERATION);
311  (i)->FastGetSolutionStepValue(ACCELERATION, 0) = VolumeAcceleration;
312  (i)->FastGetSolutionStepValue(VELOCITY, 0) += VolumeAcceleration * rCurrentProcessInfo[DELTA_TIME];
313  }
314  }
315 
316  const double timeInterval = rCurrentProcessInfo[DELTA_TIME];
317  unsigned int timeStep = rCurrentProcessInfo[STEP];
318  if (timeStep == 1)
319  {
320  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) = 0;
321  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1) = 0;
322  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0) = 0;
323  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 1) = 0;
324  }
325  else
326  {
327  double &CurrentPressure = (i)->FastGetSolutionStepValue(PRESSURE, 0);
328  double &PreviousPressure = (i)->FastGetSolutionStepValue(PRESSURE, 1);
329  double &CurrentPressureVelocity = (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0);
330  double &CurrentPressureAcceleration = (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0);
331 
332  CurrentPressureAcceleration = CurrentPressureVelocity / timeInterval;
333 
334  CurrentPressureVelocity = (CurrentPressure - PreviousPressure) / timeInterval;
335 
336  CurrentPressureAcceleration += -CurrentPressureVelocity / timeInterval;
337  }
338  }
339  }
340 
341  void CalculateAccelerations()
342  {
343  ModelPart &rModelPart = BaseType::GetModelPart();
344  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
345 
346  for (ModelPart::NodeIterator i = rModelPart.NodesBegin();
347  i != rModelPart.NodesEnd(); ++i)
348  {
349 
350  array_1d<double, 3> &CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 0);
351  array_1d<double, 3> &PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 1);
352 
353  array_1d<double, 3> &CurrentAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION, 0);
354  array_1d<double, 3> &PreviousAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION, 1);
355 
356  /* if((i)->IsNot(ISOLATED) || (i)->Is(SOLID)){ */
357  if ((i)->IsNot(ISOLATED) && ((i)->IsNot(RIGID) || (i)->Is(SOLID)))
358  {
359  UpdateAccelerations(CurrentAcceleration, CurrentVelocity, PreviousAcceleration, PreviousVelocity);
360  }
361  else if ((i)->Is(RIGID))
362  {
363  array_1d<double, 3> Zeros(3, 0.0);
364  (i)->FastGetSolutionStepValue(ACCELERATION, 0) = Zeros;
365  (i)->FastGetSolutionStepValue(ACCELERATION, 1) = Zeros;
366  }
367  else
368  {
369  (i)->FastGetSolutionStepValue(PRESSURE, 0) = 0.0;
370  (i)->FastGetSolutionStepValue(PRESSURE, 1) = 0.0;
371  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) = 0.0;
372  (i)->FastGetSolutionStepValue(PRESSURE_VELOCITY, 1) = 0.0;
373  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 0) = 0.0;
374  (i)->FastGetSolutionStepValue(PRESSURE_ACCELERATION, 1) = 0.0;
375  if ((i)->SolutionStepsDataHas(VOLUME_ACCELERATION))
376  {
377  array_1d<double, 3> &VolumeAcceleration = (i)->FastGetSolutionStepValue(VOLUME_ACCELERATION);
378  (i)->FastGetSolutionStepValue(ACCELERATION, 0) = VolumeAcceleration;
379  (i)->FastGetSolutionStepValue(VELOCITY, 0) += VolumeAcceleration * rCurrentProcessInfo[DELTA_TIME];
380  }
381  }
382  }
383  }
384 
385  inline void UpdateAccelerations(array_1d<double, 3> &CurrentAcceleration,
386  const array_1d<double, 3> &CurrentVelocity,
387  array_1d<double, 3> &PreviousAcceleration,
388  const array_1d<double, 3> &PreviousVelocity)
389  {
390  ModelPart &rModelPart = BaseType::GetModelPart();
391  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
392  double Dt = rCurrentProcessInfo[DELTA_TIME];
393  noalias(CurrentAcceleration) = 2.0 * (CurrentVelocity - PreviousVelocity) / Dt - PreviousAcceleration;
394  }
395 
396  virtual void CalculateDisplacementsAndPorosity()
397  {
398  ModelPart &rModelPart = BaseType::GetModelPart();
399  ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
400  const double TimeStep = rCurrentProcessInfo[DELTA_TIME];
401 
402  for (ModelPart::NodeIterator i = rModelPart.NodesBegin();
403  i != rModelPart.NodesEnd(); ++i)
404  {
405  if ((i)->IsNot(PFEMFlags::EULERIAN_INLET))
406 
407  {
408  array_1d<double, 3> &CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 0);
409  array_1d<double, 3> &PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY, 1);
410 
411  array_1d<double, 3> &CurrentDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT, 0);
412  array_1d<double, 3> &PreviousDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT, 1);
413 
414  /* if( i->IsFixed(DISPLACEMENT_X) == false ) */
415  CurrentDisplacement[0] = 0.5 * TimeStep * (CurrentVelocity[0] + PreviousVelocity[0]) + PreviousDisplacement[0];
416 
417  /* if( i->IsFixed(DISPLACEMENT_Y) == false ) */
418  CurrentDisplacement[1] = 0.5 * TimeStep * (CurrentVelocity[1] + PreviousVelocity[1]) + PreviousDisplacement[1];
419 
420  /* if( i->IsFixed(DISPLACEMENT_Z) == false ) */
421  CurrentDisplacement[2] = 0.5 * TimeStep * (CurrentVelocity[2] + PreviousVelocity[2]) + PreviousDisplacement[2];
422 
423  // currentFluidFractionRate = (currentFluidFraction - previousFluidFraction)/TimeStep;
424  }
425  }
426  }
427 
428  virtual void UpdateStressStrain() {}
429 
430  virtual void Clear() override {}
431 
435 
436  virtual void SetEchoLevel(int Level) override
437  {
438  BaseType::SetEchoLevel(Level);
439  }
440 
444 
448 
450  std::string Info() const override
451  {
452  std::stringstream buffer;
453  buffer << "VPStrategy";
454  return buffer.str();
455  }
456 
458  void PrintInfo(std::ostream &rOStream) const override
459  {
460  rOStream << "VPStrategy";
461  }
462 
464  void PrintData(std::ostream &rOStream) const override
465  {
466  }
467 
471 
473 
474  protected:
477 
481 
485 
489 
493 
495 
499  virtual bool SolveMomentumIteration(unsigned int it, unsigned int maxIt, bool &fixedTimeStep, double &velocityNorm)
500  {
501  return false;
502  }
503 
504  virtual bool SolveContinuityIteration(unsigned int it, unsigned int maxIt, double &NormP)
505  {
506  return false;
507  }
508 
509  void ComputeErrorL2Norm(double tensilStressSign) // tensilStressSign = 1.0 for FIC, tensilStressSign = -1.0 for FS
510  {
511  ModelPart &rModelPart = BaseType::GetModelPart();
512  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
513  const double currentTime = rCurrentProcessInfo[TIME];
514  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
515 
516  long double sumErrorL2Velocity = 0;
517  long double sumErrorL2VelocityX = 0;
518  long double sumErrorL2VelocityY = 0;
519  long double sumErrorL2Pressure = 0;
520  long double sumErrorL2TauXX = 0;
521  long double sumErrorL2TauYY = 0;
522  long double sumErrorL2TauXY = 0;
523 
524 #pragma omp parallel
525  {
526  ModelPart::ElementIterator ElemBegin;
527  ModelPart::ElementIterator ElemEnd;
528  OpenMPUtils::PartitionedIterators(rModelPart.Elements(), ElemBegin, ElemEnd);
529  for (ModelPart::ElementIterator itElem = ElemBegin; itElem != ElemEnd; ++itElem)
530  {
531 
532  Element::GeometryType &geometry = itElem->GetGeometry();
533  long double nodalArea = 0;
534 
535  if (dimension == 2)
536  {
537  nodalArea = geometry.Area() / 3.0;
538  }
539  else if (dimension == 3)
540  {
541  nodalArea = geometry.Volume() * 0.25;
542  }
543 
544  long double bariPosX = 0;
545  long double bariPosY = 0;
546 
547  long double eleErrorL2Velocity = 0;
548  long double eleErrorL2VelocityX = 0;
549  long double eleErrorL2VelocityY = 0;
550  long double eleErrorL2Pressure = 0;
551 
552  // ShapeFunctionDerivativesArrayType DN_DX;
553  Matrix NContainer;
554  NContainer = geometry.ShapeFunctionsValues(GeometryData::IntegrationMethod::GI_GAUSS_1);
555 
556  const Vector &N = row(NContainer, 0);
557  const unsigned int NumNodes = geometry.size();
558 
559  double elementalPressure = N[0] * geometry(0)->FastGetSolutionStepValue(PRESSURE);
560  double elementalVelocityX = N[0] * geometry(0)->FastGetSolutionStepValue(VELOCITY_X);
561  double elementalVelocityY = N[0] * geometry(0)->FastGetSolutionStepValue(VELOCITY_Y);
562  ;
563 
564  for (unsigned int i = 1; i < NumNodes; i++)
565  {
566  elementalPressure += N[i] * geometry(i)->FastGetSolutionStepValue(PRESSURE);
567  elementalVelocityX += N[i] * geometry(i)->FastGetSolutionStepValue(VELOCITY_X);
568  elementalVelocityY += N[i] * geometry(i)->FastGetSolutionStepValue(VELOCITY_Y);
569  }
570 
571  for (unsigned int i = 0; i < geometry.size(); i++)
572  {
573 
574  const long double nodalPosX = geometry(i)->X();
575  const long double nodalPosY = geometry(i)->Y();
576 
577  bariPosX += nodalPosX / 3.0;
578  bariPosY += nodalPosY / 3.0;
579  }
580 
581  const long double posX = bariPosX;
582  const long double posY = bariPosY;
583  long double expectedVelocityX = pow(posX, 2) * (1.0 - posX) * (1.0 - posX) * (2.0 * posY - 6.0 * pow(posY, 2) + 4.0 * pow(posY, 3));
584  long double expectedVelocityY = -pow(posY, 2) * (1.0 - posY) * (1.0 - posY) * (2.0 * posX - 6.0 * pow(posX, 2) + 4.0 * pow(posX, 3));
585  long double expectedPressure = -tensilStressSign * posX * (1.0 - posX);
586 
587  eleErrorL2VelocityX = elementalVelocityX - expectedVelocityX;
588  eleErrorL2VelocityY = elementalVelocityY - expectedVelocityY;
589  eleErrorL2Pressure = elementalPressure - expectedPressure;
590 
591  sumErrorL2VelocityX += pow(eleErrorL2VelocityX, 2) * geometry.Area();
592  sumErrorL2VelocityY += pow(eleErrorL2VelocityY, 2) * geometry.Area();
593  sumErrorL2Pressure += pow(eleErrorL2Pressure, 2) * geometry.Area();
594 
595  const long double tauXX = 0; // itElem->GetValue(ELEMENTAL_DEVIATORIC_STRESS_XX);
596  const long double tauYY = 0; // itElem->GetValue(ELEMENTAL_DEVIATORIC_STRESS_YY);
597  const long double tauXY = 0; // itElem->GetValue(ELEMENTAL_DEVIATORIC_STRESS_XY);
598 
599  long double expectedTauXX = 2.0 * (-4.0 * (1.0 - bariPosX) * bariPosX * (-1.0 + 2.0 * bariPosX) * bariPosY * (1.0 - 3.0 * bariPosY + 2.0 * pow(bariPosY, 2)));
600  long double expectedTauYY = 2.0 * (4.0 * bariPosX * (1.0 - 3.0 * bariPosX + 2.0 * pow(bariPosX, 2)) * (1.0 - bariPosY) * bariPosY * (-1.0 + 2.0 * bariPosY));
601  long double expectedTauXY = (2.0 * (1.0 - 6.0 * bariPosY + 6.0 * pow(bariPosY, 2)) * (1.0 - bariPosX) * (1.0 - bariPosX) * pow(bariPosX, 2) - 2.0 * (1.0 - 6.0 * bariPosX + 6.0 * pow(bariPosX, 2)) * (1.0 - bariPosY) * (1 - bariPosY) * pow(bariPosY, 2));
602 
603  long double nodalErrorTauXX = tauXX - expectedTauXX;
604  long double nodalErrorTauYY = tauYY - expectedTauYY;
605  long double nodalErrorTauXY = tauXY - expectedTauXY;
606 
607  sumErrorL2TauXX += pow(nodalErrorTauXX, 2) * geometry.Area();
608  sumErrorL2TauYY += pow(nodalErrorTauYY, 2) * geometry.Area();
609  sumErrorL2TauXY += pow(nodalErrorTauXY, 2) * geometry.Area();
610  }
611  }
612 
613  long double errorL2Velocity = sqrt(sumErrorL2Velocity);
614  long double errorL2VelocityX = sqrt(sumErrorL2VelocityX);
615  long double errorL2VelocityY = sqrt(sumErrorL2VelocityY);
616  long double errorL2Pressure = sqrt(sumErrorL2Pressure);
617  long double errorL2TauXX = sqrt(sumErrorL2TauXX);
618  long double errorL2TauYY = sqrt(sumErrorL2TauYY);
619  long double errorL2TauXY = sqrt(sumErrorL2TauXY);
620 
621  std::ofstream myfileVelocity;
622  myfileVelocity.open("errorL2VelocityFile.txt", std::ios::app);
623  myfileVelocity << currentTime << "\t" << errorL2Velocity << "\n";
624  myfileVelocity.close();
625 
626  std::ofstream myfileVelocityX;
627  myfileVelocityX.open("errorL2VelocityXFile.txt", std::ios::app);
628  myfileVelocityX << currentTime << "\t" << errorL2VelocityX << "\n";
629  myfileVelocityX.close();
630 
631  std::ofstream myfileVelocityY;
632  myfileVelocityY.open("errorL2VelocityYFile.txt", std::ios::app);
633  myfileVelocityY << currentTime << "\t" << errorL2VelocityY << "\n";
634  myfileVelocityY.close();
635 
636  std::ofstream myfilePressure;
637  myfilePressure.open("errorL2PressureFile.txt", std::ios::app);
638  myfilePressure << currentTime << "\t" << errorL2Pressure << "\n";
639  myfilePressure.close();
640 
641  std::ofstream myfileTauXX;
642  myfileTauXX.open("errorL2TauXXFile.txt", std::ios::app);
643  myfileTauXX << currentTime << "\t" << errorL2TauXX << "\n";
644  myfileTauXX.close();
645 
646  std::ofstream myfileTauYY;
647  myfileTauYY.open("errorL2TauYYFile.txt", std::ios::app);
648  myfileTauYY << currentTime << "\t" << errorL2TauYY << "\n";
649  myfileTauYY.close();
650 
651  std::ofstream myfileTauXY;
652  myfileTauXY.open("errorL2TauXYFile.txt", std::ios::app);
653  myfileTauXY << currentTime << "\t" << errorL2TauXY << "\n";
654  myfileTauXY.close();
655  }
656 
657  void ComputeErrorL2NormCasePoiseuille()
658  {
659  ModelPart &rModelPart = BaseType::GetModelPart();
660  const ProcessInfo &rCurrentProcessInfo = rModelPart.GetProcessInfo();
661  const double currentTime = rCurrentProcessInfo[TIME];
662  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
663 
664  double sumErrorL2VelocityTheta = 0;
665  double sumErrorL2TauTheta = 0;
666 
667  double r_in = 0.2;
668  double R_out = 0.5;
669  double kappa = r_in / R_out;
670  double omega = 0.5;
671  double viscosity = 100.0;
672 
673 #pragma omp parallel
674  {
675  ModelPart::ElementIterator ElemBegin;
676  ModelPart::ElementIterator ElemEnd;
677  OpenMPUtils::PartitionedIterators(rModelPart.Elements(), ElemBegin, ElemEnd);
678  for (ModelPart::ElementIterator itElem = ElemBegin; itElem != ElemEnd; ++itElem)
679  {
680 
681  Element::GeometryType &geometry = itElem->GetGeometry();
682  long double nodalArea = 0;
683 
684  if (dimension == 2)
685  {
686  nodalArea = geometry.Area() / 3.0;
687  }
688  else if (dimension == 3)
689  {
690  nodalArea = geometry.Volume() * 0.25;
691  }
692 
693  long double bariPosX = 0;
694  long double bariPosY = 0;
695 
696  long double eleErrorL2Velocity = 0;
697  long double eleErrorL2VelocityX = 0;
698  long double eleErrorL2VelocityY = 0;
699  long double eleErrorL2Pressure = 0;
700 
701  Matrix NContainer;
702  NContainer = geometry.ShapeFunctionsValues(GeometryData::IntegrationMethod::GI_GAUSS_1);
703 
704  const Vector &N = row(NContainer, 0);
705  // itElem->EvaluateInPoint(elementalPressure,PRESSURE,N);
706  const unsigned int NumNodes = geometry.size();
707 
708  double elementalPressure = N[0] * geometry(0)->FastGetSolutionStepValue(PRESSURE);
709  double elementalVelocityX = N[0] * geometry(0)->FastGetSolutionStepValue(VELOCITY_X);
710  double elementalVelocityY = N[0] * geometry(0)->FastGetSolutionStepValue(VELOCITY_Y);
711  ;
712 
713  for (unsigned int i = 1; i < NumNodes; i++)
714  {
715  elementalPressure += N[i] * geometry(i)->FastGetSolutionStepValue(PRESSURE);
716  elementalVelocityX += N[i] * geometry(i)->FastGetSolutionStepValue(VELOCITY_X);
717  elementalVelocityY += N[i] * geometry(i)->FastGetSolutionStepValue(VELOCITY_Y);
718  }
719 
720  for (unsigned int i = 0; i < geometry.size(); i++)
721  {
722 
723  // index = i*dimension;
724  const long double nodalPosX = geometry(i)->X();
725  const long double nodalPosY = geometry(i)->Y();
726 
727  bariPosX += nodalPosX / 3.0;
728  bariPosY += nodalPosY / 3.0;
729  }
730 
731  const long double posX = bariPosX;
732  const long double posY = bariPosY;
733  const double rPos = sqrt(pow(posX, 2) + pow(posY, 2));
734  const double cosalfa = posX / rPos;
735  const double sinalfa = posY / rPos;
736  const double sin2alfa = 2.0 * cosalfa * sinalfa;
737  const double cos2alfa = 1.0 - 2.0 * pow(sinalfa, 2);
738 
739  double expectedVelocityTheta = pow(kappa, 2) * omega * R_out / (1.0 - pow(kappa, 2)) * (R_out / rPos - rPos / R_out);
740  double computedVelocityTheta = sqrt(pow(elementalVelocityX, 2) + pow(elementalVelocityY, 2));
741  double nodalErrorVelocityTheta = computedVelocityTheta - expectedVelocityTheta;
742 
743  const long double tauXX = 0; // itElem->GetValue(ELEMENTAL_DEVIATORIC_STRESS_XX);
744  const long double tauYY = 0; // itElem->GetValue(ELEMENTAL_DEVIATORIC_STRESS_YY);
745  const long double tauXY = 0; // itElem->GetValue(ELEMENTAL_DEVIATORIC_STRESS_XY);
746 
747  double expectedTauTheta = (2.0 * viscosity * pow(kappa, 2) * omega * pow(R_out, 2)) / (1.0 - pow(kappa, 2)) / pow(rPos, 2);
748  double computedTauTheta = (tauXX - tauYY) * sin2alfa / 2.0 - tauXY * cos2alfa;
749  double nodalErrorTauTheta = computedTauTheta - expectedTauTheta;
750 
751  sumErrorL2VelocityTheta += pow(nodalErrorVelocityTheta, 2) * geometry.Area();
752  sumErrorL2TauTheta += pow(nodalErrorTauTheta, 2) * geometry.Area();
753  }
754  }
755 
756  double errorL2VelocityTheta = sqrt(sumErrorL2VelocityTheta);
757  double errorL2TauTheta = sqrt(sumErrorL2TauTheta);
758 
759  std::ofstream myfileVelocity;
760  myfileVelocity.open("errorL2Poiseuille.txt", std::ios::app);
761  myfileVelocity << currentTime << "\t" << errorL2VelocityTheta << "\t" << errorL2TauTheta << "\n";
762  myfileVelocity.close();
763  }
764 
765  double ComputeVelocityNorm()
766  {
767  ModelPart &rModelPart = BaseType::GetModelPart();
768  const int n_nodes = rModelPart.NumberOfNodes();
769 
770  double NormV = 0.00;
771 
772 #pragma omp parallel for reduction(+ \
773  : NormV)
774  for (int i_node = 0; i_node < n_nodes; ++i_node)
775  {
776  const auto it_node = rModelPart.NodesBegin() + i_node;
777  const auto &r_vel = it_node->FastGetSolutionStepValue(VELOCITY);
778  for (unsigned int d = 0; d < 3; ++d)
779  {
780  NormV += r_vel[d] * r_vel[d];
781  }
782  }
783  NormV = BaseType::GetModelPart().GetCommunicator().GetDataCommunicator().SumAll(NormV);
784  NormV = sqrt(NormV);
785 
786  const double zero_tol = 1.0e-12;
787  if (NormV < zero_tol)
788  NormV = 1.00;
789 
790  return NormV;
791  }
792 
793  double ComputePressureNorm()
794  {
795  ModelPart &rModelPart = BaseType::GetModelPart();
796  const int n_nodes = rModelPart.NumberOfNodes();
797 
798  double NormP = 0.00;
799 
800 #pragma omp parallel for reduction(+ \
801  : NormP)
802  for (int i_node = 0; i_node < n_nodes; ++i_node)
803  {
804  const auto it_node = rModelPart.NodesBegin() + i_node;
805  const double Pr = it_node->FastGetSolutionStepValue(PRESSURE);
806  NormP += Pr * Pr;
807  }
808  NormP = BaseType::GetModelPart().GetCommunicator().GetDataCommunicator().SumAll(NormP);
809  NormP = sqrt(NormP);
810 
811  const double zero_tol = 1.0e-12;
812  if (NormP < zero_tol)
813  NormP = 1.00;
814 
815  return NormP;
816  }
817 
818  virtual bool CheckVelocityConvergence(const double NormDv, double &errorNormDv)
819  {
820  return false;
821  }
822 
823  virtual bool CheckPressureConvergence(const double NormDp, double &errorNormDp, double &NormP)
824  {
825  return false;
826  }
827 
828  virtual bool FixTimeStepMomentum(const double DvErrorNorm, bool &fixedTimeStep)
829  {
830  return false;
831  }
832 
833  virtual bool CheckMomentumConvergence(const double DvErrorNorm, bool &fixedTimeStep)
834  {
835  return false;
836  }
837 
838  virtual bool FixTimeStepContinuity(const double DvErrorNorm, bool &fixedTimeStep)
839  {
840  return false;
841  }
842 
843  virtual bool CheckContinuityConvergence(const double DvErrorNorm, bool &fixedTimeStep)
844  {
845  return false;
846  }
847 
851 
855 
859 
861 
864 
868 
869  // Fractional step index.
870  /* 1 : Momentum step (calculate fractional step velocity)
871  * 2-3 : Unused (reserved for componentwise calculation of frac step velocity)
872  * 4 : Pressure step
873  * 5 : Computation of projections
874  * 6 : End of step velocity
875  */
876  // unsigned int mStepId;
877 
881 
885 
886  virtual void InitializeStrategy(SolverSettingsType &rSolverConfig)
887  {
888  KRATOS_TRY;
889  KRATOS_CATCH("");
890  }
891 
895 
899 
903 
905  VPStrategy &operator=(VPStrategy const &rOther) {}
906 
908  VPStrategy(VPStrategy const &rOther) {}
909 
911 
912  };
913 
917 
919 
921 
922 } // namespace Kratos.
923 
924 #endif // KRATOS_V_P_STRATEGY_H
boundary_normals_calculation_utilities.hpp
Kratos::Communicator::GetDataCommunicator
virtual const DataCommunicator & GetDataCommunicator() const
Definition: communicator.cpp:340
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::size
SizeType size() const
Definition: geometry.h:518
Kratos::Geometry::ShapeFunctionsValues
const Matrix & ShapeFunctionsValues() const
Definition: geometry.h:3393
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::size
size_type size() const
Definition: global_pointers_vector.h:307
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::ElementIterator
MeshType::ElementIterator ElementIterator
Definition: model_part.h:174
Kratos::ModelPart::NumberOfNodes
SizeType NumberOfNodes(IndexType ThisIndex=0) const
Definition: model_part.h:341
Kratos::ModelPart::NodesEnd
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
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::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
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::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::MoveMesh
virtual void MoveMesh()
This function is designed to move the mesh.
Definition: solving_strategy.h:330
Kratos::TwoStepVPSolverSettings
Helper class to define solution strategies for TwoStepVPStrategy.
Definition: solver_settings.h:57
Kratos::VPStrategy
Definition: v_p_strategy.h:59
Kratos::VPStrategy::VPStrategy
VPStrategy(ModelPart &rModelPart, SolverSettingsType &rSolverConfig)
Definition: v_p_strategy.h:73
Kratos::VPStrategy::SetEchoLevel
virtual void SetEchoLevel(int Level) override
This sets the level of echo for the solving strategy.
Definition: v_p_strategy.h:436
Kratos::VPStrategy::Info
std::string Info() const override
Turn back information as a string.
Definition: v_p_strategy.h:450
Kratos::VPStrategy::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(VPStrategy)
Kratos::VPStrategy::FixTimeStepContinuity
virtual bool FixTimeStepContinuity(const double DvErrorNorm, bool &fixedTimeStep)
Definition: v_p_strategy.h:838
Kratos::VPStrategy::VPStrategy
VPStrategy(VPStrategy const &rOther)
Copy constructor.
Definition: v_p_strategy.h:908
Kratos::VPStrategy::PrintData
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: v_p_strategy.h:464
Kratos::VPStrategy::PrintInfo
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: v_p_strategy.h:458
Kratos::VPStrategy::CalculateAccelerations
void CalculateAccelerations()
Definition: v_p_strategy.h:341
Kratos::VPStrategy::CalculateTemporalVariables
virtual void CalculateTemporalVariables()
Definition: v_p_strategy.h:274
Kratos::VPStrategy::CalculatePressureVelocity
void CalculatePressureVelocity()
Definition: v_p_strategy.h:224
Kratos::VPStrategy::CalculatePressureAcceleration
void CalculatePressureAcceleration()
Definition: v_p_strategy.h:249
Kratos::VPStrategy::CheckVelocityConvergence
virtual bool CheckVelocityConvergence(const double NormDv, double &errorNormDv)
Definition: v_p_strategy.h:818
Kratos::VPStrategy::CheckContinuityConvergence
virtual bool CheckContinuityConvergence(const double DvErrorNorm, bool &fixedTimeStep)
Definition: v_p_strategy.h:843
Kratos::VPStrategy::operator=
VPStrategy & operator=(VPStrategy const &rOther)
Assignment operator.
Definition: v_p_strategy.h:905
Kratos::VPStrategy::CheckPressureConvergence
virtual bool CheckPressureConvergence(const double NormDp, double &errorNormDp, double &NormP)
Definition: v_p_strategy.h:823
Kratos::VPStrategy::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: v_p_strategy.h:385
Kratos::VPStrategy::ComputeVelocityNorm
double ComputeVelocityNorm()
Definition: v_p_strategy.h:765
Kratos::VPStrategy::Clear
virtual void Clear() override
Clears the internal storage.
Definition: v_p_strategy.h:430
Kratos::VPStrategy::CheckMomentumConvergence
virtual bool CheckMomentumConvergence(const double DvErrorNorm, bool &fixedTimeStep)
Definition: v_p_strategy.h:833
Kratos::VPStrategy::FixTimeStepMomentum
virtual bool FixTimeStepMomentum(const double DvErrorNorm, bool &fixedTimeStep)
Definition: v_p_strategy.h:828
Kratos::VPStrategy::CalculateDisplacementsAndPorosity
virtual void CalculateDisplacementsAndPorosity()
Definition: v_p_strategy.h:396
Kratos::VPStrategy::FinalizeSolutionStep
virtual void FinalizeSolutionStep() override
Performs all the required operations that should be done (for each step) after solving the solution s...
Definition: v_p_strategy.h:103
Kratos::VPStrategy::InitializeStrategy
virtual void InitializeStrategy(SolverSettingsType &rSolverConfig)
Definition: v_p_strategy.h:886
Kratos::VPStrategy::VPStrategy
VPStrategy(ModelPart &rModelPart, typename TLinearSolver::Pointer pVelocityLinearSolver, typename TLinearSolver::Pointer pPressureLinearSolver, bool ReformDofSet=true, unsigned int DomainSize=2)
Definition: v_p_strategy.h:80
Kratos::VPStrategy::SolverSettingsType
TwoStepVPSolverSettings< TSparseSpace, TDenseSpace, TLinearSolver > SolverSettingsType
Definition: v_p_strategy.h:67
Kratos::VPStrategy::ComputeErrorL2NormCasePoiseuille
void ComputeErrorL2NormCasePoiseuille()
Definition: v_p_strategy.h:657
Kratos::VPStrategy::UpdateTopology
void UpdateTopology(ModelPart &rModelPart, unsigned int echoLevel)
Definition: v_p_strategy.h:107
Kratos::VPStrategy::ComputeErrorL2Norm
void ComputeErrorL2Norm(double tensilStressSign)
Definition: v_p_strategy.h:509
Kratos::VPStrategy::UpdateStressStrain
virtual void UpdateStressStrain()
Definition: v_p_strategy.h:428
Kratos::VPStrategy::SolveContinuityIteration
virtual bool SolveContinuityIteration(unsigned int it, unsigned int maxIt, double &NormP)
Definition: v_p_strategy.h:504
Kratos::VPStrategy::SolveSolutionStep
virtual bool SolveSolutionStep() override
Solves the current step. This function returns true if a solution has been found, false otherwise.
Definition: v_p_strategy.h:98
Kratos::VPStrategy::ComputePressureNorm
double ComputePressureNorm()
Definition: v_p_strategy.h:793
Kratos::VPStrategy::~VPStrategy
virtual ~VPStrategy()
Destructor.
Definition: v_p_strategy.h:91
Kratos::VPStrategy::SolveMomentumIteration
virtual bool SolveMomentumIteration(unsigned int it, unsigned int maxIt, bool &fixedTimeStep, double &velocityNorm)
Calculate the coefficients for time iteration.
Definition: v_p_strategy.h:499
Kratos::VPStrategy::BaseType
SolvingStrategy< TSparseSpace, TDenseSpace > BaseType
Definition: v_p_strategy.h:65
Kratos::VPStrategy::SetBlockedAndIsolatedFlags
void SetBlockedAndIsolatedFlags()
Definition: v_p_strategy.h:116
Kratos::VPStrategy::InitializeSolutionStep
virtual void InitializeSolutionStep() override
Performs all the required operations that should be done (for each step) before solving the solution ...
Definition: v_p_strategy.h:105
Kratos::VPStrategy::Check
virtual int Check() override
Function to perform expensive checks.
Definition: v_p_strategy.h:93
Kratos::array_1d< double, 3 >
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
deprecated_variables.h
mesher_utilities.hpp
model_part.h
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::row
AMatrix::MatrixRow< const TExpressionType > row(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression, std::size_t RowIndex)
Definition: amatrix_interface.h:649
edgebased_var.viscosity
float viscosity
Definition: edgebased_var.py:8
face_heat.Dt
Dt
Definition: face_heat.py:78
generate_total_lagrangian_mixed_volumetric_strain_element.n_nodes
int n_nodes
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:15
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
sensitivityMatrix.N
N
Definition: sensitivityMatrix.py:29
tensors::i
integer i
Definition: TensorModule.f:17
openmp_utils.h
pfem_fluid_dynamics_application_variables.h
solving_strategy.h
Generated on Sat Feb 3 2024 02:10:04 for KratosMultiphysics by doxygen 1.9.1