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_residualbased_elimination_builder_and_solver_continuity.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Riccardo Rossi, Alessandro Franci
11 //
12 //
13 
14 #if !defined(KRATOS_NODAL_RESIDUAL_BASED_ELIMINATION_BUILDER_AND_SOLVER_CONTINUITY)
15 #define KRATOS_NODAL_RESIDUAL_BASED_ELIMINATION_BUILDER_AND_SOLVER_CONTINUITY
16 
17 /* System includes */
18 #include <set>
19 #ifdef _OPENMP
20 #include <omp.h>
21 #endif
22 
23 /* External includes */
24 // #define USE_GOOGLE_HASH
25 #ifdef USE_GOOGLE_HASH
26 #include "sparsehash/dense_hash_set" //included in external libraries
27 #else
28 #include <unordered_set>
29 #endif
30 
31 /* Project includes */
32 #include "utilities/timer.h"
33 #include "includes/define.h"
34 #include "includes/key_hash.h"
35 #include "solving_strategies/builder_and_solvers/builder_and_solver.h"
36 #include "includes/model_part.h"
37 
38 #include "pfem_fluid_dynamics_application_variables.h"
39 
40 namespace Kratos
41 {
42 
45 
49 
53 
57 
61 
74  template <class TSparseSpace,
75  class TDenseSpace, //= DenseSpace<double>,
76  class TLinearSolver //= LinearSolver<TSparseSpace,TDenseSpace>
77  >
78  class NodalResidualBasedEliminationBuilderAndSolverContinuity
79  : public BuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>
80  {
81  public:
84  KRATOS_CLASS_POINTER_DEFINITION(NodalResidualBasedEliminationBuilderAndSolverContinuity);
85 
86  typedef BuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver> BaseType;
87 
88  typedef typename BaseType::TSchemeType TSchemeType;
89 
90  typedef typename BaseType::TDataType TDataType;
91 
92  typedef typename BaseType::DofsArrayType DofsArrayType;
93 
94  typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
95 
96  typedef typename BaseType::TSystemVectorType TSystemVectorType;
97 
98  typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
99 
100  typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
101 
102  typedef typename BaseType::TSystemMatrixPointerType TSystemMatrixPointerType;
103  typedef typename BaseType::TSystemVectorPointerType TSystemVectorPointerType;
104 
105  typedef Node NodeType;
106 
107  typedef typename BaseType::NodesArrayType NodesArrayType;
108  typedef typename BaseType::ElementsArrayType ElementsArrayType;
109  typedef typename BaseType::ConditionsArrayType ConditionsArrayType;
110 
111  typedef typename BaseType::ElementsContainerType ElementsContainerType;
112 
113  typedef Vector VectorType;
114 
118 
121  NodalResidualBasedEliminationBuilderAndSolverContinuity(
122  typename TLinearSolver::Pointer pNewLinearSystemSolver)
123  : BuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>(pNewLinearSystemSolver)
124  {
125  // KRATOS_INFO("NodalResidualBasedEliminationBuilderAndSolverContinuity") << "Using the standard builder and solver " << std::endl;
126  }
127 
130  ~NodalResidualBasedEliminationBuilderAndSolverContinuity() override
131  {
132  }
133 
137 
141 
142  void BuildNodally(
143  typename TSchemeType::Pointer pScheme,
144  ModelPart &rModelPart,
145  TSystemMatrixType &A,
146  TSystemVectorType &b)
147  {
148  KRATOS_TRY
149 
150  KRATOS_ERROR_IF(!pScheme) << "No scheme provided!" << std::endl;
151 
152  // contributions to the continuity equation system
153  LocalSystemMatrixType LHS_Contribution = LocalSystemMatrixType(0, 0);
154  LocalSystemVectorType RHS_Contribution = LocalSystemVectorType(0);
155 
156  Element::EquationIdVectorType EquationId;
157  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
158 
159  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
160  const double timeInterval = CurrentProcessInfo[DELTA_TIME];
161  const double deviatoric_threshold = 0.1;
162  double pressure = 0;
163  double deltaPressure = 0;
164  double meanMeshSize = 0;
165  double characteristicLength = 0;
166  double density = 0;
167  double nodalVelocityNorm = 0;
168  double tauStab = 0;
169  double dNdXi = 0;
170  double dNdYi = 0;
171  double dNdZi = 0;
172  double dNdXj = 0;
173  double dNdYj = 0;
174  double dNdZj = 0;
175  unsigned int firstRow = 0;
176  unsigned int firstCol = 0;
177 
178  /* #pragma omp parallel */
179  {
180  ModelPart::NodeIterator NodesBegin;
181  ModelPart::NodeIterator NodesEnd;
182  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
183 
184  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
185  {
186 
187  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
188  const unsigned int neighSize = neighb_nodes.size() + 1;
189 
190  if (neighSize > 1)
191  {
192 
193  const double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
194  noalias(LHS_Contribution) = ZeroMatrix(neighSize, neighSize);
195  noalias(RHS_Contribution) = ZeroVector(neighSize);
196 
197  if (EquationId.size() != neighSize)
198  EquationId.resize(neighSize, false);
199 
200  double deviatoricCoeff = 0;
201  this->GetDeviatoricCoefficientForFluid(rModelPart, itNode, deviatoricCoeff);
202 
203  if (deviatoricCoeff > deviatoric_threshold)
204  {
205  deviatoricCoeff = deviatoric_threshold;
206  }
207 
208  double volumetricCoeff = timeInterval * itNode->FastGetSolutionStepValue(BULK_MODULUS);
209 
210  deltaPressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
211 
212  LHS_Contribution(0, 0) += nodalVolume / volumetricCoeff;
213  RHS_Contribution[0] += -deltaPressure * nodalVolume / volumetricCoeff;
214 
215  RHS_Contribution[0] += itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) * nodalVolume;
216 
217  const unsigned int xDofPos = itNode->GetDofPosition(PRESSURE);
218  EquationId[0] = itNode->GetDof(PRESSURE, xDofPos).EquationId();
219 
220  for (unsigned int i = 0; i < neighb_nodes.size(); i++)
221  {
222  EquationId[i + 1] = neighb_nodes[i].GetDof(PRESSURE, xDofPos).EquationId();
223  }
224 
225  firstRow = 0;
226  firstCol = 0;
227  meanMeshSize = itNode->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE);
228  characteristicLength = 1.0 * meanMeshSize;
229  density = itNode->FastGetSolutionStepValue(DENSITY);
230 
231  /* double tauStab=1.0/(8.0*deviatoricCoeff/(meanMeshSize*meanMeshSize)+2.0*density/timeInterval); */
232 
233  if (dimension == 2)
234  {
235  nodalVelocityNorm = sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X) * itNode->FastGetSolutionStepValue(VELOCITY_X) +
236  itNode->FastGetSolutionStepValue(VELOCITY_Y) * itNode->FastGetSolutionStepValue(VELOCITY_Y));
237  }
238  else if (dimension == 3)
239  {
240  nodalVelocityNorm = sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X) * itNode->FastGetSolutionStepValue(VELOCITY_X) +
241  itNode->FastGetSolutionStepValue(VELOCITY_Y) * itNode->FastGetSolutionStepValue(VELOCITY_Y) +
242  itNode->FastGetSolutionStepValue(VELOCITY_Z) * itNode->FastGetSolutionStepValue(VELOCITY_Z));
243  }
244 
245  tauStab = 1.0 * (characteristicLength * characteristicLength * timeInterval) / (density * nodalVelocityNorm * timeInterval * characteristicLength + density * characteristicLength * characteristicLength + 8.0 * deviatoricCoeff * timeInterval);
246  itNode->FastGetSolutionStepValue(NODAL_TAU) = tauStab;
247 
248  LHS_Contribution(0, 0) += +nodalVolume * tauStab * density / (volumetricCoeff * timeInterval);
249  RHS_Contribution[0] += -nodalVolume * tauStab * density / (volumetricCoeff * timeInterval) * (deltaPressure - itNode->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) * timeInterval);
250 
251  if (itNode->Is(FREE_SURFACE))
252  {
253  // // double nodalFreesurfaceArea=itNode->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA);
254  // /* LHS_Contribution(0,0) += + 2.0 * tauStab * nodalFreesurfaceArea / meanMeshSize; */
255  // /* RHS_Contribution[0] += - 2.0 * tauStab * nodalFreesurfaceArea / meanMeshSize * itNode->FastGetSolutionStepValue(PRESSURE,0); */
256  LHS_Contribution(0, 0) += +4.0 * tauStab * nodalVolume / (meanMeshSize * meanMeshSize);
257  RHS_Contribution[0] += -4.0 * tauStab * nodalVolume / (meanMeshSize * meanMeshSize) * itNode->FastGetSolutionStepValue(PRESSURE, 0);
258 
259  const array_1d<double, 3> &Normal = itNode->FastGetSolutionStepValue(NORMAL);
260  Vector &SpatialDefRate = itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
261  array_1d<double, 3> nodalAcceleration = 0.5 * (itNode->FastGetSolutionStepValue(VELOCITY, 0) - itNode->FastGetSolutionStepValue(VELOCITY, 1)) / timeInterval - itNode->FastGetSolutionStepValue(ACCELERATION, 1);
262  /* nodalAcceleration= (itNode->FastGetSolutionStepValue(VELOCITY,0)-itNode->FastGetSolutionStepValue(VELOCITY,1))/timeInterval; */
263 
264  double nodalNormalAcceleration = 0;
265  double nodalNormalProjDefRate = 0;
266  if (dimension == 2)
267  {
268  nodalNormalProjDefRate = Normal[0] * SpatialDefRate[0] * Normal[0] + Normal[1] * SpatialDefRate[1] * Normal[1] + 2 * Normal[0] * SpatialDefRate[2] * Normal[1];
269  /* nodalNormalAcceleration=Normal[0]*itNode->FastGetSolutionStepValue(ACCELERATION_X,1) + Normal[1]*itNode->FastGetSolutionStepValue(ACCELERATION_Y,1); */
270  // nodalNormalAcceleration=(0.5*(itNode->FastGetSolutionStepValue(VELOCITY_X,0)-itNode->FastGetSolutionStepValue(VELOCITY_X,1))/timeInterval+0.5*itNode->FastGetSolutionStepValue(ACCELERATION_X,1))*Normal[0] +
271  // (0.5*(itNode->FastGetSolutionStepValue(VELOCITY_Y,0)-itNode->FastGetSolutionStepValue(VELOCITY_Y,1))/timeInterval+0.5*itNode->FastGetSolutionStepValue(ACCELERATION_Y,1))*Normal[1];
272  nodalNormalAcceleration = Normal[0] * nodalAcceleration[0] + Normal[1] * nodalAcceleration[1];
273  }
274  else if (dimension == 3)
275  {
276  nodalNormalProjDefRate = Normal[0] * SpatialDefRate[0] * Normal[0] + Normal[1] * SpatialDefRate[1] * Normal[1] + Normal[2] * SpatialDefRate[2] * Normal[2] +
277  2 * Normal[0] * SpatialDefRate[3] * Normal[1] + 2 * Normal[0] * SpatialDefRate[4] * Normal[2] + 2 * Normal[1] * SpatialDefRate[5] * Normal[2];
278 
279  /* nodalNormalAcceleration=Normal[0]*itNode->FastGetSolutionStepValue(ACCELERATION_X) + Normal[1]*itNode->FastGetSolutionStepValue(ACCELERATION_Y) + Normal[2]*itNode->FastGetSolutionStepValue(ACCELERATION_Z); */
280  /* nodalNormalAcceleration=Normal[0]*nodalAcceleration[0] + Normal[1]*nodalAcceleration[1] + Normal[2]*nodalAcceleration[2]; */
281  }
282  // RHS_Contribution[0] += tauStab * (density*nodalNormalAcceleration - 4.0*deviatoricCoeff*nodalNormalProjDefRate/meanMeshSize) * nodalFreesurfaceArea;
283  double accelerationContribution = 2.0 * density * nodalNormalAcceleration / meanMeshSize;
284  double deviatoricContribution = 8.0 * deviatoricCoeff * nodalNormalProjDefRate / (meanMeshSize * meanMeshSize);
285 
286  RHS_Contribution[0] += 1.0 * tauStab * (accelerationContribution - deviatoricContribution) * nodalVolume;
287  }
288 
289  array_1d<double, 3> &VolumeAcceleration = itNode->FastGetSolutionStepValue(VOLUME_ACCELERATION);
290 
291  // double posX= itNode->X();
292 
293  // double posY= itNode->Y();
294 
295  // double coeffX =(12.0-24.0*posY)*pow(posX,4);
296 
297  // coeffX += (-24.0+48.0*posY)*pow(posX,3);
298 
299  // coeffX += (-48.0*posY+72.0*pow(posY,2)-48.0*pow(posY,3)+12.0)*pow(posX,2);
300 
301  // coeffX += (-2.0+24.0*posY-72.0*pow(posY,2)+48.0*pow(posY,3))*posX;
302 
303  // coeffX += 1.0-4.0*posY+12.0*pow(posY,2)-8.0*pow(posY,3);
304 
305  // double coeffY =(8.0-48.0*posY+48.0*pow(posY,2))*pow(posX,3);
306 
307  // coeffY += (-12.0+72.0*posY-72.0*pow(posY,2))*pow(posX,2);
308 
309  // coeffY += (4.0-24.0*posY+48.0*pow(posY,2)-48.0*pow(posY,3)+24.0*pow(posY,4))*posX;
310 
311  // coeffY += -12.0*pow(posY,2)+24.0*pow(posY,3)-12.0*pow(posY,4);
312 
313  for (unsigned int i = 0; i < neighSize; i++)
314  {
315 
316  dNdXi = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstCol];
317  dNdYi = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstCol + 1];
318 
319  if (i != 0)
320  {
321  // i==0 of EquationIs has been already filled with the master node (that is not included in neighb_nodes). The next is stored for i+1
322  EquationId[i] = neighb_nodes[i - 1].GetDof(PRESSURE, xDofPos).EquationId();
323  // at i==0 density and volume acceleration are taken from the master node
324  density = neighb_nodes[i - 1].FastGetSolutionStepValue(DENSITY);
325  // VolumeAcceleration = neighb_nodes[i-1].FastGetSolutionStepValue(VOLUME_ACCELERATION);
326 
327  // // posX= neighb_nodes[i-1].X();
328 
329  // // posY= neighb_nodes[i-1].Y();
330 
331  // // coeffX =(12.0-24.0*posY)*pow(posX,4);
332 
333  // // coeffX += (-24.0+48.0*posY)*pow(posX,3);
334 
335  // // coeffX += (-48.0*posY+72.0*pow(posY,2)-48.0*pow(posY,3)+12.0)*pow(posX,2);
336 
337  // // coeffX += (-2.0+24.0*posY-72.0*pow(posY,2)+48.0*pow(posY,3))*posX;
338 
339  // // coeffX += 1.0-4.0*posY+12.0*pow(posY,2)-8.0*pow(posY,3);
340 
341  // // coeffY =(8.0-48.0*posY+48.0*pow(posY,2))*pow(posX,3);
342 
343  // // coeffY += (-12.0+72.0*posY-72.0*pow(posY,2))*pow(posX,2);
344 
345  // // coeffY += (4.0-24.0*posY+48.0*pow(posY,2)-48.0*pow(posY,3)+24.0*pow(posY,4))*posX;
346 
347  // // coeffY += -12.0*pow(posY,2)+24.0*pow(posY,3)-12.0*pow(posY,4);
348  }
349 
350  if (dimension == 2)
351  {
352  // RHS_Contribution[i] += - tauStab * density * (dNdXi* VolumeAcceleration[0]*coeffX + dNdYi* VolumeAcceleration[1]*coeffY) * nodalVolume;
353  RHS_Contribution[i] += -tauStab * density * (dNdXi * VolumeAcceleration[0] + dNdYi * VolumeAcceleration[1]) * nodalVolume;
354  }
355  else if (dimension == 3)
356  {
357  dNdZi = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstCol + 2];
358  RHS_Contribution[i] += -tauStab * density * (dNdXi * VolumeAcceleration[0] + dNdYi * VolumeAcceleration[1] + dNdZi * VolumeAcceleration[2]) * nodalVolume;
359  }
360 
361  firstRow = 0;
362 
363  for (unsigned int j = 0; j < neighSize; j++)
364  {
365  dNdXj = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstRow];
366  dNdYj = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstRow + 1];
367  // double Vx=itNode->FastGetSolutionStepValue(VELOCITY_X,0);
368  // double Vy=itNode->FastGetSolutionStepValue(VELOCITY_Y,0);
369  if (j != 0)
370  {
371  pressure = neighb_nodes[j - 1].FastGetSolutionStepValue(PRESSURE, 0);
372  // Vx= neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_X,0);
373  // Vy= neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Y,0);
374  // meanMeshSize=neighb_nodes[j-1].FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE);
375  // characteristicLength=2.0*meanMeshSize;
376  // density=neighb_nodes[j-1].FastGetSolutionStepValue(DENSITY);
377  // if(dimension==2){
378  // nodalVelocityNorm= sqrt(neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_X)*neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_X) +
379  // neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Y)*neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Y));
380  // }else if(dimension==3){
381  // nodalVelocityNorm=sqrt(neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_X)*neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_X) +
382  // neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Y)*neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Y) +
383  // neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Z)*neighb_nodes[j-1].FastGetSolutionStepValue(VELOCITY_Z));
384  // }
385  }
386  else
387  {
388  pressure = itNode->FastGetSolutionStepValue(PRESSURE, 0);
389  // meanMeshSize=itNode->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE);
390  // characteristicLength=2.0*meanMeshSize;
391  // density=itNode->FastGetSolutionStepValue(DENSITY);
392  // if(dimension==2){
393  // nodalVelocityNorm= sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X)*itNode->FastGetSolutionStepValue(VELOCITY_X) +
394  // itNode->FastGetSolutionStepValue(VELOCITY_Y)*itNode->FastGetSolutionStepValue(VELOCITY_Y));
395  // }else if(dimension==3){
396  // nodalVelocityNorm=sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X)*itNode->FastGetSolutionStepValue(VELOCITY_X) +
397  // itNode->FastGetSolutionStepValue(VELOCITY_Y)*itNode->FastGetSolutionStepValue(VELOCITY_Y) +
398  // itNode->FastGetSolutionStepValue(VELOCITY_Z)*itNode->FastGetSolutionStepValue(VELOCITY_Z));
399  // }
400  }
401 
402  // tauStab= 1.0 * (characteristicLength * characteristicLength * timeInterval) / ( density * nodalVelocityNorm * timeInterval * characteristicLength + density * characteristicLength * characteristicLength + 8.0 * deviatoricCoeff * timeInterval );
403 
404  if (dimension == 2)
405  {
406  // // ////////////////// Laplacian term for LHS
407  LHS_Contribution(i, j) += +tauStab * (dNdXi * dNdXj + dNdYi * dNdYj) * nodalVolume;
408  // // ////////////////// Laplacian term L_ij*P_j for RHS
409  RHS_Contribution[i] += -tauStab * (dNdXi * dNdXj + dNdYi * dNdYj) * nodalVolume * pressure;
410 
411  // RHS_Contribution[i] += (dNdXj*Vx + dNdYj*Vy)*nodalVolume/3.0;
412  // LHS_Contribution(i,j)+= nodalVolume/volumetricCoeff/(1.0+double(neighSize));
413  // if(i==j){
414  // RHS_Contribution[i] += (-deltaPressure/volumetricCoeff )*nodalVolume;
415  // }
416  }
417  else if (dimension == 3)
418  {
419  dNdZj = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstRow + 2];
421  LHS_Contribution(i, j) += +tauStab * (dNdXi * dNdXj + dNdYi * dNdYj + dNdZi * dNdZj) * nodalVolume;
423  RHS_Contribution[i] += -tauStab * (dNdXi * dNdXj + dNdYi * dNdYj + dNdZi * dNdZj) * nodalVolume * pressure;
424  }
425  firstRow += dimension;
426  }
427 
428  firstCol += dimension;
429  }
430 
431 #ifdef _OPENMP
432  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId, mlock_array);
433 #else
434  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
435 #endif
436  }
437  }
438  }
439 
440  KRATOS_CATCH("")
441  }
442 
443  void GetDeviatoricCoefficientForFluid(ModelPart &rModelPart, ModelPart::NodeIterator itNode, double &deviatoricCoefficient)
444  {
445  const double tolerance = 1e-12;
446 
447  if (rModelPart.GetNodalSolutionStepVariablesList().Has(STATIC_FRICTION)) // mu(I)-rheology
448  {
449  const double static_friction = itNode->FastGetSolutionStepValue(STATIC_FRICTION);
450  const double dynamic_friction = itNode->FastGetSolutionStepValue(DYNAMIC_FRICTION);
451  const double delta_friction = dynamic_friction - static_friction;
452  const double inertial_number_zero = itNode->FastGetSolutionStepValue(INERTIAL_NUMBER_ZERO);
453  const double grain_diameter = itNode->FastGetSolutionStepValue(GRAIN_DIAMETER);
454  const double grain_density = itNode->FastGetSolutionStepValue(GRAIN_DENSITY);
455  const double regularization_coeff = itNode->FastGetSolutionStepValue(REGULARIZATION_COEFFICIENT);
456 
457  const double theta = 0.5;
458  double mean_pressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) * theta + itNode->FastGetSolutionStepValue(PRESSURE, 1) * (1 - theta);
459 
460  double pressure_tolerance = -1.0e-07;
461  if (mean_pressure > pressure_tolerance)
462  {
463  mean_pressure = pressure_tolerance;
464  }
465 
466  const double equivalent_strain_rate = itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE);
467  const double exponent = -equivalent_strain_rate / regularization_coeff;
468  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);
469 
470  if (std::fabs(equivalent_strain_rate) > tolerance)
471  {
472  const double first_viscous_term = static_friction * (1 - std::exp(exponent)) / equivalent_strain_rate;
473  deviatoricCoefficient = (first_viscous_term + second_viscous_term) * std::fabs(mean_pressure);
474  }
475  else
476  {
477  deviatoricCoefficient = 1.0; // this is for the first iteration and first time step
478  }
479  }
480  else if (rModelPart.GetNodalSolutionStepVariablesList().Has(INTERNAL_FRICTION_ANGLE)) // frictiional viscoplastic model
481  {
482  const double dynamic_viscosity = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
483  const double friction_angle = itNode->FastGetSolutionStepValue(INTERNAL_FRICTION_ANGLE);
484  const double cohesion = itNode->FastGetSolutionStepValue(COHESION);
485  const double adaptive_exponent = itNode->FastGetSolutionStepValue(ADAPTIVE_EXPONENT);
486 
487  const double theta = 0.5;
488  double mean_pressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) * theta + itNode->FastGetSolutionStepValue(PRESSURE, 1) * (1 - theta);
489 
490  double pressure_tolerance = -1.0e-07;
491  if (mean_pressure > pressure_tolerance)
492  {
493  mean_pressure = pressure_tolerance;
494  }
495 
496  const double equivalent_strain_rate = itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE);
497 
498  // Ensuring that the case of equivalent_strain_rate = 0 is not problematic
499  if (std::fabs(equivalent_strain_rate) > tolerance)
500  {
501  const double friction_angle_rad = friction_angle * Globals::Pi / 180.0;
502  const double tanFi = std::tan(friction_angle_rad);
503  double regularization = 1.0 - std::exp(-adaptive_exponent * equivalent_strain_rate);
504  deviatoricCoefficient = dynamic_viscosity + regularization * ((cohesion + tanFi * fabs(mean_pressure)) / equivalent_strain_rate);
505  }
506  else
507  {
508  deviatoricCoefficient = dynamic_viscosity;
509  }
510  }
511  else if (rModelPart.GetNodalSolutionStepVariablesList().Has(YIELD_SHEAR)) // bingham model
512  {
513  const double yieldShear = itNode->FastGetSolutionStepValue(YIELD_SHEAR);
514  const double equivalentStrainRate = itNode->FastGetSolutionStepValue(NODAL_EQUIVALENT_STRAIN_RATE);
515  const double adaptiveExponent = itNode->FastGetSolutionStepValue(ADAPTIVE_EXPONENT);
516  const double exponent = -adaptiveExponent * equivalentStrainRate;
517  deviatoricCoefficient = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
518  if (std::abs(equivalentStrainRate) > tolerance)
519  {
520  deviatoricCoefficient += (yieldShear / equivalentStrainRate) * (1 - exp(exponent));
521  }
522  }
523  else if (rModelPart.GetNodalSolutionStepVariablesList().Has(DYNAMIC_VISCOSITY))
524  {
525  deviatoricCoefficient = itNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
526  }
527  itNode->FastGetSolutionStepValue(DEVIATORIC_COEFFICIENT) = deviatoricCoefficient;
528  }
529 
530  void BuildNodallyUnlessLaplacian(
531  typename TSchemeType::Pointer pScheme,
532  ModelPart &rModelPart,
533  TSystemMatrixType &A,
534  TSystemVectorType &b)
535  {
536  KRATOS_TRY
537 
538  KRATOS_ERROR_IF(!pScheme) << "No scheme provided!" << std::endl;
539 
540  // contributions to the continuity equation system
541  LocalSystemMatrixType LHS_Contribution = LocalSystemMatrixType(0, 0);
542  LocalSystemVectorType RHS_Contribution = LocalSystemVectorType(0);
543 
544  Element::EquationIdVectorType EquationId;
545  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
546 
547  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
548  const double timeInterval = CurrentProcessInfo[DELTA_TIME];
549  const double deviatoric_threshold = 0.1;
550  double deltaPressure = 0;
551  double meanMeshSize = 0;
552  double characteristicLength = 0;
553  double density = 0;
554  double nodalVelocityNorm = 0;
555  double tauStab = 0;
556  double dNdXi = 0;
557  double dNdYi = 0;
558  double dNdZi = 0;
559  unsigned int firstCol = 0;
560 
561  /* #pragma omp parallel */
562  {
563  ModelPart::NodeIterator NodesBegin;
564  ModelPart::NodeIterator NodesEnd;
565  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
566 
567  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
568  {
569 
570  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
571  const unsigned int neighSize = neighb_nodes.size() + 1;
572 
573  if (neighSize > 1)
574  {
575 
576  const double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
577 
578  noalias(LHS_Contribution) = ZeroMatrix(neighSize, neighSize);
579  noalias(RHS_Contribution) = ZeroVector(neighSize);
580 
581  if (EquationId.size() != neighSize)
582  EquationId.resize(neighSize, false);
583 
584  double deviatoricCoeff = 0;
585  this->GetDeviatoricCoefficientForFluid(rModelPart, itNode, deviatoricCoeff);
586 
587  if (deviatoricCoeff > deviatoric_threshold)
588  {
589  deviatoricCoeff = deviatoric_threshold;
590  }
591 
592  double volumetricCoeff = timeInterval * itNode->FastGetSolutionStepValue(BULK_MODULUS);
593 
594  deltaPressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
595 
596  LHS_Contribution(0, 0) += nodalVolume / volumetricCoeff;
597 
598  RHS_Contribution[0] += -deltaPressure * nodalVolume / volumetricCoeff;
599 
600  RHS_Contribution[0] += itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) * nodalVolume;
601 
602  const unsigned int xDofPos = itNode->GetDofPosition(PRESSURE);
603 
604  EquationId[0] = itNode->GetDof(PRESSURE, xDofPos).EquationId();
605 
606  for (unsigned int i = 0; i < neighb_nodes.size(); i++)
607  {
608  EquationId[i + 1] = neighb_nodes[i].GetDof(PRESSURE, xDofPos).EquationId();
609  }
610 
611  firstCol = 0;
612  meanMeshSize = itNode->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE);
613  characteristicLength = 1.0 * meanMeshSize;
614  density = itNode->FastGetSolutionStepValue(DENSITY);
615 
616  /* double tauStab=1.0/(8.0*deviatoricCoeff/(meanMeshSize*meanMeshSize)+2.0*density/timeInterval); */
617 
618  if (dimension == 2)
619  {
620  nodalVelocityNorm = sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X) * itNode->FastGetSolutionStepValue(VELOCITY_X) +
621  itNode->FastGetSolutionStepValue(VELOCITY_Y) * itNode->FastGetSolutionStepValue(VELOCITY_Y));
622  }
623  else if (dimension == 3)
624  {
625  nodalVelocityNorm = sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X) * itNode->FastGetSolutionStepValue(VELOCITY_X) +
626  itNode->FastGetSolutionStepValue(VELOCITY_Y) * itNode->FastGetSolutionStepValue(VELOCITY_Y) +
627  itNode->FastGetSolutionStepValue(VELOCITY_Z) * itNode->FastGetSolutionStepValue(VELOCITY_Z));
628  }
629 
630  tauStab = 1.0 * (characteristicLength * characteristicLength * timeInterval) /
631  (density * nodalVelocityNorm * timeInterval * characteristicLength + density * characteristicLength * characteristicLength + 8.0 * deviatoricCoeff * timeInterval);
632 
633  itNode->FastGetSolutionStepValue(NODAL_TAU) = tauStab;
634 
635  LHS_Contribution(0, 0) += +nodalVolume * tauStab * density / (volumetricCoeff * timeInterval);
636  RHS_Contribution[0] += -nodalVolume * tauStab * density / (volumetricCoeff * timeInterval) * (deltaPressure - itNode->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) * timeInterval);
637 
638  if (itNode->Is(FREE_SURFACE))
639  {
640  // // double nodalFreesurfaceArea=itNode->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA);
641  // /* LHS_Contribution(0,0) += + 2.0 * tauStab * nodalFreesurfaceArea / meanMeshSize; */
642  // /* RHS_Contribution[0] += - 2.0 * tauStab * nodalFreesurfaceArea / meanMeshSize * itNode->FastGetSolutionStepValue(PRESSURE,0); */
643  LHS_Contribution(0, 0) += +4.0 * tauStab * nodalVolume / (meanMeshSize * meanMeshSize);
644  RHS_Contribution[0] += -4.0 * tauStab * nodalVolume / (meanMeshSize * meanMeshSize) * itNode->FastGetSolutionStepValue(PRESSURE, 0);
645 
646  array_1d<double, 3> &Normal = itNode->FastGetSolutionStepValue(NORMAL);
647  Vector &SpatialDefRate = itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
648  array_1d<double, 3> nodalAcceleration = 0.5 * (itNode->FastGetSolutionStepValue(VELOCITY, 0) - itNode->FastGetSolutionStepValue(VELOCITY, 1)) / timeInterval - itNode->FastGetSolutionStepValue(ACCELERATION, 1);
649  /* nodalAcceleration= (itNode->FastGetSolutionStepValue(VELOCITY,0)-itNode->FastGetSolutionStepValue(VELOCITY,1))/timeInterval; */
650 
651  double nodalNormalAcceleration = 0;
652  double nodalNormalProjDefRate = 0;
653  if (dimension == 2)
654  {
655  nodalNormalProjDefRate = Normal[0] * SpatialDefRate[0] * Normal[0] + Normal[1] * SpatialDefRate[1] * Normal[1] + 2 * Normal[0] * SpatialDefRate[2] * Normal[1];
656  /* nodalNormalAcceleration=Normal[0]*itNode->FastGetSolutionStepValue(ACCELERATION_X,1) + Normal[1]*itNode->FastGetSolutionStepValue(ACCELERATION_Y,1); */
657  // nodalNormalAcceleration=(0.5*(itNode->FastGetSolutionStepValue(VELOCITY_X,0)-itNode->FastGetSolutionStepValue(VELOCITY_X,1))/timeInterval+0.5*itNode->FastGetSolutionStepValue(ACCELERATION_X,1))*Normal[0] +
658  // (0.5*(itNode->FastGetSolutionStepValue(VELOCITY_Y,0)-itNode->FastGetSolutionStepValue(VELOCITY_Y,1))/timeInterval+0.5*itNode->FastGetSolutionStepValue(ACCELERATION_Y,1))*Normal[1];
659  nodalNormalAcceleration = Normal[0] * nodalAcceleration[0] + Normal[1] * nodalAcceleration[1];
660  }
661  else if (dimension == 3)
662  {
663  nodalNormalProjDefRate = Normal[0] * SpatialDefRate[0] * Normal[0] + Normal[1] * SpatialDefRate[1] * Normal[1] + Normal[2] * SpatialDefRate[2] * Normal[2] +
664  2 * Normal[0] * SpatialDefRate[3] * Normal[1] + 2 * Normal[0] * SpatialDefRate[4] * Normal[2] + 2 * Normal[1] * SpatialDefRate[5] * Normal[2];
665 
666  /* nodalNormalAcceleration=Normal[0]*itNode->FastGetSolutionStepValue(ACCELERATION_X) + Normal[1]*itNode->FastGetSolutionStepValue(ACCELERATION_Y) + Normal[2]*itNode->FastGetSolutionStepValue(ACCELERATION_Z); */
667  /* nodalNormalAcceleration=Normal[0]*nodalAcceleration[0] + Normal[1]*nodalAcceleration[1] + Normal[2]*nodalAcceleration[2]; */
668  }
669  // RHS_Contribution[0] += tauStab * (density*nodalNormalAcceleration - 4.0*deviatoricCoeff*nodalNormalProjDefRate/meanMeshSize) * nodalFreesurfaceArea;
670  double accelerationContribution = 2.0 * density * nodalNormalAcceleration / meanMeshSize;
671  double deviatoricContribution = 8.0 * deviatoricCoeff * nodalNormalProjDefRate / (meanMeshSize * meanMeshSize);
672 
673  RHS_Contribution[0] += 1.0 * tauStab * (accelerationContribution - deviatoricContribution) * nodalVolume;
674  }
675 
676  array_1d<double, 3> &VolumeAcceleration = itNode->FastGetSolutionStepValue(VOLUME_ACCELERATION);
677 
678  // double posX= itNode->X();
679 
680  // double posY= itNode->Y();
681 
682  // double coeffX =(12.0-24.0*posY)*pow(posX,4);
683 
684  // coeffX += (-24.0+48.0*posY)*pow(posX,3);
685 
686  // coeffX += (-48.0*posY+72.0*pow(posY,2)-48.0*pow(posY,3)+12.0)*pow(posX,2);
687 
688  // coeffX += (-2.0+24.0*posY-72.0*pow(posY,2)+48.0*pow(posY,3))*posX;
689 
690  // coeffX += 1.0-4.0*posY+12.0*pow(posY,2)-8.0*pow(posY,3);
691 
692  // double coeffY =(8.0-48.0*posY+48.0*pow(posY,2))*pow(posX,3);
693 
694  // coeffY += (-12.0+72.0*posY-72.0*pow(posY,2))*pow(posX,2);
695 
696  // coeffY += (4.0-24.0*posY+48.0*pow(posY,2)-48.0*pow(posY,3)+24.0*pow(posY,4))*posX;
697 
698  // coeffY += -12.0*pow(posY,2)+24.0*pow(posY,3)-12.0*pow(posY,4);
699 
700  for (unsigned int i = 0; i < neighSize; i++)
701  {
702 
703  dNdXi = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstCol];
704  dNdYi = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstCol + 1];
705 
706  if (i != 0)
707  {
708  // i==0 of EquationIs has been already filled with the master node (that is not included in neighb_nodes). The next is stored for i+1
709  EquationId[i] = neighb_nodes[i - 1].GetDof(PRESSURE, xDofPos).EquationId();
710  // at i==0 density and volume acceleration are taken from the master node
711  density = neighb_nodes[i - 1].FastGetSolutionStepValue(DENSITY);
712 
713  // // VolumeAcceleration = neighb_nodes[i-1].FastGetSolutionStepValue(VOLUME_ACCELERATION);
714 
715  // // posX= neighb_nodes[i-1].X();
716 
717  // // posY= neighb_nodes[i-1].Y();
718 
719  // // coeffX =(12.0-24.0*posY)*pow(posX,4);
720 
721  // // coeffX += (-24.0+48.0*posY)*pow(posX,3);
722 
723  // // coeffX += (-48.0*posY+72.0*pow(posY,2)-48.0*pow(posY,3)+12.0)*pow(posX,2);
724 
725  // // coeffX += (-2.0+24.0*posY-72.0*pow(posY,2)+48.0*pow(posY,3))*posX;
726 
727  // // coeffX += 1.0-4.0*posY+12.0*pow(posY,2)-8.0*pow(posY,3);
728 
729  // // coeffY =(8.0-48.0*posY+48.0*pow(posY,2))*pow(posX,3);
730 
731  // // coeffY += (-12.0+72.0*posY-72.0*pow(posY,2))*pow(posX,2);
732 
733  // // coeffY += (4.0-24.0*posY+48.0*pow(posY,2)-48.0*pow(posY,3)+24.0*pow(posY,4))*posX;
734 
735  // // coeffY += -12.0*pow(posY,2)+24.0*pow(posY,3)-12.0*pow(posY,4);
736  }
737 
738  if (dimension == 2)
739  {
740  // RHS_Contribution[i] += - tauStab * density * (dNdXi* VolumeAcceleration[0]*coeffX + dNdYi* VolumeAcceleration[1]*coeffY) * nodalVolume;
741  RHS_Contribution[i] += -tauStab * density * (dNdXi * VolumeAcceleration[0] + dNdYi * VolumeAcceleration[1]) * nodalVolume;
742  }
743  else if (dimension == 3)
744  {
745  dNdZi = itNode->FastGetSolutionStepValue(NODAL_SFD_NEIGHBOURS)[firstCol + 2];
746  RHS_Contribution[i] += -tauStab * density * (dNdXi * VolumeAcceleration[0] + dNdYi * VolumeAcceleration[1] + dNdZi * VolumeAcceleration[2]) * nodalVolume;
747  }
748 
749  firstCol += dimension;
750  }
751 
752 #ifdef _OPENMP
753  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId, mlock_array);
754 #else
755  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
756 #endif
757  }
758  }
759  }
760 
761  KRATOS_CATCH("")
762  }
763 
764  void BuildNodallyNoVolumetricStabilizedTerms(
765  typename TSchemeType::Pointer pScheme,
766  ModelPart &rModelPart,
767  TSystemMatrixType &A,
768  TSystemVectorType &b)
769  {
770  KRATOS_TRY
771 
772  KRATOS_ERROR_IF(!pScheme) << "No scheme provided!" << std::endl;
773 
774  // contributions to the continuity equation system
775  LocalSystemMatrixType LHS_Contribution = LocalSystemMatrixType(0, 0);
776  LocalSystemVectorType RHS_Contribution = LocalSystemVectorType(0);
777 
778  Element::EquationIdVectorType EquationId;
779  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
780 
781  const unsigned int dimension = rModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
782  const double timeInterval = CurrentProcessInfo[DELTA_TIME];
783  const double deviatoric_threshold = 0.1;
784  double deltaPressure = 0;
785  double meanMeshSize = 0;
786  double characteristicLength = 0;
787  double density = 0;
788  double nodalVelocityNorm = 0;
789  double tauStab = 0;
790 
791  /* #pragma omp parallel */
792  {
793  ModelPart::NodeIterator NodesBegin;
794  ModelPart::NodeIterator NodesEnd;
795  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
796 
797  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
798  {
799 
800  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
801  const unsigned int neighSize = neighb_nodes.size() + 1;
802 
803  if (neighSize > 1)
804  {
805 
806  const double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
807 
808  noalias(LHS_Contribution) = ZeroMatrix(neighSize, neighSize);
809  noalias(RHS_Contribution) = ZeroVector(neighSize);
810 
811  if (EquationId.size() != neighSize)
812  EquationId.resize(neighSize, false);
813 
814  double deviatoricCoeff = 0;
815  this->GetDeviatoricCoefficientForFluid(rModelPart, itNode, deviatoricCoeff);
816 
817  if (deviatoricCoeff > deviatoric_threshold)
818  {
819  deviatoricCoeff = deviatoric_threshold;
820  }
821 
822  double volumetricCoeff = timeInterval * itNode->FastGetSolutionStepValue(BULK_MODULUS);
823 
824  deltaPressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
825 
826  LHS_Contribution(0, 0) += nodalVolume / volumetricCoeff;
827 
828  RHS_Contribution[0] += -deltaPressure * nodalVolume / volumetricCoeff;
829 
830  RHS_Contribution[0] += itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) * nodalVolume;
831 
832  const unsigned int xDofPos = itNode->GetDofPosition(PRESSURE);
833 
834  EquationId[0] = itNode->GetDof(PRESSURE, xDofPos).EquationId();
835 
836  for (unsigned int i = 0; i < neighb_nodes.size(); i++)
837  {
838  EquationId[i + 1] = neighb_nodes[i].GetDof(PRESSURE, xDofPos).EquationId();
839  }
840 
841  meanMeshSize = itNode->FastGetSolutionStepValue(NODAL_MEAN_MESH_SIZE);
842  characteristicLength = 1.0 * meanMeshSize;
843  density = itNode->FastGetSolutionStepValue(DENSITY);
844 
845  /* double tauStab=1.0/(8.0*deviatoricCoeff/(meanMeshSize*meanMeshSize)+2.0*density/timeInterval); */
846 
847  if (dimension == 2)
848  {
849  nodalVelocityNorm = sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X) * itNode->FastGetSolutionStepValue(VELOCITY_X) +
850  itNode->FastGetSolutionStepValue(VELOCITY_Y) * itNode->FastGetSolutionStepValue(VELOCITY_Y));
851  }
852  else if (dimension == 3)
853  {
854  nodalVelocityNorm = sqrt(itNode->FastGetSolutionStepValue(VELOCITY_X) * itNode->FastGetSolutionStepValue(VELOCITY_X) +
855  itNode->FastGetSolutionStepValue(VELOCITY_Y) * itNode->FastGetSolutionStepValue(VELOCITY_Y) +
856  itNode->FastGetSolutionStepValue(VELOCITY_Z) * itNode->FastGetSolutionStepValue(VELOCITY_Z));
857  }
858 
859  tauStab = 1.0 * (characteristicLength * characteristicLength * timeInterval) /
860  (density * nodalVelocityNorm * timeInterval * characteristicLength + density * characteristicLength * characteristicLength + 8.0 * deviatoricCoeff * timeInterval);
861 
862  itNode->FastGetSolutionStepValue(NODAL_TAU) = tauStab;
863 
864  LHS_Contribution(0, 0) += +nodalVolume * tauStab * density / (volumetricCoeff * timeInterval);
865  RHS_Contribution[0] += -nodalVolume * tauStab * density / (volumetricCoeff * timeInterval) * (deltaPressure - itNode->FastGetSolutionStepValue(PRESSURE_VELOCITY, 0) * timeInterval);
866 
867  if (itNode->Is(FREE_SURFACE))
868  {
869  // // double nodalFreesurfaceArea=itNode->FastGetSolutionStepValue(NODAL_FREESURFACE_AREA);
870  // /* LHS_Contribution(0,0) += + 2.0 * tauStab * nodalFreesurfaceArea / meanMeshSize; */
871  // /* RHS_Contribution[0] += - 2.0 * tauStab * nodalFreesurfaceArea / meanMeshSize * itNode->FastGetSolutionStepValue(PRESSURE,0); */
872 
873  LHS_Contribution(0, 0) += +4.0 * tauStab * nodalVolume / (meanMeshSize * meanMeshSize);
874  RHS_Contribution[0] += -4.0 * tauStab * nodalVolume / (meanMeshSize * meanMeshSize) * itNode->FastGetSolutionStepValue(PRESSURE, 0);
875 
876  array_1d<double, 3> &Normal = itNode->FastGetSolutionStepValue(NORMAL);
877  Vector &SpatialDefRate = itNode->FastGetSolutionStepValue(NODAL_SPATIAL_DEF_RATE);
878  array_1d<double, 3> nodalAcceleration = 0.5 * (itNode->FastGetSolutionStepValue(VELOCITY, 0) - itNode->FastGetSolutionStepValue(VELOCITY, 1)) / timeInterval - itNode->FastGetSolutionStepValue(ACCELERATION, 1);
879  /* nodalAcceleration= (itNode->FastGetSolutionStepValue(VELOCITY,0)-itNode->FastGetSolutionStepValue(VELOCITY,1))/timeInterval; */
880 
881  double nodalNormalAcceleration = 0;
882  double nodalNormalProjDefRate = 0;
883  if (dimension == 2)
884  {
885  nodalNormalProjDefRate = Normal[0] * SpatialDefRate[0] * Normal[0] + Normal[1] * SpatialDefRate[1] * Normal[1] + 2 * Normal[0] * SpatialDefRate[2] * Normal[1];
886  /* nodalNormalAcceleration=Normal[0]*itNode->FastGetSolutionStepValue(ACCELERATION_X,1) + Normal[1]*itNode->FastGetSolutionStepValue(ACCELERATION_Y,1); */
887  // nodalNormalAcceleration=(0.5*(itNode->FastGetSolutionStepValue(VELOCITY_X,0)-itNode->FastGetSolutionStepValue(VELOCITY_X,1))/timeInterval+0.5*itNode->FastGetSolutionStepValue(ACCELERATION_X,1))*Normal[0] +
888  // (0.5*(itNode->FastGetSolutionStepValue(VELOCITY_Y,0)-itNode->FastGetSolutionStepValue(VELOCITY_Y,1))/timeInterval+0.5*itNode->FastGetSolutionStepValue(ACCELERATION_Y,1))*Normal[1];
889  nodalNormalAcceleration = Normal[0] * nodalAcceleration[0] + Normal[1] * nodalAcceleration[1];
890  }
891  else if (dimension == 3)
892  {
893  nodalNormalProjDefRate = Normal[0] * SpatialDefRate[0] * Normal[0] + Normal[1] * SpatialDefRate[1] * Normal[1] + Normal[2] * SpatialDefRate[2] * Normal[2] +
894  2 * Normal[0] * SpatialDefRate[3] * Normal[1] + 2 * Normal[0] * SpatialDefRate[4] * Normal[2] + 2 * Normal[1] * SpatialDefRate[5] * Normal[2];
895 
896  /* nodalNormalAcceleration=Normal[0]*itNode->FastGetSolutionStepValue(ACCELERATION_X) + Normal[1]*itNode->FastGetSolutionStepValue(ACCELERATION_Y) + Normal[2]*itNode->FastGetSolutionStepValue(ACCELERATION_Z); */
897  /* nodalNormalAcceleration=Normal[0]*nodalAcceleration[0] + Normal[1]*nodalAcceleration[1] + Normal[2]*nodalAcceleration[2]; */
898  }
899  // RHS_Contribution[0] += tauStab * (density*nodalNormalAcceleration - 4.0*deviatoricCoeff*nodalNormalProjDefRate/meanMeshSize) * nodalFreesurfaceArea;
900  double accelerationContribution = 2.0 * density * nodalNormalAcceleration / meanMeshSize;
901  double deviatoricContribution = 8.0 * deviatoricCoeff * nodalNormalProjDefRate / (meanMeshSize * meanMeshSize);
902 
903  RHS_Contribution[0] += 1.0 * tauStab * (accelerationContribution - deviatoricContribution) * nodalVolume;
904  }
905 
906 #ifdef _OPENMP
907  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId, mlock_array);
908 #else
909  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
910 #endif
911  }
912  }
913  }
914 
915  KRATOS_CATCH("")
916  }
917 
918  void BuildNodallyNotStabilized(
919  typename TSchemeType::Pointer pScheme,
920  ModelPart &rModelPart,
921  TSystemMatrixType &A,
922  TSystemVectorType &b)
923  {
924  KRATOS_TRY
925 
926  KRATOS_ERROR_IF(!pScheme) << "No scheme provided!" << std::endl;
927 
928  // contributions to the continuity equation system
929  LocalSystemMatrixType LHS_Contribution = LocalSystemMatrixType(0, 0);
930  LocalSystemVectorType RHS_Contribution = LocalSystemVectorType(0);
931 
932  Element::EquationIdVectorType EquationId;
933  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
934 
935  const double timeInterval = CurrentProcessInfo[DELTA_TIME];
936  const double deviatoric_threshold = 0.1;
937  double deltaPressure = 0;
938 
939  /* #pragma omp parallel */
940  {
941  ModelPart::NodeIterator NodesBegin;
942  ModelPart::NodeIterator NodesEnd;
943  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
944 
945  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
946  {
947 
948  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
949  const unsigned int neighSize = neighb_nodes.size() + 1;
950 
951  if (neighSize > 1)
952  {
953 
954  const double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
955 
956  noalias(LHS_Contribution) = ZeroMatrix(neighSize, neighSize);
957  noalias(RHS_Contribution) = ZeroVector(neighSize);
958 
959  if (EquationId.size() != neighSize)
960  EquationId.resize(neighSize, false);
961 
962  double deviatoricCoeff = 0;
963  this->GetDeviatoricCoefficientForFluid(rModelPart, itNode, deviatoricCoeff);
964 
965  if (deviatoricCoeff > deviatoric_threshold)
966  {
967  deviatoricCoeff = deviatoric_threshold;
968  }
969 
970  double volumetricCoeff = timeInterval * itNode->FastGetSolutionStepValue(BULK_MODULUS);
971 
972  deltaPressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
973 
974  LHS_Contribution(0, 0) += nodalVolume / volumetricCoeff;
975 
976  RHS_Contribution[0] += -deltaPressure * nodalVolume / volumetricCoeff;
977 
978  RHS_Contribution[0] += itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) * nodalVolume;
979 
980  const unsigned int xDofPos = itNode->GetDofPosition(PRESSURE);
981 
982  EquationId[0] = itNode->GetDof(PRESSURE, xDofPos).EquationId();
983 
984  for (unsigned int i = 0; i < neighb_nodes.size(); i++)
985  {
986  EquationId[i + 1] = neighb_nodes[i].GetDof(PRESSURE, xDofPos).EquationId();
987  }
988 
989 #ifdef _OPENMP
990  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId, mlock_array);
991 #else
992  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
993 #endif
994  }
995  }
996  }
997 
998  KRATOS_CATCH("")
999  }
1000 
1001  void BuildAll(
1002  typename TSchemeType::Pointer pScheme,
1003  ModelPart &rModelPart,
1004  TSystemMatrixType &A,
1005  TSystemVectorType &b)
1006  {
1007  KRATOS_TRY
1008 
1009  KRATOS_ERROR_IF(!pScheme) << "No scheme provided!" << std::endl;
1010 
1011  // contributions to the continuity equation system
1012  LocalSystemMatrixType LHS_Contribution = LocalSystemMatrixType(0, 0);
1013  LocalSystemVectorType RHS_Contribution = LocalSystemVectorType(0);
1014 
1015  Element::EquationIdVectorType EquationId;
1016 
1017  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
1018 
1019  const double timeInterval = CurrentProcessInfo[DELTA_TIME];
1020  const double deviatoric_threshold = 0.1;
1021  double deltaPressure = 0;
1022 
1023  /* #pragma omp parallel */
1024  // {
1025  ModelPart::NodeIterator NodesBegin;
1026  ModelPart::NodeIterator NodesEnd;
1027  OpenMPUtils::PartitionedIterators(rModelPart.Nodes(), NodesBegin, NodesEnd);
1028 
1029  for (ModelPart::NodeIterator itNode = NodesBegin; itNode != NodesEnd; ++itNode)
1030  {
1031 
1032  NodeWeakPtrVectorType &neighb_nodes = itNode->GetValue(NEIGHBOUR_NODES);
1033  const unsigned int neighSize = neighb_nodes.size() + 1;
1034 
1035  if (neighSize > 1)
1036  {
1037 
1038  if (LHS_Contribution.size1() != 1)
1039  LHS_Contribution.resize(1, 1, false); // false says not to preserve existing storage!!
1040 
1041  if (RHS_Contribution.size() != 1)
1042  RHS_Contribution.resize(1, false); // false says not to preserve existing storage!!
1043 
1044  noalias(LHS_Contribution) = ZeroMatrix(1, 1);
1045  noalias(RHS_Contribution) = ZeroVector(1);
1046 
1047  if (EquationId.size() != 1)
1048  EquationId.resize(1, false);
1049 
1050  double nodalVolume = itNode->FastGetSolutionStepValue(NODAL_VOLUME);
1051 
1052  if (nodalVolume > 0)
1053  { // in interface nodes not in contact with fluid elements the nodal volume is zero
1054 
1055  double deviatoricCoeff = 0;
1056  this->GetDeviatoricCoefficientForFluid(rModelPart, itNode, deviatoricCoeff);
1057 
1058  if (deviatoricCoeff > deviatoric_threshold)
1059  {
1060  deviatoricCoeff = deviatoric_threshold;
1061  }
1062 
1063  double volumetricCoeff = timeInterval * itNode->FastGetSolutionStepValue(BULK_MODULUS);
1064 
1065  deltaPressure = itNode->FastGetSolutionStepValue(PRESSURE, 0) - itNode->FastGetSolutionStepValue(PRESSURE, 1);
1066 
1067  LHS_Contribution(0, 0) += nodalVolume / volumetricCoeff;
1068 
1069  RHS_Contribution[0] += -deltaPressure * nodalVolume / volumetricCoeff;
1070 
1071  RHS_Contribution[0] += itNode->GetSolutionStepValue(NODAL_VOLUMETRIC_DEF_RATE) * nodalVolume;
1072  }
1073 
1074  const unsigned int xDofPos = itNode->GetDofPosition(PRESSURE);
1075 
1076  EquationId[0] = itNode->GetDof(PRESSURE, xDofPos).EquationId();
1077 
1078 #ifdef _OPENMP
1079  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId, mlock_array);
1080 #else
1081  Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
1082 #endif
1083  }
1084  //}
1085  }
1086 
1087  // }
1088 
1089  ElementsArrayType &pElements = rModelPart.Elements();
1090  int number_of_threads = ParallelUtilities::GetNumThreads();
1091 
1092 #ifdef _OPENMP
1093  int A_size = A.size1();
1094 
1095  // creating an array of lock variables of the size of the system matrix
1096  std::vector<omp_lock_t> lock_array(A.size1());
1097 
1098  for (int i = 0; i < A_size; i++)
1099  omp_init_lock(&lock_array[i]);
1100 #endif
1101 
1102  DenseVector<unsigned int> element_partition;
1103  CreatePartition(number_of_threads, pElements.size(), element_partition);
1104  if (this->GetEchoLevel() > 0)
1105  {
1106  KRATOS_WATCH(number_of_threads);
1107  KRATOS_WATCH(element_partition);
1108  }
1109 
1110 #pragma omp parallel for firstprivate(number_of_threads) schedule(static, 1)
1111  for (int k = 0; k < number_of_threads; k++)
1112  {
1113  // contributions to the system
1114  LocalSystemMatrixType elementalLHS_Contribution = LocalSystemMatrixType(0, 0);
1115  LocalSystemVectorType elementalRHS_Contribution = LocalSystemVectorType(0);
1116 
1117  // vector containing the localization in the system of the different
1118  // terms
1119  Element::EquationIdVectorType elementalEquationId;
1120  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
1121  typename ElementsArrayType::ptr_iterator it_begin = pElements.ptr_begin() + element_partition[k];
1122  typename ElementsArrayType::ptr_iterator it_end = pElements.ptr_begin() + element_partition[k + 1];
1123 
1124  unsigned int pos = (rModelPart.Nodes().begin())->GetDofPosition(PRESSURE);
1125 
1126  // assemble all elements
1127  for (typename ElementsArrayType::ptr_iterator it = it_begin; it != it_end; ++it)
1128  {
1129  // calculate elemental contribution
1130  (*it)->CalculateLocalSystem(elementalLHS_Contribution, elementalRHS_Contribution, CurrentProcessInfo);
1131 
1132  Geometry<Node> &geom = (*it)->GetGeometry();
1133  if (elementalEquationId.size() != geom.size())
1134  elementalEquationId.resize(geom.size(), false);
1135 
1136  for (unsigned int i = 0; i < geom.size(); i++)
1137  elementalEquationId[i] = geom[i].GetDof(PRESSURE, pos).EquationId();
1138 
1139  // assemble the elemental contribution
1140 #ifdef _OPENMP
1141  this->Assemble(A, b, elementalLHS_Contribution, elementalRHS_Contribution, elementalEquationId, lock_array);
1142 #else
1143  this->Assemble(A, b, elementalLHS_Contribution, elementalRHS_Contribution, elementalEquationId);
1144 #endif
1145  }
1146  }
1147 
1148 #ifdef _OPENMP
1149  for (int i = 0; i < A_size; i++)
1150  omp_destroy_lock(&lock_array[i]);
1151 #endif
1152 
1153  KRATOS_CATCH("")
1154  }
1161  void SystemSolve(
1162  TSystemMatrixType &A,
1163  TSystemVectorType &Dx,
1164  TSystemVectorType &b) override
1165  {
1166  KRATOS_TRY
1167 
1168  double norm_b;
1169  if (TSparseSpace::Size(b) != 0)
1170  norm_b = TSparseSpace::TwoNorm(b);
1171  else
1172  norm_b = 0.00;
1173 
1174  if (norm_b != 0.00)
1175  {
1176  // do solve
1177  BaseType::mpLinearSystemSolver->Solve(A, Dx, b);
1178  }
1179  else
1180  TSparseSpace::SetToZero(Dx);
1181 
1182  // Prints informations about the current time
1183  KRATOS_INFO_IF("NodalResidualBasedEliminationBuilderAndSolverContinuity", this->GetEchoLevel() > 1) << *(BaseType::mpLinearSystemSolver) << std::endl;
1184 
1185  KRATOS_CATCH("")
1186  }
1187 
1195  void SystemSolveWithPhysics(
1196  TSystemMatrixType &A,
1197  TSystemVectorType &Dx,
1198  TSystemVectorType &b,
1199  ModelPart &rModelPart)
1200  {
1201  KRATOS_TRY
1202 
1203  double norm_b;
1204  if (TSparseSpace::Size(b) != 0)
1205  norm_b = TSparseSpace::TwoNorm(b);
1206  else
1207  norm_b = 0.00;
1208 
1209  if (norm_b != 0.00)
1210  {
1211  // provide physical data as needed
1212  if (BaseType::mpLinearSystemSolver->AdditionalPhysicalDataIsNeeded())
1213  BaseType::mpLinearSystemSolver->ProvideAdditionalData(A, Dx, b, BaseType::mDofSet, rModelPart);
1214 
1215  // do solve
1216  BaseType::mpLinearSystemSolver->Solve(A, Dx, b);
1217  }
1218  else
1219  {
1220  TSparseSpace::SetToZero(Dx);
1221  KRATOS_WARNING_IF("NodalResidualBasedEliminationBuilderAndSolverContinuity", rModelPart.GetCommunicator().MyPID() == 0) << "ATTENTION! setting the RHS to zero!" << std::endl;
1222  }
1223 
1224  // Prints informations about the current time
1225  KRATOS_INFO_IF("NodalResidualBasedEliminationBuilderAndSolverContinuity", this->GetEchoLevel() > 1 && rModelPart.GetCommunicator().MyPID() == 0) << *(BaseType::mpLinearSystemSolver) << std::endl;
1226 
1227  KRATOS_CATCH("")
1228  }
1229 
1240  void BuildAndSolve(
1241  typename TSchemeType::Pointer pScheme,
1242  ModelPart &rModelPart,
1243  TSystemMatrixType &A,
1244  TSystemVectorType &Dx,
1245  TSystemVectorType &b) override
1246  {
1247  KRATOS_TRY
1248 
1249  Timer::Start("Build");
1250 
1251  /* boost::timer c_build_time; */
1252 
1254  // BuildNodally(pScheme, rModelPart, A, b);
1256 
1257  // /////////////////////// NODAL + ELEMENTAL LAPLACIAN ///////////////////////
1258  // BuildNodallyUnlessLaplacian(pScheme, rModelPart, A, b);
1259  // Build(pScheme, rModelPart, A, b);
1260  // /////////////////////// NODAL + ELEMENTAL LAPLACIAN ///////////////////////
1261 
1263  // BuildNodallyNoVolumetricStabilizedTerms(pScheme, rModelPart, A, b);
1264  // Build(pScheme, rModelPart, A, b);
1265  // /////////////////////// NODAL + ELEMENTAL LAPLACIAN ///////////////////////
1266 
1268  // BuildNodallyNotStabilized(pScheme, rModelPart, A, b);
1269  // Build(pScheme, rModelPart, A, b);
1270 
1271  BuildAll(pScheme, rModelPart, A, b);
1273 
1275  // Build(pScheme, rModelPart, A, b);
1277 
1278  Timer::Stop("Build");
1279 
1280  // ApplyPointLoads(pScheme,rModelPart,b);
1281 
1282  // Does nothing...dirichlet conditions are naturally dealt with in defining the residual
1283  ApplyDirichletConditions(pScheme, rModelPart, A, Dx, b);
1284 
1285  KRATOS_INFO_IF("ResidualBasedBlockBuilderAndSolver", (this->GetEchoLevel() == 3)) << "Before the solution of the system"
1286  << "\nSystem Matrix = " << A << "\nUnknowns vector = " << Dx << "\nRHS vector = " << b << std::endl;
1287 
1288  /* const double start_solve = OpenMPUtils::GetCurrentTime(); */
1289  Timer::Start("Solve");
1290 
1291  /* boost::timer c_solve_time; */
1292  SystemSolveWithPhysics(A, Dx, b, rModelPart);
1293  /* std::cout << "CONTINUITY EQ: solve_time : " << c_solve_time.elapsed() << std::endl; */
1294 
1295  Timer::Stop("Solve");
1296  /* const double stop_solve = OpenMPUtils::GetCurrentTime(); */
1297 
1298  KRATOS_INFO_IF("ResidualBasedBlockBuilderAndSolver", (this->GetEchoLevel() == 3)) << "After the solution of the system"
1299  << "\nSystem Matrix = " << A << "\nUnknowns vector = " << Dx << "\nRHS vector = " << b << std::endl;
1300 
1301  KRATOS_CATCH("")
1302  }
1303 
1304  void Build(
1305  typename TSchemeType::Pointer pScheme,
1306  ModelPart &r_model_part,
1307  TSystemMatrixType &A,
1308  TSystemVectorType &b) override
1309  {
1310  KRATOS_TRY
1311  if (!pScheme)
1312  KRATOS_THROW_ERROR(std::runtime_error, "No scheme provided!", "");
1313 
1314  // getting the elements from the model
1315  ElementsArrayType &pElements = r_model_part.Elements();
1316 
1317  // //getting the array of the conditions
1318  // ConditionsArrayType& ConditionsArray = r_model_part.Conditions();
1319 
1320  // resetting to zero the vector of reactions
1321  TSparseSpace::SetToZero(*(BaseType::mpReactionsVector));
1322 
1323  // create a partition of the element array
1324  int number_of_threads = ParallelUtilities::GetNumThreads();
1325 
1326 #ifdef _OPENMP
1327  int A_size = A.size1();
1328 
1329  // creating an array of lock variables of the size of the system matrix
1330  std::vector<omp_lock_t> lock_array(A.size1());
1331 
1332  for (int i = 0; i < A_size; i++)
1333  omp_init_lock(&lock_array[i]);
1334 #endif
1335 
1336  DenseVector<unsigned int> element_partition;
1337  CreatePartition(number_of_threads, pElements.size(), element_partition);
1338  if (this->GetEchoLevel() > 0)
1339  {
1340  KRATOS_WATCH(number_of_threads);
1341  KRATOS_WATCH(element_partition);
1342  }
1343 
1344  // double start_prod = OpenMPUtils::GetCurrentTime();
1345 
1346 #pragma omp parallel for firstprivate(number_of_threads) schedule(static, 1)
1347  for (int k = 0; k < number_of_threads; k++)
1348  {
1349  // contributions to the system
1350  LocalSystemMatrixType LHS_Contribution = LocalSystemMatrixType(0, 0);
1351  LocalSystemVectorType RHS_Contribution = LocalSystemVectorType(0);
1352 
1353  // vector containing the localization in the system of the different
1354  // terms
1355  Element::EquationIdVectorType EquationId;
1356  const ProcessInfo &CurrentProcessInfo = r_model_part.GetProcessInfo();
1357  typename ElementsArrayType::ptr_iterator it_begin = pElements.ptr_begin() + element_partition[k];
1358  typename ElementsArrayType::ptr_iterator it_end = pElements.ptr_begin() + element_partition[k + 1];
1359 
1360  unsigned int pos = (r_model_part.Nodes().begin())->GetDofPosition(PRESSURE);
1361 
1362  // assemble all elements
1363  for (typename ElementsArrayType::ptr_iterator it = it_begin; it != it_end; ++it)
1364  {
1365 
1366  // calculate elemental contribution
1367  (*it)->CalculateLocalSystem(LHS_Contribution, RHS_Contribution, CurrentProcessInfo);
1368 
1369  Geometry<Node> &geom = (*it)->GetGeometry();
1370  if (EquationId.size() != geom.size())
1371  EquationId.resize(geom.size(), false);
1372 
1373  for (unsigned int i = 0; i < geom.size(); i++)
1374  EquationId[i] = geom[i].GetDof(PRESSURE, pos).EquationId();
1375 
1376  // assemble the elemental contribution
1377 #ifdef _OPENMP
1378  this->Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId, lock_array);
1379 #else
1380  this->Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
1381 #endif
1382  }
1383  }
1384 
1385  // if (this->GetEchoLevel() > 0)
1386  // {
1387  // double stop_prod = OpenMPUtils::GetCurrentTime();
1388  // std::cout << "parallel building time: " << stop_prod - start_prod << std::endl;
1389  // }
1390 
1391 #ifdef _OPENMP
1392  for (int i = 0; i < A_size; i++)
1393  omp_destroy_lock(&lock_array[i]);
1394 #endif
1395 
1396  KRATOS_CATCH("")
1397  }
1398 
1407  void SetUpDofSet(
1408  typename TSchemeType::Pointer pScheme,
1409  ModelPart &rModelPart) override
1410  {
1411  KRATOS_TRY;
1412 
1413  KRATOS_INFO_IF("NodalResidualBasedEliminationBuilderAndSolverContinuity", this->GetEchoLevel() > 1 && rModelPart.GetCommunicator().MyPID() == 0) << "Setting up the dofs" << std::endl;
1414 
1415  // Gets the array of elements from the modeler
1416  ElementsArrayType &pElements = rModelPart.Elements();
1417  const int nelements = static_cast<int>(pElements.size());
1418 
1419  Element::DofsVectorType ElementalDofList;
1420 
1421  const ProcessInfo &CurrentProcessInfo = rModelPart.GetProcessInfo();
1422 
1423  unsigned int nthreads = ParallelUtilities::GetNumThreads();
1424 
1425  // typedef boost::fast_pool_allocator< NodeType::DofType::Pointer > allocator_type;
1426  // typedef std::unordered_set < NodeType::DofType::Pointer,
1427  // DofPointerHasher,
1428  // DofPointerComparor,
1429  // allocator_type > set_type;
1430 
1431 #ifdef USE_GOOGLE_HASH
1432  typedef google::dense_hash_set<NodeType::DofType::Pointer, DofPointerHasher> set_type;
1433 #else
1434  typedef std::unordered_set<NodeType::DofType::Pointer, DofPointerHasher> set_type;
1435 #endif
1436  //
1437 
1438  std::vector<set_type> dofs_aux_list(nthreads);
1439  // std::vector<allocator_type> allocators(nthreads);
1440 
1441  for (int i = 0; i < static_cast<int>(nthreads); i++)
1442  {
1443 #ifdef USE_GOOGLE_HASH
1444  dofs_aux_list[i].set_empty_key(NodeType::DofType::Pointer());
1445 #else
1446  // dofs_aux_list[i] = set_type( allocators[i]);
1447  dofs_aux_list[i].reserve(nelements);
1448 #endif
1449  }
1450 
1451  // #pragma omp parallel for firstprivate(nelements, ElementalDofList)
1452  for (int i = 0; i < static_cast<int>(nelements); ++i)
1453  {
1454  auto it_elem = pElements.begin() + i;
1455  const IndexType this_thread_id = OpenMPUtils::ThisThread();
1456 
1457  // Gets list of Dof involved on every element
1458  pScheme->GetDofList(*it_elem, ElementalDofList, CurrentProcessInfo);
1459 
1460  dofs_aux_list[this_thread_id].insert(ElementalDofList.begin(), ElementalDofList.end());
1461  }
1462 
1463  // here we do a reduction in a tree so to have everything on thread 0
1464  unsigned int old_max = nthreads;
1465  unsigned int new_max = ceil(0.5 * static_cast<double>(old_max));
1466  while (new_max >= 1 && new_max != old_max)
1467  {
1468 
1469 #pragma omp parallel for
1470  for (int i = 0; i < static_cast<int>(new_max); i++)
1471  {
1472  if (i + new_max < old_max)
1473  {
1474  dofs_aux_list[i].insert(dofs_aux_list[i + new_max].begin(), dofs_aux_list[i + new_max].end());
1475  dofs_aux_list[i + new_max].clear();
1476  }
1477  }
1478 
1479  old_max = new_max;
1480  new_max = ceil(0.5 * static_cast<double>(old_max));
1481  }
1482 
1483  DofsArrayType Doftemp;
1484  BaseType::mDofSet = DofsArrayType();
1485 
1486  Doftemp.reserve(dofs_aux_list[0].size());
1487  for (auto it = dofs_aux_list[0].begin(); it != dofs_aux_list[0].end(); it++)
1488  {
1489  Doftemp.push_back((*it));
1490  }
1491  Doftemp.Sort();
1492 
1493  BaseType::mDofSet = Doftemp;
1494 
1495  // Throws an execption if there are no Degrees of freedom involved in the analysis
1496  KRATOS_ERROR_IF(BaseType::mDofSet.size() == 0) << "No degrees of freedom!" << std::endl;
1497 
1498  BaseType::mDofSetIsInitialized = true;
1499 
1500  KRATOS_INFO_IF("NodalResidualBasedEliminationBuilderAndSolverContinuity", this->GetEchoLevel() > 2 && rModelPart.GetCommunicator().MyPID() == 0) << "Finished setting up the dofs" << std::endl;
1501 
1502 #ifdef _OPENMP
1503  if (mlock_array.size() != 0)
1504  {
1505  for (int i = 0; i < static_cast<int>(mlock_array.size()); i++)
1506  omp_destroy_lock(&mlock_array[i]);
1507  }
1508 
1509  mlock_array.resize(BaseType::mDofSet.size());
1510 
1511  for (int i = 0; i < static_cast<int>(mlock_array.size()); i++)
1512  omp_init_lock(&mlock_array[i]);
1513 #endif
1514 
1515  // If reactions are to be calculated, we check if all the dofs have reactions defined
1516  // This is tobe done only in debug mode
1517 #ifdef KRATOS_DEBUG
1518  if (BaseType::GetCalculateReactionsFlag())
1519  {
1520  for (auto dof_iterator = BaseType::mDofSet.begin(); dof_iterator != BaseType::mDofSet.end(); ++dof_iterator)
1521  {
1522  KRATOS_ERROR_IF_NOT(dof_iterator->HasReaction()) << "Reaction variable not set for the following : " << std::endl
1523  << "Node : " << dof_iterator->Id() << std::endl
1524  << "Dof : " << (*dof_iterator) << std::endl
1525  << "Not possible to calculate reactions." << std::endl;
1526  }
1527  }
1528 #endif
1529 
1530  KRATOS_CATCH("");
1531  }
1532 
1537  void SetUpSystem(
1538  ModelPart &rModelPart) override
1539  {
1540  // Set equation id for degrees of freedom
1541  // the free degrees of freedom are positioned at the beginning of the system,
1542  // while the fixed one are at the end (in opposite order).
1543  //
1544  // that means that if the EquationId is greater than "mEquationSystemSize"
1545  // the pointed degree of freedom is restrained
1546  //
1547  int free_id = 0;
1548  int fix_id = BaseType::mDofSet.size();
1549 
1550  for (typename DofsArrayType::iterator dof_iterator = BaseType::mDofSet.begin(); dof_iterator != BaseType::mDofSet.end(); ++dof_iterator)
1551  if (dof_iterator->IsFixed())
1552  dof_iterator->SetEquationId(--fix_id);
1553  else
1554  dof_iterator->SetEquationId(free_id++);
1555 
1556  BaseType::mEquationSystemSize = fix_id;
1557  }
1558 
1559  //**************************************************************************
1560  //**************************************************************************
1561 
1562  void ResizeAndInitializeVectors(
1563  typename TSchemeType::Pointer pScheme,
1564  TSystemMatrixPointerType &pA,
1565  TSystemVectorPointerType &pDx,
1566  TSystemVectorPointerType &pb,
1567  ModelPart &rModelPart) override
1568  {
1569  KRATOS_TRY
1570 
1571  /* boost::timer c_contruct_matrix; */
1572 
1573  if (pA == NULL) // if the pointer is not initialized initialize it to an empty matrix
1574  {
1575  TSystemMatrixPointerType pNewA = TSystemMatrixPointerType(new TSystemMatrixType(0, 0));
1576  pA.swap(pNewA);
1577  }
1578  if (pDx == NULL) // if the pointer is not initialized initialize it to an empty matrix
1579  {
1580  TSystemVectorPointerType pNewDx = TSystemVectorPointerType(new TSystemVectorType(0));
1581  pDx.swap(pNewDx);
1582  }
1583  if (pb == NULL) // if the pointer is not initialized initialize it to an empty matrix
1584  {
1585  TSystemVectorPointerType pNewb = TSystemVectorPointerType(new TSystemVectorType(0));
1586  pb.swap(pNewb);
1587  }
1588  if (BaseType::mpReactionsVector == NULL) // if the pointer is not initialized initialize it to an empty matrix
1589  {
1590  TSystemVectorPointerType pNewReactionsVector = TSystemVectorPointerType(new TSystemVectorType(0));
1591  BaseType::mpReactionsVector.swap(pNewReactionsVector);
1592  }
1593 
1594  TSystemMatrixType &A = *pA;
1595  TSystemVectorType &Dx = *pDx;
1596  TSystemVectorType &b = *pb;
1597 
1598  // resizing the system vectors and matrix
1599  if (A.size1() == 0 || BaseType::GetReshapeMatrixFlag() == true) // if the matrix is not initialized
1600  {
1601  A.resize(BaseType::mEquationSystemSize, BaseType::mEquationSystemSize, false);
1602  ConstructMatrixStructure(pScheme, A, rModelPart);
1603  }
1604  else
1605  {
1606  if (A.size1() != BaseType::mEquationSystemSize || A.size2() != BaseType::mEquationSystemSize)
1607  {
1608  KRATOS_WATCH("it should not come here!!!!!!!! ... this is SLOW");
1609  KRATOS_ERROR << "The equation system size has changed during the simulation. This is not permited." << std::endl;
1610  A.resize(BaseType::mEquationSystemSize, BaseType::mEquationSystemSize, true);
1611  ConstructMatrixStructure(pScheme, A, rModelPart);
1612  }
1613  }
1614  if (Dx.size() != BaseType::mEquationSystemSize)
1615  Dx.resize(BaseType::mEquationSystemSize, false);
1616  if (b.size() != BaseType::mEquationSystemSize)
1617  b.resize(BaseType::mEquationSystemSize, false);
1618 
1619  // if needed resize the vector for the calculation of reactions
1620  if (BaseType::mCalculateReactionsFlag == true)
1621  {
1622  unsigned int ReactionsVectorSize = BaseType::mDofSet.size();
1623  if (BaseType::mpReactionsVector->size() != ReactionsVectorSize)
1624  BaseType::mpReactionsVector->resize(ReactionsVectorSize, false);
1625  }
1626  /* std::cout << "CONTINUITY EQ: contruct_matrix : " << c_contruct_matrix.elapsed() << std::endl; */
1627 
1628  KRATOS_CATCH("")
1629  }
1630 
1631  //**************************************************************************
1632  //**************************************************************************
1633 
1645  void ApplyDirichletConditions(
1646  typename TSchemeType::Pointer pScheme,
1647  ModelPart &rModelPart,
1648  TSystemMatrixType &A,
1649  TSystemVectorType &Dx,
1650  TSystemVectorType &b) override
1651  {
1652  }
1653 
1657  void Clear() override
1658  {
1659  this->mDofSet = DofsArrayType();
1660 
1661  if (this->mpReactionsVector != NULL)
1662  TSparseSpace::Clear((this->mpReactionsVector));
1663  // this->mReactionsVector = TSystemVectorType();
1664 
1665  this->mpLinearSystemSolver->Clear();
1666 
1667  KRATOS_INFO_IF("NodalResidualBasedEliminationBuilderAndSolverContinuity", this->GetEchoLevel() > 1) << "Clear Function called" << std::endl;
1668  }
1669 
1677  int Check(ModelPart &rModelPart) override
1678  {
1679  KRATOS_TRY
1680 
1681  return 0;
1682  KRATOS_CATCH("");
1683  }
1684 
1688 
1692 
1696 
1698 
1699  protected:
1702 
1706 
1710 
1714 
1715  void Assemble(
1716  TSystemMatrixType &A,
1717  TSystemVectorType &b,
1718  const LocalSystemMatrixType &LHS_Contribution,
1719  const LocalSystemVectorType &RHS_Contribution,
1720  const Element::EquationIdVectorType &EquationId
1721 #ifdef _OPENMP
1722  ,
1723  std::vector<omp_lock_t> &lock_array
1724 #endif
1725  )
1726  {
1727  unsigned int local_size = LHS_Contribution.size1();
1728 
1729  for (unsigned int i_local = 0; i_local < local_size; i_local++)
1730  {
1731  unsigned int i_global = EquationId[i_local];
1732 
1733  if (i_global < BaseType::mEquationSystemSize)
1734  {
1735 #ifdef _OPENMP
1736  omp_set_lock(&lock_array[i_global]);
1737 #endif
1738  b[i_global] += RHS_Contribution(i_local);
1739  for (unsigned int j_local = 0; j_local < local_size; j_local++)
1740  {
1741  unsigned int j_global = EquationId[j_local];
1742  if (j_global < BaseType::mEquationSystemSize)
1743  {
1744  A(i_global, j_global) += LHS_Contribution(i_local, j_local);
1745  }
1746  }
1747 #ifdef _OPENMP
1748  omp_unset_lock(&lock_array[i_global]);
1749 #endif
1750  }
1751  // note that assembly on fixed rows is not performed here
1752  }
1753  }
1754 
1755  //**************************************************************************
1756 
1757  virtual void ConstructMatrixStructure(
1758  typename TSchemeType::Pointer pScheme,
1759  TSystemMatrixType &A,
1760  ModelPart &rModelPart)
1761  {
1762  // filling with zero the matrix (creating the structure)
1763  Timer::Start("MatrixStructure");
1764 
1765  const std::size_t equation_size = BaseType::mEquationSystemSize;
1766 
1767  std::vector<std::unordered_set<std::size_t>> indices(equation_size);
1768 
1769 #pragma omp parallel for firstprivate(equation_size)
1770  for (int iii = 0; iii < static_cast<int>(equation_size); iii++)
1771  {
1772  indices[iii].reserve(40);
1773  }
1774 
1775  Element::EquationIdVectorType ids(3, 0);
1776 
1777 #pragma omp parallel firstprivate(ids)
1778  {
1779  // The process info
1780  ProcessInfo &r_current_process_info = rModelPart.GetProcessInfo();
1781 
1782  // We repeat the same declaration for each thead
1783  std::vector<std::unordered_set<std::size_t>> temp_indexes(equation_size);
1784 
1785 #pragma omp for
1786  for (int index = 0; index < static_cast<int>(equation_size); ++index)
1787  temp_indexes[index].reserve(30);
1788 
1789  // Getting the size of the array of elements from the model
1790  const int number_of_elements = static_cast<int>(rModelPart.Elements().size());
1791 
1792  // Element initial iterator
1793  const auto el_begin = rModelPart.ElementsBegin();
1794 
1795 // We iterate over the elements
1796 #pragma omp for schedule(guided, 512) nowait
1797  for (int i_elem = 0; i_elem < number_of_elements; ++i_elem)
1798  {
1799  auto it_elem = el_begin + i_elem;
1800  pScheme->EquationId(*it_elem, ids, r_current_process_info);
1801 
1802  for (auto &id_i : ids)
1803  {
1804  if (id_i < BaseType::mEquationSystemSize)
1805  {
1806  auto &row_indices = temp_indexes[id_i];
1807  for (auto &id_j : ids)
1808  if (id_j < BaseType::mEquationSystemSize)
1809  row_indices.insert(id_j);
1810  }
1811  }
1812  }
1813 
1814  // Getting the size of the array of the conditions
1815  const int number_of_conditions = static_cast<int>(rModelPart.Conditions().size());
1816 
1817  // Condition initial iterator
1818  const auto cond_begin = rModelPart.ConditionsBegin();
1819 
1820 // We iterate over the conditions
1821 #pragma omp for schedule(guided, 512) nowait
1822  for (int i_cond = 0; i_cond < number_of_conditions; ++i_cond)
1823  {
1824  auto it_cond = cond_begin + i_cond;
1825  pScheme->EquationId(*it_cond, ids, r_current_process_info);
1826  for (auto &id_i : ids)
1827  {
1828  if (id_i < BaseType::mEquationSystemSize)
1829  {
1830  auto &row_indices = temp_indexes[id_i];
1831  for (auto &id_j : ids)
1832  if (id_j < BaseType::mEquationSystemSize)
1833  row_indices.insert(id_j);
1834  }
1835  }
1836  }
1837 
1838 // Merging all the temporal indexes
1839 #pragma omp critical
1840  {
1841  for (int i = 0; i < static_cast<int>(temp_indexes.size()); ++i)
1842  {
1843  indices[i].insert(temp_indexes[i].begin(), temp_indexes[i].end());
1844  }
1845  }
1846  }
1847 
1848  // count the row sizes
1849  unsigned int nnz = 0;
1850  for (unsigned int i = 0; i < indices.size(); i++)
1851  nnz += indices[i].size();
1852 
1853  A = boost::numeric::ublas::compressed_matrix<double>(indices.size(), indices.size(), nnz);
1854 
1855  double *Avalues = A.value_data().begin();
1856  std::size_t *Arow_indices = A.index1_data().begin();
1857  std::size_t *Acol_indices = A.index2_data().begin();
1858 
1859  // filling the index1 vector - DO NOT MAKE PARALLEL THE FOLLOWING LOOP!
1860  Arow_indices[0] = 0;
1861  for (int i = 0; i < static_cast<int>(A.size1()); i++)
1862  Arow_indices[i + 1] = Arow_indices[i] + indices[i].size();
1863 
1864 #pragma omp parallel for
1865  for (int i = 0; i < static_cast<int>(A.size1()); i++)
1866  {
1867  const unsigned int row_begin = Arow_indices[i];
1868  const unsigned int row_end = Arow_indices[i + 1];
1869  unsigned int k = row_begin;
1870  for (auto it = indices[i].begin(); it != indices[i].end(); it++)
1871  {
1872  Acol_indices[k] = *it;
1873  Avalues[k] = 0.0;
1874  k++;
1875  }
1876 
1877  std::sort(&Acol_indices[row_begin], &Acol_indices[row_end]);
1878  }
1879 
1880  A.set_filled(indices.size() + 1, nnz);
1881 
1882  Timer::Stop("MatrixStructure");
1883  }
1884 
1885  void AssembleLHS(
1886  TSystemMatrixType &A,
1887  LocalSystemMatrixType &LHS_Contribution,
1888  Element::EquationIdVectorType &EquationId)
1889  {
1890  unsigned int local_size = LHS_Contribution.size1();
1891 
1892  for (unsigned int i_local = 0; i_local < local_size; i_local++)
1893  {
1894  unsigned int i_global = EquationId[i_local];
1895  if (i_global < BaseType::mEquationSystemSize)
1896  {
1897  for (unsigned int j_local = 0; j_local < local_size; j_local++)
1898  {
1899  unsigned int j_global = EquationId[j_local];
1900  if (j_global < BaseType::mEquationSystemSize)
1901  A(i_global, j_global) += LHS_Contribution(i_local, j_local);
1902  }
1903  }
1904  }
1905  }
1906 
1910 
1914 
1918 
1920 
1921  private:
1924 
1928 #ifdef _OPENMP
1929  std::vector<omp_lock_t> mlock_array;
1930 #endif
1934 
1938 
1939  inline void AddUnique(std::vector<std::size_t> &v, const std::size_t &candidate)
1940  {
1941  std::vector<std::size_t>::iterator i = v.begin();
1942  std::vector<std::size_t>::iterator endit = v.end();
1943  while (i != endit && (*i) != candidate)
1944  {
1945  i++;
1946  }
1947  if (i == endit)
1948  {
1949  v.push_back(candidate);
1950  }
1951  }
1952  inline void CreatePartition(unsigned int number_of_threads, const int number_of_rows, DenseVector<unsigned int> &partitions)
1953  {
1954  partitions.resize(number_of_threads + 1);
1955  int partition_size = number_of_rows / number_of_threads;
1956  partitions[0] = 0;
1957  partitions[number_of_threads] = number_of_rows;
1958  for (unsigned int i = 1; i < number_of_threads; i++)
1959  partitions[i] = partitions[i - 1] + partition_size;
1960  }
1961 
1962  void AssembleRHS(
1963  TSystemVectorType &b,
1964  const LocalSystemVectorType &RHS_Contribution,
1965  const Element::EquationIdVectorType &EquationId)
1966  {
1967  unsigned int local_size = RHS_Contribution.size();
1968 
1969  if (BaseType::mCalculateReactionsFlag == false)
1970  {
1971  for (unsigned int i_local = 0; i_local < local_size; i_local++)
1972  {
1973  const unsigned int i_global = EquationId[i_local];
1974 
1975  if (i_global < BaseType::mEquationSystemSize) // free dof
1976  {
1977  // ASSEMBLING THE SYSTEM VECTOR
1978  double &b_value = b[i_global];
1979  const double &rhs_value = RHS_Contribution[i_local];
1980 
1981 #pragma omp atomic
1982  b_value += rhs_value;
1983  }
1984  }
1985  }
1986  else
1987  {
1988  TSystemVectorType &ReactionsVector = *BaseType::mpReactionsVector;
1989  for (unsigned int i_local = 0; i_local < local_size; i_local++)
1990  {
1991  const unsigned int i_global = EquationId[i_local];
1992 
1993  if (i_global < BaseType::mEquationSystemSize) // free dof
1994  {
1995  // ASSEMBLING THE SYSTEM VECTOR
1996  double &b_value = b[i_global];
1997  const double &rhs_value = RHS_Contribution[i_local];
1998 
1999 #pragma omp atomic
2000  b_value += rhs_value;
2001  }
2002  else // fixed dof
2003  {
2004  double &b_value = ReactionsVector[i_global - BaseType::mEquationSystemSize];
2005  const double &rhs_value = RHS_Contribution[i_local];
2006 
2007 #pragma omp atomic
2008  b_value += rhs_value;
2009  }
2010  }
2011  }
2012  }
2013 
2014  //**************************************************************************
2015 
2016  void AssembleLHS_CompleteOnFreeRows(
2017  TSystemMatrixType &A,
2018  LocalSystemMatrixType &LHS_Contribution,
2019  Element::EquationIdVectorType &EquationId)
2020  {
2021  unsigned int local_size = LHS_Contribution.size1();
2022  for (unsigned int i_local = 0; i_local < local_size; i_local++)
2023  {
2024  unsigned int i_global = EquationId[i_local];
2025  if (i_global < BaseType::mEquationSystemSize)
2026  {
2027  for (unsigned int j_local = 0; j_local < local_size; j_local++)
2028  {
2029  int j_global = EquationId[j_local];
2030  A(i_global, j_global) += LHS_Contribution(i_local, j_local);
2031  }
2032  }
2033  }
2034  }
2035 
2039 
2043 
2047 
2051 
2053 
2054  }; /* Class NodalResidualBasedEliminationBuilderAndSolverContinuity */
2055 
2057 
2060 
2062 
2063 } /* namespace Kratos.*/
2064 
2065 #endif /* KRATOS_NODAL_RESIDUAL_BASED_ELIMINATION_BUILDER_AND_SOLVER defined */
builder_and_solver.h
Kratos::BuilderAndSolver
Current class provides an implementation for the base builder and solving operations.
Definition: builder_and_solver.h:64
Kratos::BuilderAndSolver::IndexType
std::size_t IndexType
Definition of the index type.
Definition: builder_and_solver.h:76
Kratos::BuilderAndSolver::mDofSetIsInitialized
bool mDofSetIsInitialized
If the matrix is reshaped each step.
Definition: builder_and_solver.h:743
Kratos::BuilderAndSolver::TSystemVectorType
TSparseSpace::VectorType TSystemVectorType
Definition of the vector size.
Definition: builder_and_solver.h:85
Kratos::BuilderAndSolver::NodesArrayType
ModelPart::NodesContainerType NodesArrayType
The containers of the entities.
Definition: builder_and_solver.h:109
Kratos::BuilderAndSolver::GetCalculateReactionsFlag
bool GetCalculateReactionsFlag() const
This method returns the flag mCalculateReactionsFlag.
Definition: builder_and_solver.h:184
Kratos::BuilderAndSolver::TSystemMatrixType
TSparseSpace::MatrixType TSystemMatrixType
Definition of the sparse matrix.
Definition: builder_and_solver.h:82
Kratos::BuilderAndSolver::ConditionsArrayType
ModelPart::ConditionsContainerType ConditionsArrayType
Definition: builder_and_solver.h:111
Kratos::BuilderAndSolver::mpReactionsVector
TSystemVectorPointerType mpReactionsVector
Definition: builder_and_solver.h:751
Kratos::BuilderAndSolver::GetReshapeMatrixFlag
bool GetReshapeMatrixFlag() const
This method returns the flag mReshapeMatrixFlag.
Definition: builder_and_solver.h:220
Kratos::BuilderAndSolver::LocalSystemMatrixType
TDenseSpace::MatrixType LocalSystemMatrixType
The local matrix definition.
Definition: builder_and_solver.h:94
Kratos::BuilderAndSolver::TSystemVectorPointerType
TSparseSpace::VectorPointerType TSystemVectorPointerType
Definition of the pointer to the vector.
Definition: builder_and_solver.h:91
Kratos::BuilderAndSolver::TDataType
TSparseSpace::DataType TDataType
Definition of the data type.
Definition: builder_and_solver.h:79
Kratos::BuilderAndSolver::mpLinearSystemSolver
TLinearSolver::Pointer mpLinearSystemSolver
Definition: builder_and_solver.h:737
Kratos::BuilderAndSolver::mDofSet
DofsArrayType mDofSet
Pointer to the linear solver.
Definition: builder_and_solver.h:739
Kratos::BuilderAndSolver::LocalSystemVectorType
TDenseSpace::VectorType LocalSystemVectorType
The local vector definition.
Definition: builder_and_solver.h:97
Kratos::BuilderAndSolver::TSystemMatrixPointerType
TSparseSpace::MatrixPointerType TSystemMatrixPointerType
Definition of the pointer to the sparse matrix.
Definition: builder_and_solver.h:88
Kratos::BuilderAndSolver::mCalculateReactionsFlag
bool mCalculateReactionsFlag
Flag taking care if the dof set was initialized ot not.
Definition: builder_and_solver.h:745
Kratos::BuilderAndSolver::GetEchoLevel
int GetEchoLevel() const
It returns the echo level.
Definition: builder_and_solver.h:674
Kratos::BuilderAndSolver::mEquationSystemSize
unsigned int mEquationSystemSize
Flag taking in account if it is needed or not to calculate the reactions.
Definition: builder_and_solver.h:747
Kratos::BuilderAndSolver::ElementsArrayType
ModelPart::ElementsContainerType ElementsArrayType
Definition: builder_and_solver.h:110
Kratos::Communicator::MyPID
virtual int MyPID() const
Definition: communicator.cpp:91
Kratos::Dof::Pointer
Dof * Pointer
Pointer definition of Dof.
Definition: dof.h:93
Kratos::Element::DofsVectorType
std::vector< DofType::Pointer > DofsVectorType
Definition: element.h:100
Kratos::Element::EquationIdVectorType
std::vector< std::size_t > EquationIdVectorType
Definition: element.h:98
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::size
SizeType size() const
Definition: geometry.h:518
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::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::ElementsBegin
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
Kratos::ModelPart::ConditionsBegin
ConditionIterator ConditionsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1361
Kratos::ModelPart::GetCommunicator
Communicator & GetCommunicator()
Definition: model_part.h:1821
Kratos::ModelPart::Conditions
ConditionsContainerType & Conditions(IndexType ThisIndex=0)
Definition: model_part.h:1381
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::NodeIterator
MeshType::NodeIterator NodeIterator
Definition: model_part.h:134
Kratos::ModelPart::Elements
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
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::NodalResidualBasedEliminationBuilderAndSolverContinuity::LocalSystemMatrixType
BaseType::LocalSystemMatrixType LocalSystemMatrixType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:100
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::ElementsArrayType
BaseType::ElementsArrayType ElementsArrayType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:108
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::SetUpSystem
void SetUpSystem(ModelPart &rModelPart) override
Organises the dofset in order to speed up the building phase.
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1537
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::NodeType
Node NodeType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:105
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::SystemSolveWithPhysics
void SystemSolveWithPhysics(TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b, ModelPart &rModelPart)
This is a call to the linear system solver (taking into account some physical particularities of the ...
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1195
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::AssembleLHS
void AssembleLHS(TSystemMatrixType &A, LocalSystemMatrixType &LHS_Contribution, Element::EquationIdVectorType &EquationId)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1885
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::TSystemMatrixPointerType
BaseType::TSystemMatrixPointerType TSystemMatrixPointerType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:102
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::NodalResidualBasedEliminationBuilderAndSolverContinuity
NodalResidualBasedEliminationBuilderAndSolverContinuity(typename TLinearSolver::Pointer pNewLinearSystemSolver)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:121
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::TSystemMatrixType
BaseType::TSystemMatrixType TSystemMatrixType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:94
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BuildNodallyUnlessLaplacian
void BuildNodallyUnlessLaplacian(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &b)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:530
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::TSystemVectorType
BaseType::TSystemVectorType TSystemVectorType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:96
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::VectorType
Vector VectorType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:113
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::GetDeviatoricCoefficientForFluid
void GetDeviatoricCoefficientForFluid(ModelPart &rModelPart, ModelPart::NodeIterator itNode, double &deviatoricCoefficient)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:443
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BuildNodallyNotStabilized
void BuildNodallyNotStabilized(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &b)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:918
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BuildAndSolve
void BuildAndSolve(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
Function to perform the building and solving phase at the same time.
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1240
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::SetUpDofSet
void SetUpDofSet(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart) override
Builds the list of the DofSets involved in the problem by "asking" to each element and condition its ...
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1407
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::ConstructMatrixStructure
virtual void ConstructMatrixStructure(typename TSchemeType::Pointer pScheme, TSystemMatrixType &A, ModelPart &rModelPart)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1757
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::ConditionsArrayType
BaseType::ConditionsArrayType ConditionsArrayType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:109
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::Check
int Check(ModelPart &rModelPart) override
This function is designed to be called once to perform all the checks needed on the input provided....
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1677
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::Clear
void Clear() override
This function is intended to be called at the end of the solution step to clean up memory storage not...
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1657
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BuildNodally
void BuildNodally(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &b)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:142
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::TDataType
BaseType::TDataType TDataType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:90
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BuildAll
void BuildAll(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &b)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1001
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BuildNodallyNoVolumetricStabilizedTerms
void BuildNodallyNoVolumetricStabilizedTerms(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &b)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:764
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(NodalResidualBasedEliminationBuilderAndSolverContinuity)
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::~NodalResidualBasedEliminationBuilderAndSolverContinuity
~NodalResidualBasedEliminationBuilderAndSolverContinuity() override
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:130
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::TSystemVectorPointerType
BaseType::TSystemVectorPointerType TSystemVectorPointerType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:103
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::NodesArrayType
BaseType::NodesArrayType NodesArrayType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:107
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::Assemble
void Assemble(TSystemMatrixType &A, TSystemVectorType &b, const LocalSystemMatrixType &LHS_Contribution, const LocalSystemVectorType &RHS_Contribution, const Element::EquationIdVectorType &EquationId)
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1715
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::TSchemeType
BaseType::TSchemeType TSchemeType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:88
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::ApplyDirichletConditions
void ApplyDirichletConditions(typename TSchemeType::Pointer pScheme, ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
Applies the dirichlet conditions. This operation may be very heavy or completely unexpensive dependin...
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1645
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::ElementsContainerType
BaseType::ElementsContainerType ElementsContainerType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:111
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::Build
void Build(typename TSchemeType::Pointer pScheme, ModelPart &r_model_part, TSystemMatrixType &A, TSystemVectorType &b) override
equivalent (but generally faster) then performing BuildLHS and BuildRHS
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1304
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::BaseType
BuilderAndSolver< TSparseSpace, TDenseSpace, TLinearSolver > BaseType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:86
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::DofsArrayType
BaseType::DofsArrayType DofsArrayType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:92
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::LocalSystemVectorType
BaseType::LocalSystemVectorType LocalSystemVectorType
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:98
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::ResizeAndInitializeVectors
void ResizeAndInitializeVectors(typename TSchemeType::Pointer pScheme, TSystemMatrixPointerType &pA, TSystemVectorPointerType &pDx, TSystemVectorPointerType &pb, ModelPart &rModelPart) override
This method initializes and resizes the system of equations.
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1562
Kratos::NodalResidualBasedEliminationBuilderAndSolverContinuity::SystemSolve
void SystemSolve(TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b) override
This is a call to the linear system solver.
Definition: nodal_residualbased_elimination_builder_and_solver_continuity.h:1161
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::ParallelUtilities::GetNumThreads
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
Kratos::PointerVectorSet
A sorted associative container similar to an STL set, but uses a vector to store pointers to its data...
Definition: pointer_vector_set.h:72
Kratos::PointerVectorSet::iterator
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: pointer_vector_set.h:95
Kratos::PointerVectorSet::size
size_type size() const
Returns the number of elements in the container.
Definition: pointer_vector_set.h:502
Kratos::PointerVectorSet::Sort
void Sort()
Sort the elements in the set.
Definition: pointer_vector_set.h:753
Kratos::PointerVectorSet::push_back
void push_back(TPointerType x)
Adds a pointer to the end of the set.
Definition: pointer_vector_set.h:544
Kratos::PointerVectorSet::reserve
void reserve(int reservedsize)
Reserves memory for a specified number of elements.
Definition: pointer_vector_set.h:733
Kratos::PointerVectorSet::end
iterator end()
Returns an iterator pointing to the end of the container.
Definition: pointer_vector_set.h:314
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::Scheme
This class provides the implementation of the basic tasks that are needed by the solution strategy.
Definition: scheme.h:56
Kratos::Timer::Start
static void Start(std::string const &rIntervalName)
This method starts the timer meassures.
Definition: timer.cpp:109
Kratos::Timer::Stop
static void Stop(std::string const &rIntervalName)
This method stops the timer meassures.
Definition: timer.cpp:125
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_WATCH
#define KRATOS_WATCH(variable)
Definition: define.h:806
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
KRATOS_ERROR
#define KRATOS_ERROR
Definition: exception.h:161
KRATOS_ERROR_IF_NOT
#define KRATOS_ERROR_IF_NOT(conditional)
Definition: exception.h:163
KRATOS_ERROR_IF
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
KRATOS_INFO_IF
#define KRATOS_INFO_IF(label, conditional)
Definition: logger.h:251
KRATOS_WARNING_IF
#define KRATOS_WARNING_IF(label, conditional)
Definition: logger.h:266
key_hash.h
model_part.h
Kratos::Globals::Pi
constexpr double Pi
Definition: global_variables.h:25
Kratos::Python::TwoNorm
double TwoNorm(SparseSpaceType &dummy, SparseSpaceType::VectorType &x)
Definition: add_strategies_to_python.cpp:164
Kratos::Python::Size
TSpaceType::IndexType Size(TSpaceType &dummy, typename TSpaceType::VectorType const &rV)
Definition: add_strategies_to_python.cpp:111
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::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
face_heat.density
float density
Definition: face_heat.py:56
generate_convection_diffusion_explicit_element.v
v
Definition: generate_convection_diffusion_explicit_element.py:114
generate_total_lagrangian_mixed_volumetric_strain_element.local_size
int local_size
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:17
generate_total_lagrangian_mixed_volumetric_strain_element.b
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
isotropic_damage_automatic_differentiation.dimension
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
sensitivityMatrix.A
A
Definition: sensitivityMatrix.py:70
tensors::i
integer i
Definition: TensorModule.f:17
test_pureconvectionsolver_benchmarking.pressure
float pressure
Definition: test_pureconvectionsolver_benchmarking.py:101
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
pfem_fluid_dynamics_application_variables.h
Generated on Sat Feb 3 2024 02:08:58 for KratosMultiphysics by doxygen 1.9.1