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.
mapping_variables_2D_utilities.hpp
Go to the documentation of this file.
1 
2 // | / |
3 // ' / __| _` | __| _ \ __|
4 // . \ | ( | | ( |\__ `
5 // _|\_\_| \__,_|\__|\___/ ____/
6 // Multi-Physics
7 //
8 // License: BSD License
9 // Kratos default license: kratos/license.txt
10 //
11 // Main authors: Joaquin Irazabal
12 //
13 
14 
15 #if !defined(KRATOS_MAPPING_VARIABLES_2D_UTILITIES )
16 #define KRATOS_MAPPING_VARIABLES_2D_UTILITIES
17 
18 // System includes
19 #include <fstream>
20 #include <iostream>
21 #include <cmath>
22 
23 // Project includes
24 #include "geometries/geometry.h"
25 #include "includes/define.h"
26 #include "includes/model_part.h"
27 #include "includes/kratos_parameters.h"
28 #include "utilities/openmp_utils.h"
29 
30 // Application includes
31 #include "dam_application_variables.h"
32 
33 namespace Kratos
34 {
35 
36 class MappingVariables2DUtilities
37 {
38 
39 protected:
40 
42 
43  struct UtilityVariables
44  {
45  double X_max, X_min, Y_max, Y_min;
46  int NRows, NColumns;
47  double RowSize, ColumnSize;
48  };
49 
50 
51 public:
52 
53  KRATOS_CLASS_POINTER_DEFINITION( MappingVariables2DUtilities );
54 
56  MappingVariables2DUtilities() {}
57 
59 
61  virtual ~MappingVariables2DUtilities() {}
62 
64 
65  void MappingThermalModelParts (ModelPart& rModelPartOld, ModelPart& rModelPartNew, bool add_temperature, bool add_reference_temperature)
66  {
67  // Define necessary variables
68  UtilityVariables AuxVariables;
69  this->InitializeMapping(AuxVariables,rModelPartNew);
70  this->NodalThermalVariablesMapping(AuxVariables,rModelPartOld,rModelPartNew,add_temperature,add_reference_temperature);
71  }
72 
73  void MappingMechanicalModelParts (ModelPart& rModelPartOld, ModelPart& rModelPartNew, bool add_displacement, bool add_stress)
74  {
75  // Define necessary variables
76  UtilityVariables AuxVariables;
77  this->InitializeMapping(AuxVariables,rModelPartNew);
78  this->NodalMechanicalVariablesMapping(AuxVariables,rModelPartOld,rModelPartNew,add_displacement,add_stress);
79  }
80 
81 protected:
82 
84 
86 
88 
89  void InitializeMapping(
90  UtilityVariables& rAuxVariables,
91  ModelPart& rModelPartNew)
92  {
93  this->ComputeCellMatrixDimensions(rAuxVariables,rModelPartNew);
94  }
95 
97 
98  void NodalThermalVariablesMapping(
99  const UtilityVariables& AuxVariables,
100  ModelPart& rModelPartOld,
101  ModelPart& rModelPartNew,
102  bool add_temperature,
103  bool add_reference_temperature)
104  {
105  // Define ElementOld Cell matrix
106  std::vector< std::vector< std::vector<Element::Pointer> > > ElementOldCellMatrix;
107  ElementOldCellMatrix.resize(AuxVariables.NRows);
108  for(int i = 0; i < AuxVariables.NRows; i++) ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
109 
110  // Locate Old Elments in cells
111  double X_me;
112  double Y_me;
113  int PointsNumber;
114 
115  int NElems = static_cast<int>(rModelPartOld.Elements().size());
116  ModelPart::ElementsContainerType::iterator el_begin = rModelPartOld.ElementsBegin();
117 
118  #pragma omp parallel for private(X_me,Y_me,PointsNumber)
119  for(int i = 0; i < NElems; i++)
120  {
121  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
122 
123  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
124  double X_right = X_left;
125  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
126  double Y_bot = Y_top;
127  PointsNumber = itElemOld->GetGeometry().PointsNumber();
128 
129  for(int j = 1; j < PointsNumber; j++)
130  {
131  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
132  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
133 
134  if(X_me > X_right) X_right = X_me;
135  else if(X_me < X_left) X_left = X_me;
136  if(Y_me > Y_top) Y_top = Y_me;
137  else if(Y_me < Y_bot) Y_bot = Y_me;
138  }
139 
140  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
141  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
142  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
143  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
144 
145  if(Column_left < 0) Column_left = 0;
146  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
147  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
148  else if(Column_right < 0) Column_right = 0;
149 
150  if(Row_top < 0) Row_top = 0;
151  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
152  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
153  else if(Row_bot < 0) Row_bot = 0;
154 
155  for(int k = Row_top; k <= Row_bot; k++)
156  {
157  for(int l = Column_left; l <= Column_right; l++)
158  {
159  #pragma omp critical
160  {
161  ElementOldCellMatrix[k][l].push_back((*(itElemOld.base())));
162  }
163  }
164  }
165 
166  for(int k = Row_top; k <= Row_bot; k++)
167  {
168  for(int l = Column_left; l <= Column_right; l++)
169  {
170  #pragma omp critical
171  {
172  ElementOldCellMatrix[k][l].push_back((*(itElemOld.base())));
173  }
174  }
175  }
176  }
177 
178  // Locate new nodes inside old elements and interpolate nodal variables
179  const int NNodes = static_cast<int>(rModelPartNew.Nodes().size());
180  ModelPart::NodesContainerType::iterator node_begin = rModelPartNew.NodesBegin();
181 
182  array_1d<double,3> GlobalCoordinates;
183  array_1d<double,3> LocalCoordinates;
184 
185  #pragma omp parallel for private(X_me,Y_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
186  for(int i = 0; i < NNodes; i++)
187  {
188  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
189 
190  X_me = itNodeNew->X0();
191  Y_me = itNodeNew->Y0();
192 
193  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
194  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
195 
196  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
197  else if(Column < 0) Column = 0;
198  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
199  else if(Row < 0) Row = 0;
200 
201  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
202  bool IsInside = false;
203  Element::Pointer pElementOld;
204 
205  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
206  {
207  pElementOld = ElementOldCellMatrix[Row][Column][m];
208  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
209  if(IsInside) break;
210  }
211 
212  if(IsInside==false)
213  {
214  for(int i = 0; i < 12; i++)
215  {
216  for(int j = 1; j<10; j++)
217  {
218  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
219  {
220  pElementOld = ElementOldCellMatrix[Row][Column][m];
221  double tol = j * pow(10.0, (i-12));
222  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,tol);
223  if(IsInside) break;
224  }
225  if(IsInside) break;
226  }
227  if(IsInside) break;
228  }
229  }
230 
231  if(IsInside == false)
232  {
233  std::cout << "None of the old elements contains node " << itNodeNew->Id() << std::endl;
234  continue;
235  }
236 
237  PointsNumber = pElementOld->GetGeometry().PointsNumber();
238  Vector ShapeFunctionsValuesVector(PointsNumber);
239  Vector NodalVariableVector(PointsNumber);
240 
241  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
242 
243  // Interpolation of nodal variables
244  if (add_temperature)
245  {
246  double max_temperature = -10e12;
247  double min_temperature = 10e12;
248 
249  for(int j = 0; j < PointsNumber; j++)
250  {
251  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(TEMPERATURE);
252  if (NodalVariableVector[j] > max_temperature) max_temperature = NodalVariableVector[j];
253  if (NodalVariableVector[j] < min_temperature) min_temperature = NodalVariableVector[j];
254  }
255  itNodeNew->FastGetSolutionStepValue(TEMPERATURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
256  if (itNodeNew->FastGetSolutionStepValue(TEMPERATURE) > max_temperature) itNodeNew->FastGetSolutionStepValue(TEMPERATURE) = max_temperature;
257  if (itNodeNew->FastGetSolutionStepValue(TEMPERATURE) < min_temperature) itNodeNew->FastGetSolutionStepValue(TEMPERATURE) = min_temperature;
258  }
259 
260  if (add_reference_temperature)
261  {
262  double max_reference_temperature = -10e12;
263  double min_reference_temperature = 10e12;
264 
265  for(int j = 0; j < PointsNumber; j++)
266  {
267  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE);
268  if (NodalVariableVector[j] > max_reference_temperature) max_reference_temperature = NodalVariableVector[j];
269  if (NodalVariableVector[j] < min_reference_temperature) min_reference_temperature = NodalVariableVector[j];
270  }
271  itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
272  if (itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) > max_reference_temperature) itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) = max_reference_temperature;
273  if (itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) < min_reference_temperature) itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) = min_reference_temperature;
274  }
275  }
276  }
277 
278  void NodalMechanicalVariablesMapping(
279  const UtilityVariables& AuxVariables,
280  ModelPart& rModelPartOld,
281  ModelPart& rModelPartNew,
282  bool add_displacement,
283  bool add_stress)
284  {
285  // Define ElementOld Cell matrix
286  std::vector< std::vector< std::vector<Element::Pointer> > > ElementOldCellMatrix;
287  ElementOldCellMatrix.resize(AuxVariables.NRows);
288  for(int i = 0; i < AuxVariables.NRows; i++) ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
289 
290  // Locate Old Elments in cells
291  double X_me;
292  double Y_me;
293  int PointsNumber;
294 
295  int NElems = static_cast<int>(rModelPartOld.Elements().size());
296  ModelPart::ElementsContainerType::iterator el_begin = rModelPartOld.ElementsBegin();
297 
298  #pragma omp parallel for private(X_me,Y_me,PointsNumber)
299  for(int i = 0; i < NElems; i++)
300  {
301  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
302 
303  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
304  double X_right = X_left;
305  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
306  double Y_bot = Y_top;
307  PointsNumber = itElemOld->GetGeometry().PointsNumber();
308 
309  for(int j = 1; j < PointsNumber; j++)
310  {
311  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
312  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
313 
314  if(X_me > X_right) X_right = X_me;
315  else if(X_me < X_left) X_left = X_me;
316  if(Y_me > Y_top) Y_top = Y_me;
317  else if(Y_me < Y_bot) Y_bot = Y_me;
318  }
319 
320  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
321  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
322  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
323  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
324 
325  if(Column_left < 0) Column_left = 0;
326  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
327  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
328  else if(Column_right < 0) Column_right = 0;
329 
330  if(Row_top < 0) Row_top = 0;
331  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
332  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
333  else if(Row_bot < 0) Row_bot = 0;
334 
335  for(int k = Row_top; k <= Row_bot; k++)
336  {
337  for(int l = Column_left; l <= Column_right; l++)
338  {
339  #pragma omp critical
340  {
341  ElementOldCellMatrix[k][l].push_back((*(itElemOld.base())));
342  }
343  }
344  }
345  }
346 
347  // Locate new nodes inside old elements and interpolate nodal variables
348  const int NNodes = static_cast<int>(rModelPartNew.Nodes().size());
349  ModelPart::NodesContainerType::iterator node_begin = rModelPartNew.NodesBegin();
350 
351  array_1d<double,3> GlobalCoordinates;
352  array_1d<double,3> LocalCoordinates;
353 
354  #pragma omp parallel for private(X_me,Y_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
355  for(int i = 0; i < NNodes; i++)
356  {
357  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
358 
359  X_me = itNodeNew->X0();
360  Y_me = itNodeNew->Y0();
361 
362  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
363  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
364 
365  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
366  else if(Column < 0) Column = 0;
367  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
368  else if(Row < 0) Row = 0;
369 
370  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
371  bool IsInside = false;
372  Element::Pointer pElementOld;
373 
374  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
375  {
376  pElementOld = ElementOldCellMatrix[Row][Column][m];
377  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
378  if(IsInside) break;
379  }
380 
381  if(IsInside==false)
382  {
383  for(int i = 0; i < 12; i++)
384  {
385  for(int j = 1; j<10; j++)
386  {
387  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
388  {
389  pElementOld = ElementOldCellMatrix[Row][Column][m];
390  double tol = j * pow(10.0, (i-12));
391  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,tol);
392  if(IsInside) break;
393  }
394  if(IsInside) break;
395  }
396  if(IsInside) break;
397  }
398  }
399 
400  if(IsInside == false)
401  {
402  std::cout << "None of the old elements contains node " << itNodeNew->Id() << std::endl;
403  if (add_stress)
404  {
405  itNodeNew->FastGetSolutionStepValue(NODAL_CAUCHY_STRESS_TENSOR) = ZeroMatrix(2,2);
406  itNodeNew->FastGetSolutionStepValue(INITIAL_NODAL_CAUCHY_STRESS_TENSOR) = ZeroMatrix(2,2);
407  }
408  continue;
409  }
410 
411  PointsNumber = pElementOld->GetGeometry().PointsNumber();
412  Vector ShapeFunctionsValuesVector(PointsNumber);
413  Vector NodalVariableVector(PointsNumber);
414 
415  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
416 
417  // Interpolation of nodal variables
418  if (add_displacement)
419  {
420  double max_displacement_x = -10e12;
421  double min_displacement_x = 10e12;
422  for(int j = 0; j < PointsNumber; j++)
423  {
424  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_X);
425  if (NodalVariableVector[j] > max_displacement_x) max_displacement_x = NodalVariableVector[j];
426  if (NodalVariableVector[j] < min_displacement_x) min_displacement_x = NodalVariableVector[j];
427  }
428  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
429  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) > max_displacement_x) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = max_displacement_x;
430  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) < min_displacement_x) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = min_displacement_x;
431 
432  double max_displacement_y = -10e12;
433  double min_displacement_y = 10e12;
434  for(int j = 0; j < PointsNumber; j++)
435  {
436  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_Y);
437  if (NodalVariableVector[j] > max_displacement_y) max_displacement_y = NodalVariableVector[j];
438  if (NodalVariableVector[j] < min_displacement_y) min_displacement_y = NodalVariableVector[j];
439  }
440  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
441  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) > max_displacement_y) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = max_displacement_y;
442  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) < min_displacement_y) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = min_displacement_y;
443  }
444 
445  if (add_stress)
446  {
447  Matrix NodalStress = ZeroMatrix(2,2);
448 
449  Matrix MaxNodalStress = ZeroMatrix(2,2);
450  Matrix MinNodalStress = ZeroMatrix(2,2);
451 
452  for(int k = 0; k < 2; k++)
453  {
454  for(int l = 0; l < 2; l++)
455  {
456  MaxNodalStress(k,l) = -10e12;
457  MinNodalStress(k,l) = 10e12;
458  for(int j = 0; j < PointsNumber; j++)
459  {
460  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(NODAL_CAUCHY_STRESS_TENSOR)(k,l);
461  if (NodalVariableVector[j] > MaxNodalStress(k,l)) MaxNodalStress(k,l) = NodalVariableVector[j];
462  if (NodalVariableVector[j] < MinNodalStress(k,l)) MinNodalStress(k,l) = NodalVariableVector[j];
463  }
464  NodalStress(k,l) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
465  if (NodalStress(k,l) > MaxNodalStress(k,l)) NodalStress(k,l) = MaxNodalStress(k,l);
466  if (NodalStress(k,l) < MinNodalStress(k,l)) NodalStress(k,l) = MinNodalStress(k,l);
467  }
468  }
469 
470  itNodeNew->FastGetSolutionStepValue(INITIAL_NODAL_CAUCHY_STRESS_TENSOR) = NodalStress;
471  itNodeNew->FastGetSolutionStepValue(NODAL_CAUCHY_STRESS_TENSOR) = NodalStress;
472  }
473  }
474  }
475 
476 private:
477 
479 
480  void ComputeCellMatrixDimensions(
481  UtilityVariables& rAuxVariables,
482  ModelPart& rModelPart)
483  {
484  // Compute X and Y limits of the current geometry
485  unsigned int NumThreads = ParallelUtilities::GetNumThreads();
486  std::vector<double> X_max_partition(NumThreads);
487  std::vector<double> X_min_partition(NumThreads);
488  std::vector<double> Y_max_partition(NumThreads);
489  std::vector<double> Y_min_partition(NumThreads);
490 
491  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
492  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
493 
494  #pragma omp parallel
495  {
496  int k = OpenMPUtils::ThisThread();
497  X_max_partition[k] = node_begin->X();
498  X_min_partition[k] = X_max_partition[k];
499  Y_max_partition[k] = node_begin->Y();
500  Y_min_partition[k] = Y_max_partition[k];
501 
502  #pragma omp for
503  for(int i = 0; i < NNodes; i++)
504  {
505  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
506 
507  double X_me = itNode->X();
508  double Y_me = itNode->Y();
509 
510  if( X_me > X_max_partition[k] ) X_max_partition[k] = X_me;
511  else if( X_me < X_min_partition[k] ) X_min_partition[k] = X_me;
512 
513  if( Y_me > Y_max_partition[k] ) Y_max_partition[k] = Y_me;
514  else if( Y_me < Y_min_partition[k] ) Y_min_partition[k] = Y_me;
515  }
516  }
517 
518  rAuxVariables.X_max = X_max_partition[0];
519  rAuxVariables.X_min = X_min_partition[0];
520  rAuxVariables.Y_max = Y_max_partition[0];
521  rAuxVariables.Y_min = Y_min_partition[0];
522 
523  for(unsigned int i=1; i < NumThreads; i++)
524  {
525  if(X_max_partition[i] > rAuxVariables.X_max) rAuxVariables.X_max = X_max_partition[i];
526  if(X_min_partition[i] < rAuxVariables.X_min) rAuxVariables.X_min = X_min_partition[i];
527  if(Y_max_partition[i] > rAuxVariables.Y_max) rAuxVariables.Y_max = Y_max_partition[i];
528  if(Y_min_partition[i] < rAuxVariables.Y_min) rAuxVariables.Y_min = Y_min_partition[i];
529  }
530 
531  // Calculate Average Element Length
532  double AverageElementLength = 0.0;
533 
534  int NElems = static_cast<int>(rModelPart.Elements().size());
535  ModelPart::ElementsContainerType::iterator el_begin = rModelPart.ElementsBegin();
536 
537  #pragma omp parallel for reduction(+:AverageElementLength)
538  for(int i = 0; i < NElems; i++)
539  {
540  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
541 
542  AverageElementLength += itElem->GetGeometry().Length();
543  }
544 
545  AverageElementLength = AverageElementLength/NElems;
546 
547  // Compute FracturePoints CellMatrix dimensions
548  rAuxVariables.NRows = int((rAuxVariables.Y_max-rAuxVariables.Y_min)/(AverageElementLength));
549  rAuxVariables.NColumns = int((rAuxVariables.X_max-rAuxVariables.X_min)/(AverageElementLength));
550  if(rAuxVariables.NRows <= 0) rAuxVariables.NRows = 1;
551  if(rAuxVariables.NColumns <= 0) rAuxVariables.NColumns = 1;
552  rAuxVariables.RowSize = (rAuxVariables.Y_max-rAuxVariables.Y_min)/rAuxVariables.NRows;
553  rAuxVariables.ColumnSize = (rAuxVariables.X_max-rAuxVariables.X_min)/rAuxVariables.NColumns;
554  }
555 
556 }; // Class MappingVariables2DUtilities
557 
558 } // namespace Kratos.
559 
560 #endif /* KRATOS_MAPPING_VARIABLES_2D_UTILITIES defined */
Kratos::MappingVariables2DUtilities
Definition: mapping_variables_2D_utilities.hpp:37
Kratos::MappingVariables2DUtilities::InitializeMapping
void InitializeMapping(UtilityVariables &rAuxVariables, ModelPart &rModelPartNew)
Member Variables.
Definition: mapping_variables_2D_utilities.hpp:89
Kratos::MappingVariables2DUtilities::NodalMechanicalVariablesMapping
void NodalMechanicalVariablesMapping(const UtilityVariables &AuxVariables, ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_displacement, bool add_stress)
Definition: mapping_variables_2D_utilities.hpp:278
Kratos::MappingVariables2DUtilities::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(MappingVariables2DUtilities)
Kratos::MappingVariables2DUtilities::NodalThermalVariablesMapping
void NodalThermalVariablesMapping(const UtilityVariables &AuxVariables, ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_temperature, bool add_reference_temperature)
Definition: mapping_variables_2D_utilities.hpp:98
Kratos::MappingVariables2DUtilities::MappingVariables2DUtilities
MappingVariables2DUtilities()
Constructor.
Definition: mapping_variables_2D_utilities.hpp:56
Kratos::MappingVariables2DUtilities::~MappingVariables2DUtilities
virtual ~MappingVariables2DUtilities()
Destructor.
Definition: mapping_variables_2D_utilities.hpp:61
Kratos::MappingVariables2DUtilities::MappingThermalModelParts
void MappingThermalModelParts(ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_temperature, bool add_reference_temperature)
Definition: mapping_variables_2D_utilities.hpp:65
Kratos::MappingVariables2DUtilities::MappingMechanicalModelParts
void MappingMechanicalModelParts(ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_displacement, bool add_stress)
Definition: mapping_variables_2D_utilities.hpp:73
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::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
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::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::ParallelUtilities::GetNumThreads
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
Kratos::array_1d< double, 3 >
dam_application_variables.h
kratos_parameters.h
model_part.h
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::inner_prod
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::int
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
hinsberg_optimization.tol
int tol
Definition: hinsberg_optimization.py:138
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
run_marine_rain_substepping.m
int m
Definition: run_marine_rain_substepping.py:8
tensors::i
integer i
Definition: TensorModule.f:17
tensors::l
integer l
Definition: TensorModule.f:17
openmp_utils.h
Kratos::MappingVariables2DUtilities::UtilityVariables
Basic Structs for the utility -----------------------------------------------------------------------...
Definition: mapping_variables_2D_utilities.hpp:44
Kratos::MappingVariables2DUtilities::UtilityVariables::NColumns
int NColumns
Definition: mapping_variables_2D_utilities.hpp:46
Kratos::MappingVariables2DUtilities::UtilityVariables::Y_max
double Y_max
Definition: mapping_variables_2D_utilities.hpp:45
Kratos::MappingVariables2DUtilities::UtilityVariables::X_max
double X_max
Definition: mapping_variables_2D_utilities.hpp:45
Kratos::MappingVariables2DUtilities::UtilityVariables::X_min
double X_min
Definition: mapping_variables_2D_utilities.hpp:45
Kratos::MappingVariables2DUtilities::UtilityVariables::ColumnSize
double ColumnSize
Definition: mapping_variables_2D_utilities.hpp:47
Kratos::MappingVariables2DUtilities::UtilityVariables::Y_min
double Y_min
Definition: mapping_variables_2D_utilities.hpp:45
Kratos::MappingVariables2DUtilities::UtilityVariables::RowSize
double RowSize
Definition: mapping_variables_2D_utilities.hpp:47
Kratos::MappingVariables2DUtilities::UtilityVariables::NRows
int NRows
Definition: mapping_variables_2D_utilities.hpp:46
Generated on Sat Feb 3 2024 01:59:56 for KratosMultiphysics by doxygen 1.9.1