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.
initial_stress_2D_utilities.hpp
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: Ignasi de Pouplana
11 //
12 
13 #if !defined(KRATOS_INITIAL_STRESS_2D_UTILITIES )
14 #define KRATOS_INITIAL_STRESS_2D_UTILITIES
15 
16 // System includes
17 #include <fstream>
18 #include <iostream>
19 #include <cmath>
20 
21 // Project includes
22 #include "geometries/geometry.h"
23 #include "includes/define.h"
24 #include "includes/model_part.h"
26 #include "utilities/openmp_utils.h"
28 #include "utilities/math_utils.h"
29 
30 // Application includes
32 
33 namespace Kratos
34 {
35 
37 {
38 
39 protected:
40 
42 
44  {
45  double X_max, X_min, Y_max, Y_min;
48  };
49 
50 public:
51 
53 
56 
58 
61 
63 
64  void TransferInitialStresses (ModelPart& rInitialModelPart, ModelPart& rCurrentModelPart)
65  {
66  // Define necessary variables
67  UtilityVariables AuxVariables;
68 
69  this->ComputeCellMatrixDimensions(AuxVariables,rCurrentModelPart);
70 
71  this->NodalVariablesMapping(AuxVariables,rInitialModelPart,rCurrentModelPart);
72  }
73 
75 
76  void SaveInitialStresses (Parameters& rParameters, ModelPart& rCurrentModelPart)
77  {
78  std::fstream initial_stresses_mdpa;
79  std::string initial_stress_mdpa_path = rParameters["initial_input_filename"].GetString() + ".mdpa";
80 
81  initial_stresses_mdpa.open(initial_stress_mdpa_path.c_str(), std::fstream::out | std::fstream::trunc);
82  initial_stresses_mdpa.precision(15);
83 
84  initial_stresses_mdpa << "Begin Properties 0" << std::endl;
85  initial_stresses_mdpa << "End Properties" << std::endl << std::endl;
86 
87  const int NNodes = static_cast<int>(rCurrentModelPart.Nodes().size());
88  ModelPart::NodesContainerType::iterator node_begin = rCurrentModelPart.NodesBegin();
89  initial_stresses_mdpa << "Begin Nodes" << std::endl;
90  // #pragma omp parallel for
91  for(int i = 0; i < NNodes; i++) {
92  ModelPart::NodesContainerType::iterator it_node = node_begin + i;
93  initial_stresses_mdpa << it_node->Id() << " " <<
94  it_node->X0() << " " <<
95  it_node->Y0() << " " <<
96  it_node->Z0() << std::endl;
97  }
98  initial_stresses_mdpa << "End Nodes" << std::endl << std::endl;
99 
100  this->WriteLinearElements(initial_stresses_mdpa,rCurrentModelPart);
101 
102  this->CalculateNodalStresses(rCurrentModelPart);
103  Matrix InitialStressTensor(3,3);
104  initial_stresses_mdpa << "Begin NodalData INITIAL_STRESS_TENSOR" << std::endl;
105  // #pragma omp parallel for
106  for(int i = 0; i < NNodes; i++) {
107  ModelPart::NodesContainerType::iterator it_node = node_begin + i;
108  noalias(InitialStressTensor) = it_node->FastGetSolutionStepValue(NODAL_EFFECTIVE_STRESS_TENSOR);
109  initial_stresses_mdpa << it_node->Id() << " 0 [2,2] ((" <<
110  InitialStressTensor(0,0) << "," << InitialStressTensor(0,1) << "),(" <<
111  InitialStressTensor(1,0) << "," << InitialStressTensor(1,1) << "))" << std::endl;
112  }
113  initial_stresses_mdpa << "End NodalData" << std::endl;
114 
115  initial_stresses_mdpa.close();
116  }
117 
118 protected:
119 
121 
123 
125  const UtilityVariables& AuxVariables,
126  ModelPart& rInitialModelPart,
127  ModelPart& rCurrentModelPart)
128  {
129  // Define ElementOld Cell matrix
130  std::vector< std::vector< std::vector<Element::Pointer> > > ElementOldCellMatrix;
131  ElementOldCellMatrix.resize(AuxVariables.NRows);
132  for(int i = 0; i < AuxVariables.NRows; i++) ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
133 
134  // Locate Old Elments in cells
135  double X_me;
136  double Y_me;
137  int PointsNumber;
138 
139  int NElems = static_cast<int>(rInitialModelPart.Elements().size());
140  ModelPart::ElementsContainerType::iterator el_begin = rInitialModelPart.ElementsBegin();
141 
142  #pragma omp parallel for private(X_me,Y_me,PointsNumber)
143  for(int i = 0; i < NElems; i++)
144  {
145  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
146 
147  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
148  double X_right = X_left;
149  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
150  double Y_bot = Y_top;
151  PointsNumber = itElemOld->GetGeometry().PointsNumber();
152 
153  for(int j = 1; j < PointsNumber; j++)
154  {
155  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
156  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
157 
158  if(X_me > X_right) X_right = X_me;
159  else if(X_me < X_left) X_left = X_me;
160  if(Y_me > Y_top) Y_top = Y_me;
161  else if(Y_me < Y_bot) Y_bot = Y_me;
162  }
163 
164  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
165  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
166  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
167  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
168 
169  if(Column_left < 0) Column_left = 0;
170  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
171  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
172  else if(Column_right < 0) Column_right = 0;
173 
174  if(Row_top < 0) Row_top = 0;
175  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
176  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
177  else if(Row_bot < 0) Row_bot = 0;
178 
179  for(int k = Row_top; k <= Row_bot; k++)
180  {
181  for(int l = Column_left; l <= Column_right; l++)
182  {
183  #pragma omp critical
184  {
185  ElementOldCellMatrix[k][l].push_back((*(itElemOld.base())));
186  }
187  }
188  }
189  }
190 
191  // Locate current nodes inside initial elements and interpolate nodal variables
192  const int NNodes = static_cast<int>(rCurrentModelPart.Nodes().size());
193  ModelPart::NodesContainerType::iterator node_begin = rCurrentModelPart.NodesBegin();
194 
195  array_1d<double,3> GlobalCoordinates;
196  array_1d<double,3> LocalCoordinates;
197 
198  #pragma omp parallel for private(X_me,Y_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
199  for(int i = 0; i < NNodes; i++)
200  {
201  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
202 
203  X_me = itNodeNew->X0();
204  Y_me = itNodeNew->Y0();
205 
206  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
207  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
208 
209  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
210  else if(Column < 0) Column = 0;
211  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
212  else if(Row < 0) Row = 0;
213 
214  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
215  bool IsInside = false;
216  Element::Pointer pElementOld;
217 
218  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
219  {
220  pElementOld = ElementOldCellMatrix[Row][Column][m];
221  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
222  if(IsInside) break;
223  }
224  if(IsInside==false)
225  {
226  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
227  {
228  pElementOld = ElementOldCellMatrix[Row][Column][m];
229  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,1.0e-5);
230  if(IsInside) break;
231  }
232  }
233  if(IsInside == false)
234  std::cout << "ERROR!!, NONE OF THE OLD ELEMENTS CONTAINS NODE: " << itNodeNew->Id() << std::endl;
235 
236  PointsNumber = pElementOld->GetGeometry().PointsNumber();
237 
238  Vector ShapeFunctionsValuesVector(PointsNumber);
239  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
240 
241  // Interpolation of nodal variables
242  std::vector<Vector> ComponentsNodalVariableVector(3); // voigt_size
243  for(int j = 0; j < 3; j++) { // voigt_size
244  ComponentsNodalVariableVector[j].resize(PointsNumber);
245  }
246  Vector nodal_initial_stress_vector(3); // voigt_size
247  Matrix nodal_initial_stress_tensor(2,2); // dimension
248  for(int j = 0; j < PointsNumber; j++) {
249  noalias(nodal_initial_stress_tensor) = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(INITIAL_STRESS_TENSOR);
250  noalias(nodal_initial_stress_vector) = MathUtils<double>::StressTensorToVector(nodal_initial_stress_tensor);
251  for (int k= 0; k < 3; k++) { // voigt_size
252  ComponentsNodalVariableVector[k][j] = nodal_initial_stress_vector[k];
253  }
254  }
255  for (int k= 0; k < 3; k++) { // voigt_size
256  nodal_initial_stress_vector[k] = inner_prod(ShapeFunctionsValuesVector,ComponentsNodalVariableVector[k]);
257  }
258  noalias(nodal_initial_stress_tensor) = MathUtils<double>::StressVectorToTensor(nodal_initial_stress_vector);
259  Matrix& rNodalStress = itNodeNew->FastGetSolutionStepValue(INITIAL_STRESS_TENSOR);
260  if(rNodalStress.size1() != 2) // Dimension
261  rNodalStress.resize(2,2,false);
262  noalias(rNodalStress) = nodal_initial_stress_tensor;
263  }
264  }
265 
267 
268 private:
269 
271 
272  void ComputeCellMatrixDimensions(
273  UtilityVariables& rAuxVariables,
274  ModelPart& rModelPart)
275  {
276  // Compute X and Y limits of the current geometry
277  unsigned int NumThreads = ParallelUtilities::GetNumThreads();
278  std::vector<double> X_max_partition(NumThreads);
279  std::vector<double> X_min_partition(NumThreads);
280  std::vector<double> Y_max_partition(NumThreads);
281  std::vector<double> Y_min_partition(NumThreads);
282 
283  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
284  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
285 
286  #pragma omp parallel
287  {
288  int k = OpenMPUtils::ThisThread();
289 
290  X_max_partition[k] = node_begin->X();
291  X_min_partition[k] = X_max_partition[k];
292  Y_max_partition[k] = node_begin->Y();
293  Y_min_partition[k] = Y_max_partition[k];
294 
295  double X_me, Y_me;
296 
297  #pragma omp for
298  for(int i = 0; i < NNodes; i++)
299  {
300  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
301 
302  X_me = itNode->X();
303  Y_me = itNode->Y();
304 
305  if( X_me > X_max_partition[k] ) X_max_partition[k] = X_me;
306  else if( X_me < X_min_partition[k] ) X_min_partition[k] = X_me;
307 
308  if( Y_me > Y_max_partition[k] ) Y_max_partition[k] = Y_me;
309  else if( Y_me < Y_min_partition[k] ) Y_min_partition[k] = Y_me;
310  }
311  }
312 
313  rAuxVariables.X_max = X_max_partition[0];
314  rAuxVariables.X_min = X_min_partition[0];
315  rAuxVariables.Y_max = Y_max_partition[0];
316  rAuxVariables.Y_min = Y_min_partition[0];
317 
318  for(unsigned int i=1; i < NumThreads; i++)
319  {
320  if(X_max_partition[i] > rAuxVariables.X_max) rAuxVariables.X_max = X_max_partition[i];
321  if(X_min_partition[i] < rAuxVariables.X_min) rAuxVariables.X_min = X_min_partition[i];
322  if(Y_max_partition[i] > rAuxVariables.Y_max) rAuxVariables.Y_max = Y_max_partition[i];
323  if(Y_min_partition[i] < rAuxVariables.Y_min) rAuxVariables.Y_min = Y_min_partition[i];
324  }
325 
326  // Calculate Average Element Length
327  double AverageElementLength = 0.0;
328 
329  int NElems = static_cast<int>(rModelPart.Elements().size());
330  ModelPart::ElementsContainerType::iterator el_begin = rModelPart.ElementsBegin();
331 
332  #pragma omp parallel for reduction(+:AverageElementLength)
333  for(int i = 0; i < NElems; i++)
334  {
335  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
336 
337  AverageElementLength += itElem->GetGeometry().Length();
338  }
339 
340  AverageElementLength = AverageElementLength/NElems;
341 
342  // Compute FracturePoints CellMatrix dimensions
343  rAuxVariables.NRows = int((rAuxVariables.Y_max-rAuxVariables.Y_min)/(AverageElementLength));
344  rAuxVariables.NColumns = int((rAuxVariables.X_max-rAuxVariables.X_min)/(AverageElementLength));
345  if(rAuxVariables.NRows <= 0) rAuxVariables.NRows = 1;
346  if(rAuxVariables.NColumns <= 0) rAuxVariables.NColumns = 1;
347  rAuxVariables.RowSize = (rAuxVariables.Y_max-rAuxVariables.Y_min)/rAuxVariables.NRows;
348  rAuxVariables.ColumnSize = (rAuxVariables.X_max-rAuxVariables.X_min)/rAuxVariables.NColumns;
349  }
350 
352 
353  void CalculateNodalStresses(ModelPart& rCurrentModelPart) {
354 
355  // Clear nodal Stresses
356  const int NNodes = static_cast<int>(rCurrentModelPart.Nodes().size());
357  ModelPart::NodesContainerType::iterator node_begin = rCurrentModelPart.NodesBegin();
358 
359  #pragma omp parallel for
360  for(int i = 0; i < NNodes; i++)
361  {
362  ModelPart::NodesContainerType::iterator it_node = node_begin + i;
363  it_node->FastGetSolutionStepValue(NODAL_AREA) = 0.0;
364  Matrix& rNodalStress = it_node->FastGetSolutionStepValue(NODAL_EFFECTIVE_STRESS_TENSOR);
365  if(rNodalStress.size1() != 3)
366  rNodalStress.resize(3,3,false);
367  noalias(rNodalStress) = ZeroMatrix(3,3);
368  }
369 
370  // Calculate and Extrapolate Stresses
371  const ProcessInfo& r_current_process_info = rCurrentModelPart.GetProcessInfo();
372  const auto it_elem_begin = rCurrentModelPart.ElementsBegin();
373  #pragma omp parallel for
374  for(int i=0; i<static_cast<int>(rCurrentModelPart.Elements().size()); ++i) {
375  auto it_elem = it_elem_begin + i;
376  it_elem->FinalizeSolutionStep(r_current_process_info);
377  }
378 
379  // Compute smoothed Stresses
380  #pragma omp parallel for
381  for(int n = 0; n < NNodes; n++)
382  {
383  ModelPart::NodesContainerType::iterator it_node = node_begin + n;
384  const double& NodalArea = it_node->FastGetSolutionStepValue(NODAL_AREA);
385  if (NodalArea>1.0e-20)
386  {
387  const double InvNodalArea = 1.0/NodalArea;
388  Matrix& rNodalStress = it_node->FastGetSolutionStepValue(NODAL_EFFECTIVE_STRESS_TENSOR);
389  for(unsigned int i = 0; i<3; i++) // Dimension
390  {
391  for(unsigned int j = 0; j<3; j++)
392  {
393  rNodalStress(i,j) *= InvNodalArea;
394  }
395  }
396  }
397  }
398  }
399 
401 
402  void WriteLinearElements(std::fstream& rInitialStressesMdpa, ModelPart& rCurrentModelPart) {
403 
404  int NElems = static_cast<int>(rCurrentModelPart.Elements().size());
405  ModelPart::ElementsContainerType::iterator el_begin = rCurrentModelPart.ElementsBegin();
406  // NOTE: we are assuming that only one type of element is used (not quadratic elements)
407  int PointsNumber = el_begin->GetGeometry().PointsNumber();
408  if(PointsNumber == 3){
409  rInitialStressesMdpa << "Begin Elements Element2D3N" << std::endl;
410  // #pragma omp parallel for
411  for(int i = 0; i < NElems; i++)
412  {
413  ModelPart::ElementsContainerType::iterator it_elem = el_begin + i;
414  rInitialStressesMdpa << it_elem->Id() << " 0 " <<
415  it_elem->GetGeometry().GetPoint(0).Id() << " " <<
416  it_elem->GetGeometry().GetPoint(1).Id() << " " <<
417  it_elem->GetGeometry().GetPoint(2).Id() << std::endl;
418  }
419  } else {
420  rInitialStressesMdpa << "Begin Elements Element2D4N" << std::endl;
421  // #pragma omp parallel for
422  for(int i = 0; i < NElems; i++)
423  {
424  ModelPart::ElementsContainerType::iterator it_elem = el_begin + i;
425  rInitialStressesMdpa << it_elem->Id() << " 0 " <<
426  it_elem->GetGeometry().GetPoint(0).Id() << " " <<
427  it_elem->GetGeometry().GetPoint(1).Id() << " " <<
428  it_elem->GetGeometry().GetPoint(2).Id() << " " <<
429  it_elem->GetGeometry().GetPoint(3).Id() << std::endl;
430  }
431  }
432  rInitialStressesMdpa << "End Elements" << std::endl << std::endl;
433  }
434 
435 }; // Class InitialStress2DUtilities
436 
437 } // namespace Kratos.
438 
439 #endif /* KRATOS_INITIAL_STRESS_2D_UTILITIES defined */
Definition: initial_stress_2D_utilities.hpp:37
virtual ~InitialStress2DUtilities()
Destructor.
Definition: initial_stress_2D_utilities.hpp:60
void NodalVariablesMapping(const UtilityVariables &AuxVariables, ModelPart &rInitialModelPart, ModelPart &rCurrentModelPart)
Member Variables.
Definition: initial_stress_2D_utilities.hpp:124
InitialStress2DUtilities()
Constructor.
Definition: initial_stress_2D_utilities.hpp:55
void TransferInitialStresses(ModelPart &rInitialModelPart, ModelPart &rCurrentModelPart)
Definition: initial_stress_2D_utilities.hpp:64
KRATOS_CLASS_POINTER_DEFINITION(InitialStress2DUtilities)
void SaveInitialStresses(Parameters &rParameters, ModelPart &rCurrentModelPart)
Definition: initial_stress_2D_utilities.hpp:76
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
static TMatrixType StressVectorToTensor(const TVector &rStressVector)
Transforms a stess vector into a matrix. Stresses are assumed to be stored in the following way: for...
Definition: math_utils.h:1174
static TVector StressTensorToVector(const TMatrixType &rStressTensor, SizeType rSize=0)
Transforms a given symmetric Stress Tensor to Voigt Notation:
Definition: math_utils.h:1399
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
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
This class provides to Kratos a data structure for I/O based on the standard of JSON.
Definition: kratos_parameters.h:59
std::string GetString() const
This method returns the string contained in the current Parameter.
Definition: kratos_parameters.cpp:684
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
ProcessInfo
Definition: edgebased_PureConvection.py:116
out
Definition: isotropic_damage_automatic_differentiation.py:200
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
int m
Definition: run_marine_rain_substepping.py:8
integer i
Definition: TensorModule.f:17
integer l
Definition: TensorModule.f:17
Basic Structs for the utility -----------------------------------------------------------------------...
Definition: initial_stress_2D_utilities.hpp:44
int NRows
Definition: initial_stress_2D_utilities.hpp:46
int NColumns
Definition: initial_stress_2D_utilities.hpp:46
double Y_max
Definition: initial_stress_2D_utilities.hpp:45
double ColumnSize
Definition: initial_stress_2D_utilities.hpp:47
double X_min
Definition: initial_stress_2D_utilities.hpp:45
double X_max
Definition: initial_stress_2D_utilities.hpp:45
double Y_min
Definition: initial_stress_2D_utilities.hpp:45
double RowSize
Definition: initial_stress_2D_utilities.hpp:47