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_3D_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_3D_UTILITIES )
16 #define KRATOS_MAPPING_VARIABLES_3D_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 MappingVariables3DUtilities
37 {
38 
39 protected:
40 
42 
43  struct UtilityVariables
44  {
45  double X_max, X_min, Y_max, Y_min, Z_max, Z_min;
46  int NRows, NColumns, NSections;
47  double RowSize, ColumnSize, SectionSize;
48  };
49 
50 public:
51 
52  KRATOS_CLASS_POINTER_DEFINITION( MappingVariables3DUtilities );
53 
55  MappingVariables3DUtilities() {}
56 
58 
60  virtual ~MappingVariables3DUtilities() {}
61 
63 
64  void MappingThermalModelParts (ModelPart& rModelPartOld, ModelPart& rModelPartNew, bool add_temperature, bool add_reference_temperature)
65  {
66  // Define necessary variables
67  UtilityVariables AuxVariables;
68  this->InitializeMapping(AuxVariables,rModelPartNew);
69  this->NodalThermalVariablesMapping(AuxVariables,rModelPartOld,rModelPartNew,add_temperature,add_reference_temperature);
70  }
71 
72  void MappingMechanicalModelParts (ModelPart& rModelPartOld, ModelPart& rModelPartNew, bool add_displacement, bool add_stress)
73  {
74  // Define necessary variables
75  UtilityVariables AuxVariables;
76  this->InitializeMapping(AuxVariables,rModelPartNew);
77  this->NodalMechanicalVariablesMapping(AuxVariables,rModelPartOld,rModelPartNew,add_displacement,add_stress);
78  }
79 
80 protected:
81 
83 
85 
87 
88  void InitializeMapping(
89  UtilityVariables& rAuxVariables,
90  ModelPart& rModelPartNew)
91  {
92  this->ComputeCellMatrixDimensions(rAuxVariables,rModelPartNew);
93  }
94 
96 
97  void NodalThermalVariablesMapping(
98  const UtilityVariables& AuxVariables,
99  ModelPart& rModelPartOld,
100  ModelPart& rModelPartNew,
101  bool add_temperature,
102  bool add_reference_temperature)
103  {
104  // Define ElementOld Cell matrix
105  std::vector< std::vector< std::vector< std::vector<Element::Pointer> > > > ElementOldCellMatrix;
106  ElementOldCellMatrix.resize(AuxVariables.NRows);
107  for(int i = 0; i < AuxVariables.NRows; i++)
108  ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
109  for(int i = 0; i < AuxVariables.NRows; i++)
110  {
111  for(int j = 0; j < AuxVariables.NColumns; j++)
112  {
113  ElementOldCellMatrix[i][j].resize(AuxVariables.NSections);
114  }
115  }
116 
117  // Locate Old Elments in cells
118  double X_me;
119  double Y_me;
120  double Z_me;
121  int PointsNumber;
122 
123  int NElems = static_cast<int>(rModelPartOld.Elements().size());
124  ModelPart::ElementsContainerType::iterator el_begin = rModelPartOld.ElementsBegin();
125 
126  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber)
127  for(int i = 0; i < NElems; i++)
128  {
129  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
130 
131  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
132  double X_right = X_left;
133  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
134  double Y_bot = Y_top;
135  double Z_back = itElemOld->GetGeometry().GetPoint(0).Z0();
136  double Z_front = Z_back;
137  PointsNumber = itElemOld->GetGeometry().PointsNumber();
138 
139  for(int j = 1; j < PointsNumber; j++)
140  {
141  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
142  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
143  Z_me = itElemOld->GetGeometry().GetPoint(j).Z0();
144 
145  if(X_me > X_right) X_right = X_me;
146  else if(X_me < X_left) X_left = X_me;
147  if(Y_me > Y_top) Y_top = Y_me;
148  else if(Y_me < Y_bot) Y_bot = Y_me;
149  if(Z_me > Z_front) Z_front = Z_me;
150  else if(Z_me < Z_back) Z_back = Z_me;
151  }
152 
153  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
154  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
155  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
156  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
157  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
158  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
159 
160  if(Column_left < 0) Column_left = 0;
161  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
162  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
163  else if(Column_right < 0) Column_right = 0;
164  if(Row_top < 0) Row_top = 0;
165  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
166  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
167  else if(Row_bot < 0) Row_bot = 0;
168  if(Section_back < 0) Section_back = 0;
169  else if(Section_back >= AuxVariables.NSections) Section_back = AuxVariables.NSections-1;
170  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
171  else if(Section_front < 0) Section_front = 0;
172 
173  for(int s = Section_back; s <= Section_front; s++)
174  {
175  for(int k = Row_top; k <= Row_bot; k++)
176  {
177  for(int l = Column_left; l <= Column_right; l++)
178  {
179  #pragma omp critical
180  {
181  ElementOldCellMatrix[k][l][s].push_back((*(itElemOld.base())));
182  }
183  }
184  }
185  }
186  }
187 
188  // Locate new nodes inside old elements and interpolate nodal variables
189  const int NNodes = static_cast<int>(rModelPartNew.Nodes().size());
190  ModelPart::NodesContainerType::iterator node_begin = rModelPartNew.NodesBegin();
191 
192  array_1d<double,3> GlobalCoordinates;
193  array_1d<double,3> LocalCoordinates;
194 
195  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
196  for(int i = 0; i < NNodes; i++)
197  {
198  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
199 
200  X_me = itNodeNew->X0();
201  Y_me = itNodeNew->Y0();
202  Z_me = itNodeNew->Z0();
203 
204  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
205  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
206  int Section = int((Z_me-AuxVariables.Z_min)/AuxVariables.SectionSize);
207 
208  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
209  else if(Column < 0) Column = 0;
210  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
211  else if(Row < 0) Row = 0;
212  if(Section >= AuxVariables.NSections) Section = AuxVariables.NSections-1;
213  else if(Section < 0) Section = 0;
214 
215  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
216  bool IsInside = false;
217  Element::Pointer pElementOld;
218 
219  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column][Section]).size(); m++)
220  {
221  pElementOld = ElementOldCellMatrix[Row][Column][Section][m];
222  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
223  if(IsInside) break;
224  }
225 
226  if(IsInside==false)
227  {
228  for(int i = 0; i < 12; i++)
229  {
230  for(int j = 1; j<10; j++)
231  {
232  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column][Section]).size(); m++)
233  {
234  pElementOld = ElementOldCellMatrix[Row][Column][Section][m];
235  double tol = j * pow(10.0, (i-12));
236  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,tol);
237  if(IsInside) break;
238  }
239  if(IsInside) break;
240  }
241  if(IsInside) break;
242  }
243  }
244 
245  if(IsInside == false)
246  {
247  std::cout << "None of the old elements contains node " << itNodeNew->Id() << std::endl;
248  continue;
249  }
250 
251  PointsNumber = pElementOld->GetGeometry().PointsNumber();
252  Vector ShapeFunctionsValuesVector(PointsNumber);
253  Vector NodalVariableVector(PointsNumber);
254 
255  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
256 
257  // Interpolation of nodal variables
258  if (add_temperature)
259  {
260  double max_temperature = -10e12;
261  double min_temperature = 10e12;
262 
263  for(int j = 0; j < PointsNumber; j++)
264  {
265  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(TEMPERATURE);
266  if (NodalVariableVector[j] > max_temperature) max_temperature = NodalVariableVector[j];
267  if (NodalVariableVector[j] < min_temperature) min_temperature = NodalVariableVector[j];
268  }
269  itNodeNew->FastGetSolutionStepValue(TEMPERATURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
270  if (itNodeNew->FastGetSolutionStepValue(TEMPERATURE) > max_temperature) itNodeNew->FastGetSolutionStepValue(TEMPERATURE) = max_temperature;
271  if (itNodeNew->FastGetSolutionStepValue(TEMPERATURE) < min_temperature) itNodeNew->FastGetSolutionStepValue(TEMPERATURE) = min_temperature;
272  }
273 
274  if (add_reference_temperature)
275  {
276  double max_reference_temperature = -10e12;
277  double min_reference_temperature = 10e12;
278 
279  for(int j = 0; j < PointsNumber; j++)
280  {
281  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE);
282  if (NodalVariableVector[j] > max_reference_temperature) max_reference_temperature = NodalVariableVector[j];
283  if (NodalVariableVector[j] < min_reference_temperature) min_reference_temperature = NodalVariableVector[j];
284  }
285  itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
286  if (itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) > max_reference_temperature) itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) = max_reference_temperature;
287  if (itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) < min_reference_temperature) itNodeNew->FastGetSolutionStepValue(NODAL_REFERENCE_TEMPERATURE) = min_reference_temperature;
288  }
289  }
290  }
291 
292  void NodalMechanicalVariablesMapping(
293  const UtilityVariables& AuxVariables,
294  ModelPart& rModelPartOld,
295  ModelPart& rModelPartNew,
296  bool add_displacement,
297  bool add_stress)
298  {
299  // Define ElementOld Cell matrix
300  std::vector< std::vector< std::vector< std::vector<Element::Pointer> > > > ElementOldCellMatrix;
301  ElementOldCellMatrix.resize(AuxVariables.NRows);
302  for(int i = 0; i < AuxVariables.NRows; i++)
303  ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
304  for(int i = 0; i < AuxVariables.NRows; i++)
305  {
306  for(int j = 0; j < AuxVariables.NColumns; j++)
307  {
308  ElementOldCellMatrix[i][j].resize(AuxVariables.NSections);
309  }
310  }
311 
312  // Locate Old Elments in cells
313  double X_me;
314  double Y_me;
315  double Z_me;
316  int PointsNumber;
317 
318  int NElems = static_cast<int>(rModelPartOld.Elements().size());
319  ModelPart::ElementsContainerType::iterator el_begin = rModelPartOld.ElementsBegin();
320 
321  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber)
322  for(int i = 0; i < NElems; i++)
323  {
324  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
325 
326  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
327  double X_right = X_left;
328  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
329  double Y_bot = Y_top;
330  double Z_back = itElemOld->GetGeometry().GetPoint(0).Z0();
331  double Z_front = Z_back;
332  PointsNumber = itElemOld->GetGeometry().PointsNumber();
333 
334  for(int j = 1; j < PointsNumber; j++)
335  {
336  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
337  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
338  Z_me = itElemOld->GetGeometry().GetPoint(j).Z0();
339 
340  if(X_me > X_right) X_right = X_me;
341  else if(X_me < X_left) X_left = X_me;
342  if(Y_me > Y_top) Y_top = Y_me;
343  else if(Y_me < Y_bot) Y_bot = Y_me;
344  if(Z_me > Z_front) Z_front = Z_me;
345  else if(Z_me < Z_back) Z_back = Z_me;
346  }
347 
348  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
349  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
350  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
351  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
352  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
353  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
354 
355  if(Column_left < 0) Column_left = 0;
356  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
357  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
358  else if(Column_right < 0) Column_right = 0;
359  if(Row_top < 0) Row_top = 0;
360  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
361  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
362  else if(Row_bot < 0) Row_bot = 0;
363  if(Section_back < 0) Section_back = 0;
364  else if(Section_back >= AuxVariables.NSections) Section_back = AuxVariables.NSections-1;
365  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
366  else if(Section_front < 0) Section_front = 0;
367 
368  for(int s = Section_back; s <= Section_front; s++)
369  {
370  for(int k = Row_top; k <= Row_bot; k++)
371  {
372  for(int l = Column_left; l <= Column_right; l++)
373  {
374  #pragma omp critical
375  {
376  ElementOldCellMatrix[k][l][s].push_back((*(itElemOld.base())));
377  }
378  }
379  }
380  }
381  }
382 
383  // Locate new nodes inside old elements and interpolate nodal variables
384  const int NNodes = static_cast<int>(rModelPartNew.Nodes().size());
385  ModelPart::NodesContainerType::iterator node_begin = rModelPartNew.NodesBegin();
386 
387  array_1d<double,3> GlobalCoordinates;
388  array_1d<double,3> LocalCoordinates;
389 
390  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
391  for(int i = 0; i < NNodes; i++)
392  {
393  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
394 
395  X_me = itNodeNew->X0();
396  Y_me = itNodeNew->Y0();
397  Z_me = itNodeNew->Z0();
398 
399  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
400  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
401  int Section = int((Z_me-AuxVariables.Z_min)/AuxVariables.SectionSize);
402 
403  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
404  else if(Column < 0) Column = 0;
405  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
406  else if(Row < 0) Row = 0;
407  if(Section >= AuxVariables.NSections) Section = AuxVariables.NSections-1;
408  else if(Section < 0) Section = 0;
409 
410  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
411  bool IsInside = false;
412  Element::Pointer pElementOld;
413 
414  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column][Section]).size(); m++)
415  {
416  pElementOld = ElementOldCellMatrix[Row][Column][Section][m];
417  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
418  if(IsInside) break;
419  }
420 
421  if(IsInside==false)
422  {
423  for(int i = 0; i < 12; i++)
424  {
425  for(int j = 1; j<10; j++)
426  {
427  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column][Section]).size(); m++)
428  {
429  pElementOld = ElementOldCellMatrix[Row][Column][Section][m];
430  double tol = j * pow(10.0, (i-12));
431  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,tol);
432  if(IsInside) break;
433  }
434  if(IsInside) break;
435  }
436  if(IsInside) break;
437  }
438  }
439 
440  if(IsInside == false)
441  {
442  std::cout << "None of the old elements contains node " << itNodeNew->Id() << std::endl;
443  if (add_stress)
444  {
445  itNodeNew->FastGetSolutionStepValue(NODAL_CAUCHY_STRESS_TENSOR) = ZeroMatrix(3,3);
446  itNodeNew->FastGetSolutionStepValue(INITIAL_NODAL_CAUCHY_STRESS_TENSOR) = ZeroMatrix(3,3);
447  }
448  continue;
449  }
450 
451  PointsNumber = pElementOld->GetGeometry().PointsNumber();
452  Vector ShapeFunctionsValuesVector(PointsNumber);
453  Vector NodalVariableVector(PointsNumber);
454 
455  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
456 
457  // Interpolation of nodal variables
458  if (add_displacement)
459  {
460  double max_displacement_x = -10e12;
461  double min_displacement_x = 10e12;
462  for(int j = 0; j < PointsNumber; j++)
463  {
464  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_X);
465  if (NodalVariableVector[j] > max_displacement_x) max_displacement_x = NodalVariableVector[j];
466  if (NodalVariableVector[j] < min_displacement_x) min_displacement_x = NodalVariableVector[j];
467  }
468  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
469  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) > max_displacement_x) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = max_displacement_x;
470  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) < min_displacement_x) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = min_displacement_x;
471 
472  double max_displacement_y = -10e12;
473  double min_displacement_y = 10e12;
474  for(int j = 0; j < PointsNumber; j++)
475  {
476  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_Y);
477  if (NodalVariableVector[j] > max_displacement_y) max_displacement_y = NodalVariableVector[j];
478  if (NodalVariableVector[j] < min_displacement_y) min_displacement_y = NodalVariableVector[j];
479  }
480  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
481  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) > max_displacement_y) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = max_displacement_y;
482  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) < min_displacement_y) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = min_displacement_y;
483 
484  double max_displacement_z = -10e12;
485  double min_displacement_z = 10e12;
486  for(int j = 0; j < PointsNumber; j++)
487  {
488  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_Z);
489  if (NodalVariableVector[j] > max_displacement_z) max_displacement_z = NodalVariableVector[j];
490  if (NodalVariableVector[j] < min_displacement_z) min_displacement_z = NodalVariableVector[j];
491  }
492  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Z) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
493  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Z) > max_displacement_z) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Z) = max_displacement_z;
494  if (itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Z) < min_displacement_z) itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Z) = min_displacement_z;
495  }
496 
497  if (add_stress)
498  {
499  Matrix NodalStress = ZeroMatrix(3,3);
500 
501  Matrix MaxNodalStress = ZeroMatrix(3,3);
502  Matrix MinNodalStress = ZeroMatrix(3,3);
503 
504  for(int k = 0; k < 3; k++)
505  {
506  for(int l = 0; l < 3; l++)
507  {
508  MaxNodalStress(k,l) = -10e12;
509  MinNodalStress(k,l) = 10e12;
510  for(int j = 0; j < PointsNumber; j++)
511  {
512  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(NODAL_CAUCHY_STRESS_TENSOR)(k,l);
513  if (NodalVariableVector[j] > MaxNodalStress(k,l)) MaxNodalStress(k,l) = NodalVariableVector[j];
514  if (NodalVariableVector[j] < MinNodalStress(k,l)) MinNodalStress(k,l) = NodalVariableVector[j];
515  }
516  NodalStress(k,l) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
517  if (NodalStress(k,l) > MaxNodalStress(k,l)) NodalStress(k,l) = MaxNodalStress(k,l);
518  if (NodalStress(k,l) < MinNodalStress(k,l)) NodalStress(k,l) = MinNodalStress(k,l);
519  }
520  }
521 
522  itNodeNew->FastGetSolutionStepValue(INITIAL_NODAL_CAUCHY_STRESS_TENSOR) = NodalStress;
523  itNodeNew->FastGetSolutionStepValue(NODAL_CAUCHY_STRESS_TENSOR) = NodalStress;
524  }
525  }
526  }
527 
528 private:
529 
531 
532  void ComputeCellMatrixDimensions(
533  UtilityVariables& rAuxVariables,
534  ModelPart& rModelPart)
535  {
536  // Compute X, Y and Z limits of the current geometry
537  unsigned int NumThreads = ParallelUtilities::GetNumThreads();
538  std::vector<double> X_max_partition(NumThreads);
539  std::vector<double> X_min_partition(NumThreads);
540  std::vector<double> Y_max_partition(NumThreads);
541  std::vector<double> Y_min_partition(NumThreads);
542  std::vector<double> Z_max_partition(NumThreads);
543  std::vector<double> Z_min_partition(NumThreads);
544 
545  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
546  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
547 
548  #pragma omp parallel
549  {
550  int k = OpenMPUtils::ThisThread();
551  X_max_partition[k] = node_begin->X();
552  X_min_partition[k] = X_max_partition[k];
553  Y_max_partition[k] = node_begin->Y();
554  Y_min_partition[k] = Y_max_partition[k];
555  Z_max_partition[k] = node_begin->Z();
556  Z_min_partition[k] = Z_max_partition[k];
557 
558  #pragma omp for
559  for(int i = 0; i < NNodes; i++)
560  {
561  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
562 
563  double X_me = itNode->X();
564  double Y_me = itNode->Y();
565  double Z_me = itNode->Z();
566 
567  if( X_me > X_max_partition[k] ) X_max_partition[k] = X_me;
568  else if( X_me < X_min_partition[k] ) X_min_partition[k] = X_me;
569 
570  if( Y_me > Y_max_partition[k] ) Y_max_partition[k] = Y_me;
571  else if( Y_me < Y_min_partition[k] ) Y_min_partition[k] = Y_me;
572 
573  if( Z_me > Z_max_partition[k] ) Z_max_partition[k] = Z_me;
574  else if( Z_me < Z_min_partition[k] ) Z_min_partition[k] = Z_me;
575  }
576  }
577 
578  rAuxVariables.X_max = X_max_partition[0];
579  rAuxVariables.X_min = X_min_partition[0];
580  rAuxVariables.Y_max = Y_max_partition[0];
581  rAuxVariables.Y_min = Y_min_partition[0];
582  rAuxVariables.Z_max = Z_max_partition[0];
583  rAuxVariables.Z_min = Z_min_partition[0];
584 
585  for(unsigned int i=1; i < NumThreads; i++)
586  {
587  if(X_max_partition[i] > rAuxVariables.X_max) rAuxVariables.X_max = X_max_partition[i];
588  if(X_min_partition[i] < rAuxVariables.X_min) rAuxVariables.X_min = X_min_partition[i];
589  if(Y_max_partition[i] > rAuxVariables.Y_max) rAuxVariables.Y_max = Y_max_partition[i];
590  if(Y_min_partition[i] < rAuxVariables.Y_min) rAuxVariables.Y_min = Y_min_partition[i];
591  if(Z_max_partition[i] > rAuxVariables.Z_max) rAuxVariables.Z_max = Z_max_partition[i];
592  if(Z_min_partition[i] < rAuxVariables.Z_min) rAuxVariables.Z_min = Z_min_partition[i];
593  }
594 
595  // Calculate Average Element Length
596  double AverageElementLength = 0.0;
597 
598  int NElems = static_cast<int>(rModelPart.Elements().size());
599  ModelPart::ElementsContainerType::iterator el_begin = rModelPart.ElementsBegin();
600 
601  #pragma omp parallel for reduction(+:AverageElementLength)
602  for(int i = 0; i < NElems; i++)
603  {
604  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
605 
606  AverageElementLength += itElem->GetGeometry().Length();
607  }
608 
609  AverageElementLength = AverageElementLength/NElems;
610 
611  // Compute FracturePoints CellMatrix dimensions
612  rAuxVariables.NRows = int((rAuxVariables.Y_max-rAuxVariables.Y_min)/AverageElementLength);
613  rAuxVariables.NColumns = int((rAuxVariables.X_max-rAuxVariables.X_min)/AverageElementLength);
614  rAuxVariables.NSections = int((rAuxVariables.Z_max-rAuxVariables.Z_min)/AverageElementLength);
615  if(rAuxVariables.NRows <= 0) rAuxVariables.NRows = 1;
616  if(rAuxVariables.NColumns <= 0) rAuxVariables.NColumns = 1;
617  if(rAuxVariables.NSections <= 0) rAuxVariables.NSections = 1;
618  rAuxVariables.RowSize = (rAuxVariables.Y_max-rAuxVariables.Y_min)/rAuxVariables.NRows;
619  rAuxVariables.ColumnSize = (rAuxVariables.X_max-rAuxVariables.X_min)/rAuxVariables.NColumns;
620  rAuxVariables.SectionSize = (rAuxVariables.Z_max-rAuxVariables.Z_min)/rAuxVariables.NSections;
621  }
622 
623 }; // Class MappingVariables3DUtilities
624 
625 } // namespace Kratos.
626 
627 #endif /* KRATOS_MAPPING_VARIABLES_3D_UTILITIES defined */
Kratos::MappingVariables3DUtilities
Definition: mapping_variables_3D_utilities.hpp:37
Kratos::MappingVariables3DUtilities::NodalThermalVariablesMapping
void NodalThermalVariablesMapping(const UtilityVariables &AuxVariables, ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_temperature, bool add_reference_temperature)
Definition: mapping_variables_3D_utilities.hpp:97
Kratos::MappingVariables3DUtilities::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(MappingVariables3DUtilities)
Kratos::MappingVariables3DUtilities::MappingMechanicalModelParts
void MappingMechanicalModelParts(ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_displacement, bool add_stress)
Definition: mapping_variables_3D_utilities.hpp:72
Kratos::MappingVariables3DUtilities::NodalMechanicalVariablesMapping
void NodalMechanicalVariablesMapping(const UtilityVariables &AuxVariables, ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_displacement, bool add_stress)
Definition: mapping_variables_3D_utilities.hpp:292
Kratos::MappingVariables3DUtilities::MappingThermalModelParts
void MappingThermalModelParts(ModelPart &rModelPartOld, ModelPart &rModelPartNew, bool add_temperature, bool add_reference_temperature)
Definition: mapping_variables_3D_utilities.hpp:64
Kratos::MappingVariables3DUtilities::InitializeMapping
void InitializeMapping(UtilityVariables &rAuxVariables, ModelPart &rModelPartNew)
Member Variables.
Definition: mapping_variables_3D_utilities.hpp:88
Kratos::MappingVariables3DUtilities::MappingVariables3DUtilities
MappingVariables3DUtilities()
Constructor.
Definition: mapping_variables_3D_utilities.hpp:55
Kratos::MappingVariables3DUtilities::~MappingVariables3DUtilities
virtual ~MappingVariables3DUtilities()
Destructor.
Definition: mapping_variables_3D_utilities.hpp:60
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::MappingVariables3DUtilities::UtilityVariables
Basic Structs for the utility -----------------------------------------------------------------------...
Definition: mapping_variables_3D_utilities.hpp:44
Kratos::MappingVariables3DUtilities::UtilityVariables::ColumnSize
double ColumnSize
Definition: mapping_variables_3D_utilities.hpp:47
Kratos::MappingVariables3DUtilities::UtilityVariables::Z_min
double Z_min
Definition: mapping_variables_3D_utilities.hpp:45
Kratos::MappingVariables3DUtilities::UtilityVariables::Y_min
double Y_min
Definition: mapping_variables_3D_utilities.hpp:45
Kratos::MappingVariables3DUtilities::UtilityVariables::X_max
double X_max
Definition: mapping_variables_3D_utilities.hpp:45
Kratos::MappingVariables3DUtilities::UtilityVariables::SectionSize
double SectionSize
Definition: mapping_variables_3D_utilities.hpp:47
Kratos::MappingVariables3DUtilities::UtilityVariables::NColumns
int NColumns
Definition: mapping_variables_3D_utilities.hpp:46
Kratos::MappingVariables3DUtilities::UtilityVariables::NSections
int NSections
Definition: mapping_variables_3D_utilities.hpp:46
Kratos::MappingVariables3DUtilities::UtilityVariables::NRows
int NRows
Definition: mapping_variables_3D_utilities.hpp:46
Kratos::MappingVariables3DUtilities::UtilityVariables::RowSize
double RowSize
Definition: mapping_variables_3D_utilities.hpp:47
Kratos::MappingVariables3DUtilities::UtilityVariables::Z_max
double Z_max
Definition: mapping_variables_3D_utilities.hpp:45
Kratos::MappingVariables3DUtilities::UtilityVariables::Y_max
double Y_max
Definition: mapping_variables_3D_utilities.hpp:45
Kratos::MappingVariables3DUtilities::UtilityVariables::X_min
double X_min
Definition: mapping_variables_3D_utilities.hpp:45
Generated on Sat Feb 3 2024 01:59:56 for KratosMultiphysics by doxygen 1.9.1