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.
calculate_water_fraction.h
Go to the documentation of this file.
1 #if !defined(KRATOS_CALCULATE_WATER_FRACTION_UTILITY_INCLUDED )
2 #define KRATOS_CALCULATE_WATER_FRACTION_UTILITY_INCLUDED
3 
4 // System includes
5 #include <string>
6 #include <iostream>
7 #include <algorithm>
8 
9 // Project includes
10 #include "includes/define.h"
12 #include "utilities/math_utils.h"
13 #include "utilities/geometry_utilities.h"
15 #include "includes/variables.h"
16 #include "includes/model_part.h"
17 #include "includes/node.h"
18 #include "includes/element.h"
20 
21 
22 namespace Kratos
23 {
24  template< unsigned int TDim>
26  {
27  public:
28 
30 
32  : mr_model_part(model_part) //mr_model_part is saved as private variable (declared at the end of the file)
33  {
35  //std::cout << "Hello, I am the constructor of the Utility" << std::endl;
36  KRATOS_CATCH("")
37  }
38 
39 
41  {}
42  /*
43  double Calculate() //water fraction
44  {
45  KRATOS_TRY
46  //double area; //we create the needed variables
47  double sum_area=0.0;
48  double sum_water_area=0.0;
49  //double one_third=1.0/3.0;
50  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
51 
52  #pragma omp parallel for reduction(+:sum_area) reduction(+:sum_water_area)
53  for(unsigned int ii=0; ii<mr_model_part.Nodes().size(); ii++)
54  {
55  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
56  const double & nodal_area = inode->FastGetSolutionStepValue(NODAL_AREA); //resetting the temperature
57  sum_area += nodal_area;
58  if ((inode->FastGetSolutionStepValue(DISTANCE))<0.0)
59  sum_water_area += nodal_area;
60  }
61  const double water_fraction = sum_water_area / sum_area;
62  //std::cout << "Finished, the mean temperature is" << water_fraction << std::endl; //we print the result
63  return water_fraction;
64 
65  KRATOS_CATCH("")
66  }
67  */
68  double Calculate()
69  {
71 
72 
73  double sum_areas=1.0e-100;
74  //double sum_temperatures=0.0;
75  //double nodal_weight=1.0/(1.0+double(TDim));
76 
77  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
78 
79 
80  vector<unsigned int> element_partition;
81  #ifdef _OPENMP
82  int number_of_threads = omp_get_max_threads();
83  #else
84  int number_of_threads = 1;
85  #endif
86  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
87 
88  #pragma omp parallel for reduction(+:sum_areas)
89  for(int kkk=0; kkk<number_of_threads; kkk++)
90  {
91  double thread_sum_areas=0.0;
92  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
93  {
94  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
95  double Area;
96  BoundedMatrix<double, (TDim+1), TDim > DN_DX;
97  array_1d<double, (TDim+1) > N;
98  Geometry<Node >& geom = ielem->GetGeometry();
99  GeometryUtils::CalculateGeometryData(geom, DN_DX, N, Area);
100  //sum_areas+=area;
101  int negative_nodes=0;
102  int positive_nodes=0;
103  for (unsigned int k = 0; k < (TDim+1); k++)
104  {
105  if(geom[k].FastGetSolutionStepValue(DISTANCE)<0.0)
106  negative_nodes++;
107  else
108  positive_nodes++;
109  }
110  if (negative_nodes==(TDim+1))
111  thread_sum_areas+=Area;
112 
113  else if (negative_nodes>0)
114  {
115  array_1d<double,(TDim+1)> distances;
116  for (unsigned int i = 0; i < (TDim+1); i++)
117  {
118  distances[i] = geom[i].FastGetSolutionStepValue(DISTANCE);
119  }
120 
121 
122  BoundedMatrix<double,3*(TDim-1), 2> Nenriched;
123  array_1d<double,(3*(TDim-1))> volumes;
124  BoundedMatrix<double,(TDim+1), TDim > coords;
125  BoundedMatrix<double, 3*(TDim-1), (TDim+1) > Ngauss;
126  array_1d<double,(3*(TDim-1))> signs;
127  std::vector< Matrix > gauss_gradients(3*(TDim-1));
128  //fill coordinates
129 
130  //unsigned int single_triangle_node = 1;
131  for (unsigned int i = 0; i < (TDim+1); i++)
132  {
133  const array_1d<double, 3 > & xyz = geom[i].Coordinates();
134  for (unsigned int j = 0; j < TDim; j++)
135  coords(i,j)=xyz(j);
136  }
137  for (unsigned int i = 0; i < 3*(TDim-1); i++)
138  gauss_gradients[i].resize(2, TDim, false); //2 values of the 2 shape functions, and derivates in (xy) direction).
139  unsigned int ndivisions = EnrichmentUtilities::CalculateEnrichedShapeFuncions(coords, DN_DX, distances, volumes, Ngauss, signs, gauss_gradients, Nenriched);
140  for (unsigned int i=0;i!=ndivisions;i++)
141  if (signs(i)<0.0) thread_sum_areas+=volumes(i);
142 
143 
144  }
145  }
146  sum_areas = thread_sum_areas;
147 
148  }
149  return sum_areas;
150 
151  KRATOS_CATCH("")
152  }
153 
154  double CalculateWaterHeight(double x_position)
155  {
156  KRATOS_TRY
157 
158  double all_threads_water_height=-100000000.0;
159  const double tolerance=0.001;
160  const double upper_limit=x_position+tolerance;
161  const double lower_limit=x_position-tolerance;
162  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
163 
164  vector<unsigned int> node_partition;
165  #ifdef _OPENMP
166  int number_of_threads = omp_get_max_threads();
167  #else
168  int number_of_threads = 1;
169  #endif
170  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
171 
172  #pragma omp parallel for
173  for(int kkk=0; kkk<number_of_threads; kkk++)
174  {
175  double local_thread_water_height=-100000000.0;
176  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
177  {
178  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
179  double & distance = (inode->FastGetSolutionStepValue(DISTANCE));
180  if ( distance <0.0)
181  {
182  if ((inode->X())<upper_limit && (inode->X())>lower_limit && (inode->Y())>local_thread_water_height)
183  local_thread_water_height = (inode->Y());
184  }
185  //now we search for the node given certain criteria
186 
187  }
188 
189 
190  if ( local_thread_water_height > all_threads_water_height )
191  {
192  #pragma omp critical
193  {
194  if ( local_thread_water_height > all_threads_water_height ) all_threads_water_height = local_thread_water_height;
195  }
196  }
197  }
198 
199  return all_threads_water_height;
200 
201  KRATOS_CATCH("")
202  }
203 
204 
206  {
207  KRATOS_TRY
208 
209  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
210  const double delta_t = CurrentProcessInfo[DELTA_TIME];
211 
212  double all_threads_max_courant = 0.0;
213  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
214 
215  vector<unsigned int> element_partition;
216  #ifdef _OPENMP
217  int number_of_threads = omp_get_max_threads();
218  #else
219  int number_of_threads = 1;
220  #endif
221  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
222 
223  #pragma omp parallel for
224  for(int kkk=0; kkk<number_of_threads; kkk++)
225  {
226  double local_thread_max_courant = 0.0;
227  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
228  {
229  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
230  //double & distance = (inode->FastGetSolutionStepValue(DISTANCE));
231  if ((ielem->GetValue(VELOCITY_OVER_ELEM_SIZE))>local_thread_max_courant)
232  local_thread_max_courant = (ielem->GetValue(VELOCITY_OVER_ELEM_SIZE));
233  }
234 
235 
236  if ( local_thread_max_courant > all_threads_max_courant )
237  {
238  #pragma omp critical
239  {
240  if ( local_thread_max_courant > all_threads_max_courant ) all_threads_max_courant = local_thread_max_courant;
241  }
242  }
243  }
244 
245  all_threads_max_courant *= delta_t * 1.414;
246 
247  return all_threads_max_courant;
248 
249  KRATOS_CATCH("")
250  }
251 
253  {
254  KRATOS_TRY
255 
256  //using a nodal approach (faster!)
257  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
258  const double delta_t = CurrentProcessInfo[DELTA_TIME];
259 
260  double all_threads_max_courant = 0.0;
261  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
262 
263  vector<unsigned int> node_partition;
264  #ifdef _OPENMP
265  int number_of_threads = omp_get_max_threads();
266  #else
267  int number_of_threads = 1;
268  #endif
269  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
270 
271  #pragma omp parallel for
272  for(int kkk=0; kkk<number_of_threads; kkk++)
273  {
274  double local_thread_max_courant = 0.0;
275  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
276  {
277  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
278  const double & distance = (inode->FastGetSolutionStepValue(DISTANCE));
279  const double velocity = sqrt(pow(inode->FastGetSolutionStepValue(VELOCITY_X),2)+pow(inode->FastGetSolutionStepValue(VELOCITY_Y),2)+pow(inode->FastGetSolutionStepValue(VELOCITY_Z),2));
280  const double nodal_courant = (velocity*delta_t/inode->FastGetSolutionStepValue(MEAN_SIZE));
281  if(nodal_courant>local_thread_max_courant && distance < 0.0) //only for negative nodes!
282  local_thread_max_courant = nodal_courant;
283  }
284 
285 
286  if ( local_thread_max_courant > all_threads_max_courant )
287  {
288  #pragma omp critical
289  {
290  if ( local_thread_max_courant > all_threads_max_courant ) all_threads_max_courant = local_thread_max_courant;
291  }
292  }
293  }
294 
295  //all_threads_max_courant *= delta_t * 1.414;
296 
297  return all_threads_max_courant;
298 
299  KRATOS_CATCH("")
300  }
301 
302  double CalculateMeanCourant() //water fraction
303  {
304  KRATOS_TRY
305 
306  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
307  const double delta_t = CurrentProcessInfo[DELTA_TIME];
308 
309  //double area=0.0; //we create the needed variables
310  //double number_of_threads = double(OpenMPUtils::GetNumThreads());
311  double sum_courant=0.0;
312 
313  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
314 
315  vector<unsigned int> element_partition;
316  #ifdef _OPENMP
317  int number_of_threads = omp_get_max_threads();
318  #else
319  int number_of_threads = 1;
320  #endif
321  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
322 
323  #pragma omp parallel for reduction(+:sum_courant)
324  for(int kkk=0; kkk<number_of_threads; kkk++)
325  {
326  double thread_sum_courant=0.0;
327  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
328  {
329  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
330  //const double & velocity_over_elem_size = (ielem->GetValue(VELOCITY_OVER_ELEM_SIZE));
331  if ((ielem->GetValue(VELOCITY_OVER_ELEM_SIZE))>0.0)
332  thread_sum_courant += ielem->GetValue(VELOCITY_OVER_ELEM_SIZE);
333  }
334  sum_courant += thread_sum_courant;
335  }
336 
337  sum_courant *= delta_t * 1.414 / double(mr_model_part.Elements().size());
338 
339  return sum_courant;
340 
341  KRATOS_CATCH("")
342  }
343 
344  //NOW ONLY VISCOUS. but since in the first step we cannot use the pressure we just add the viscoust forces. still, lines to use pressure can be uncommented
345  double CalculateForce(int direction) //
346  {
347  KRATOS_TRY
348 
349  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
350  double viscosity = CurrentProcessInfo[VISCOSITY];
351  //double delta_t = CurrentProcessInfo[DELTA_TIME];
352  //array_1d<double,3> & gravity= CurrentProcessInfo[GRAVITY];
353 
354 
355  const array_1d<double,3> zero3 = ZeroVector(3);
356  double nodal_weight = 1.0/ (double (TDim) );
357 
358  double force=0.0;
359 
360  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
361 
362 
363  vector<unsigned int> element_partition;
364  #ifdef _OPENMP
365  int number_of_threads = omp_get_max_threads();
366  #else
367  int number_of_threads = 1;
368  #endif
369  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
370 
371  #pragma omp parallel for reduction(+:force)
372  for(int kkk=0; kkk<number_of_threads; kkk++)
373  {
374  double thread_force=0.0;
375  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
376  {
377  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
378  if (ielem->Is(ACTIVE)) //elements can be inactive to add temporary walls. fractional velocity is integrated by parts so walls are seen as having zero velocity without doing anything
379  {
380  //double Area;
381  Geometry<Node >& geom = ielem->GetGeometry();
382 
383  array_1d<unsigned int, 4 > fixed_nodes; //unordered : i position in the array might not correspond to i node of the element
384  array_1d<bool, 4 > is_node_fixed; //i position belongs to i node of the element
385  unsigned int number_of_fixed_nodes=0;
386  //bool boundary_element=false;
387  BoundedMatrix<double, 4, 3 > velocities = ZeroMatrix(4, 3);
388 
389  for (unsigned int i = 0; i < (TDim+1); i++)
390  {
391  const array_1d<double, 3 > & velocity = geom[i].FastGetSolutionStepValue(VELOCITY);
392  for (unsigned int j = 0; j < (TDim); j++)
393  velocities(i,j) = velocity[j];
394 
395  if constexpr (TDim==2)
396  {
397  if (geom[i].IsFixed(FRACT_VEL_X) && geom[i].IsFixed(FRACT_VEL_Y))
398  {
399  fixed_nodes[number_of_fixed_nodes]=i;
400  is_node_fixed[i]=true;
401  number_of_fixed_nodes++;
402  }
403  else
404  is_node_fixed[i]=false;
405  }
406  else // (TDim==3)
407  {
408  if (geom[i].IsFixed(FRACT_VEL_X) && geom[i].IsFixed(FRACT_VEL_Y) && geom[i].IsFixed(FRACT_VEL_Z) )
409  {
410  fixed_nodes[number_of_fixed_nodes]=i;
411  number_of_fixed_nodes++;
412  is_node_fixed[i]=true;
413  }
414  else
415  is_node_fixed[i]=false;
416  }
417  }
418 
419  double plane_point_distance=1.0;
420  double fixed_face_area_or_lenght=0.0;
421  array_1d<double, 3 > boundary_stress;
422  if (number_of_fixed_nodes==TDim) //it means we have cutted elements!
423  {
424  //boundary_element=true;
425  array_1d<double, 3 > normal;
426  unsigned int free_node=0;
427  if constexpr (TDim==2)
428  {
429  fixed_face_area_or_lenght = fabs(sqrt(pow((geom[fixed_nodes[1]].Y()-geom[fixed_nodes[0]].Y()),2 ) + pow( (geom[fixed_nodes[1]].X()-geom[fixed_nodes[0]].X() ),2 ) ) );
430  normal[0] = geom[fixed_nodes[1]].Y()-geom[fixed_nodes[0]].Y();
431  normal[1] = - ( geom[fixed_nodes[1]].X()-geom[fixed_nodes[0]].X() );
432  normal[2] = 0.0;
433  normal /= sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
434 
435  if (fixed_nodes[0]==0)
436  {
437  if (fixed_nodes[1]==1)
438  free_node=2;
439  else
440  free_node=1;
441  }
442  else
443  free_node=0;
444 
445  //the plane is composed by the unit normal and any of the points of fixed nodes. we will use fixed_nodes[0];
446  plane_point_distance = inner_prod( (geom[free_node].Coordinates()-geom[fixed_nodes[0]].Coordinates()) , normal);
447  //boundary_stress = geom[free_node].FastGetSolutionStepValue(VELOCITY)*viscosity/plane_point_distance;
448  if (plane_point_distance<0.0)
449  {
450  plane_point_distance*=-1.0;
451  normal *= -1.0;
452  }
453  }
454  else //(TDim==3)
455  {
456  //the area is obtained from the crossproduct of the 2 vertices:
457  MathUtils<double>::CrossProduct(normal, geom[fixed_nodes[1]].Coordinates() - geom[fixed_nodes[0]].Coordinates(), geom[fixed_nodes[2]].Coordinates() - geom[fixed_nodes[0]].Coordinates() );
458  fixed_face_area_or_lenght = 0.5 * sqrt( pow(normal[0],2) + pow(normal[1],2) + pow(normal[2],2) );
459  normal /= 2.0 * fixed_face_area_or_lenght; //this way it is a unit vector. now we must find the distance from the plane generated by the triangles to the free node:
460 
461  //fixed_face_area_or_lenght = fabs(fixed_face_area_or_lenght);
462 
463  for (unsigned int j=0; j!=(TDim+1); j++)
464  {
465  if (is_node_fixed[j]==false)
466  {
467  free_node=j;
468  break;
469  }
470  }
471 
472  //the plane is composed by the unit normal and any of the points of fixed nodes. we will use fixed_nodes[0];
473  plane_point_distance = inner_prod( (geom[free_node].Coordinates()-geom[fixed_nodes[0]].Coordinates()) , normal);
474  if (plane_point_distance<0.0)
475  normal *= -1.0;
476  {
477  plane_point_distance*=-1.0;
478  normal *= -1.0;
479  }
480  //boundary_stress = geom[free_node].FastGetSolutionStepValue(VELOCITY)*viscosity/plane_point_distance;
481  }
482 
483  boundary_stress = - geom[free_node].FastGetSolutionStepValue(VELOCITY)*viscosity/(fabs(plane_point_distance));
484  //KRATOS_WATCH(plane_point_distance)
485  //KRATOS_WATCH(boundary_stress)
486  //KRATOS_WATCH(fixed_face_area_or_lenght)
487 
488  //drag forces:
489  thread_force += boundary_stress[direction]*fixed_face_area_or_lenght; // unit density! careful!
490  //face_force+=fixed_face_area_or_lenght*normal[direction];
491  //now pressure forces:
492 
493  for (unsigned int j=0; j!=(TDim); j++) // the 2 or 3 nodes that define the fixed face:
494  {
495  /*
496  if ( (geom[fixed_nodes[j]].X())<5.0 )
497  face_force += nodal_weight*(0.5*fixed_face_area_or_lenght*normal[direction]);
498  else
499  face_force -= nodal_weight*(0.5*fixed_face_area_or_lenght*normal[direction]);
500  */
501  thread_force +=nodal_weight*(geom[fixed_nodes[j]].FastGetSolutionStepValue(PRESSURE))*fixed_face_area_or_lenght*normal[direction];
502 
503  //array_1d<double,3> & nodal_normal= (geom[fixed_nodes[j]].FastGetSolutionStepValue(NORMAL));
504  //face_force += (geom[fixed_nodes[j]].FastGetSolutionStepValue(PRESSURE))*nodal_normal[direction];
505  }
506 
507 
508 
509 
510  }
511  }
512 
513  }
514 
515  force+=thread_force;
516  }
517 
518  return force;
519 
520  KRATOS_CATCH("")
521  }
522 
523 
524 
525  protected:
526 
527 
528  private:
529 
530  ModelPart& mr_model_part;
531 
532  };
533 
534 } // namespace Kratos.
535 
536 #endif // KRATOS_CALCULATE__WATER_FRACTION_UTILITY_INCLUDED defined
Definition: calculate_water_fraction.h:26
~CalculateWaterFraction()
Definition: calculate_water_fraction.h:40
double CalculateForce(int direction)
Definition: calculate_water_fraction.h:345
KRATOS_CLASS_POINTER_DEFINITION(CalculateWaterFraction)
double CalculateMaxCourant()
Definition: calculate_water_fraction.h:205
double Calculate()
Definition: calculate_water_fraction.h:68
double CalculateMeanCourant()
Definition: calculate_water_fraction.h:302
double CalculateWaterHeight(double x_position)
Definition: calculate_water_fraction.h:154
CalculateWaterFraction(ModelPart &model_part)
Definition: calculate_water_fraction.h:31
double CalculateMaxCourantInNegativeElements()
Definition: calculate_water_fraction.h:252
static int CalculateEnrichedShapeFuncions(BoundedMatrix< double, 4, 3 > &rPoints, BoundedMatrix< double, 4, 3 > &DN_DX, array_1d< double, 4 > &rDistances, array_1d< double, 6 > &rVolumes, BoundedMatrix< double, 6, 4 > &rShapeFunctionValues, array_1d< double, 6 > &rPartitionsSign, std::vector< Matrix > &rGradientsValue, BoundedMatrix< double, 6, 2 > &NEnriched)
Definition: enrichment_utilities.h:500
Geometry base class.
Definition: geometry.h:71
static void CalculateGeometryData(const GeometryType &rGeometry, BoundedMatrix< double, 4, 3 > &rDN_DX, array_1d< double, 4 > &rN, double &rVolume)
This function is designed to compute the shape function derivatives, shape functions and volume in 3D...
Definition: geometry_utilities.h:176
Definition: amatrix_interface.h:41
static T CrossProduct(const T &a, const T &b)
Performs the vector product of the two input vectors a,b.
Definition: math_utils.h:762
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Short class definition.
Definition: array_1d.h:61
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_TRY
Definition: define.h:109
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
TABLE_NUMBER_ANGULAR_VELOCITY TABLE_NUMBER_MOMENT I33 BEAM_INERTIA_ROT_UNIT_LENGHT_Y KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BEAM_INERTIA_ROT_UNIT_LENGHT_Z) typedef std double
Definition: DEM_application_variables.h:182
float velocity
Definition: PecletTest.py:54
float viscosity
Definition: edgebased_var.py:8
model_part
Definition: face_heat.py:14
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
float delta_t
Definition: rotatingcone_PureConvectionBenchmarking.py:129
N
Definition: sensitivityMatrix.py:29
integer i
Definition: TensorModule.f:17