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.
fracture_propagation_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: Ignasi de Pouplana
12 //
13 
14 
15 #if !defined(KRATOS_PROPAGATE_FRACTURES_3D_UTILITIES )
16 #define KRATOS_PROPAGATE_FRACTURES_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 #include "utilities/parallel_utilities.h"
30 
31 // Application includes
32 #include "poromechanics_application_variables.h"
33 
34 namespace Kratos
35 {
36 
37 class FracturePropagation3DUtilities
38 {
39 
40 protected:
41 
43 
44  struct UtilityVariables
45  {
46  double X_max, X_min, Y_max, Y_min, Z_max, Z_min;
47  int NRows, NColumns, NSections;
48  double RowSize, ColumnSize, SectionSize;
49  };
50 
52 
53  struct FracturePoint
54  {
55  array_1d<double,3> Coordinates;
56  double Damage, Weight, TipDistance;
57  Element::Pointer pElement;
58  };
59 
61 
62  struct Propagation
63  {
64  int MotherFractureId;
65  array_1d<double,3> TopInitCoordinates;
66  array_1d<double,3> BotInitCoordinates;
67  array_1d<double,3> LeftInitCoordinates;
68  array_1d<double,3> RightInitCoordinates;
69  array_1d<double,3> TopEndCoordinates;
70  array_1d<double,3> BotEndCoordinates;
71  array_1d<double,3> LeftEndCoordinates;
72  array_1d<double,3> RightEndCoordinates;
73  array_1d<double,3> TipCoordinates;
74  };
75 
77 
78  struct Bifurcation
79  {
80  int MotherFractureId;
81  array_1d<double,3> TopInitCoordinates;
82  array_1d<double,3> BotInitCoordinates;
83  array_1d<double,3> LeftInitCoordinates;
84  array_1d<double,3> RightInitCoordinates;
85  array_1d<double,3> TopTopEndCoordinates;
86  array_1d<double,3> TopBotEndCoordinates;
87  array_1d<double,3> TopLeftEndCoordinates;
88  array_1d<double,3> TopRightEndCoordinates;
89  array_1d<double,3> TopTipCoordinates;
90  array_1d<double,3> BotTopEndCoordinates;
91  array_1d<double,3> BotBotEndCoordinates;
92  array_1d<double,3> BotLeftEndCoordinates;
93  array_1d<double,3> BotRightEndCoordinates;
94  array_1d<double,3> BotTipCoordinates;
95  };
96 
98 
99  struct PropagationGlobalVariables
100  {
101  std::vector< std::vector< std::vector< std::vector<FracturePoint> > > > FracturePointsCellMatrix;
102  ProcessInfo::Pointer pProcessInfo;
103  bool PropagateFractures;
104  std::vector<Propagation> PropagationVector;
105  std::vector<Bifurcation> BifurcationVector;
106  };
107 
109 
110  struct PropagationLocalVariables
111  {
112  BoundedMatrix<double,3,3> RotationMatrix;
113  array_1d<double,3> TipCoordinates;
114  array_1d<double,3> TipLocalCoordinates;
115  std::vector<FracturePoint*> TopFrontFracturePoints;
116  std::vector<FracturePoint*> BotFrontFracturePoints;
117  };
118 
120 
121  struct GaussPointOld
122  {
123  array_1d<double,3> Coordinates;
124  double StateVariable, Weight;
125  };
126 
127 public:
128 
129  KRATOS_CLASS_POINTER_DEFINITION( FracturePropagation3DUtilities );
130 
132  FracturePropagation3DUtilities() {}
133 
135 
137  virtual ~FracturePropagation3DUtilities() {}
138 
140 
141  bool CheckFracturePropagation (Parameters& rParameters, ModelPart& rModelPart, const bool& move_mesh_flag)
142  {
143  // Define necessary variables
144  UtilityVariables AuxVariables;
145  PropagationGlobalVariables PropagationData;
146 
147  this->InitializeCheckFracture(PropagationData, AuxVariables, rParameters, rModelPart, move_mesh_flag);
148 
149  //Loop for all pre-existing fractures
150  for(unsigned int i = 0; i < rParameters["fractures_list"].size(); i++)
151  {
152  this->CheckFracture(i, PropagationData, AuxVariables, rParameters);
153  }
154 
155  this->FinalizeCheckFracture(PropagationData, rParameters, rModelPart, move_mesh_flag);
156 
157  return PropagationData.PropagateFractures;
158  }
159 
161 
162  void MappingModelParts (Parameters& rParameters, ModelPart& rModelPartOld, ModelPart& rModelPartNew, const bool& move_mesh_flag)
163  {
164  // Define necessary variables
165  UtilityVariables AuxVariables;
166 
167  this->InitializeMapping(AuxVariables,rModelPartNew, move_mesh_flag);
168 
169  this->NodalVariablesMapping(AuxVariables,rParameters,rModelPartOld,rModelPartNew);
170 
171  this->GaussPointStateVariableMapping(AuxVariables,rParameters,rModelPartOld,rModelPartNew);
172 
173  if(move_mesh_flag==true)
174  this->UpdateMeshPosition(rModelPartNew);
175  }
176 
177 protected:
178 
180 
182 
184 
185  void InitializeCheckFracture(
186  PropagationGlobalVariables& rPropagationData,
187  UtilityVariables& rAuxVariables,
188  Parameters& rParameters,
189  ModelPart& rModelPart,
190  const bool& move_mesh_flag)
191  {
192  // Set PropagateFractures bool to false
193  rPropagationData.PropagateFractures = false;
194 
195  // Move mesh to the original position to work in the reference state
196  if(move_mesh_flag==true)
197  this->ResetMeshPosition(rModelPart);
198 
199  // Locate fracture points inside CellMatrix
200  this->SetFracturePoints(rPropagationData, rAuxVariables, rParameters, rModelPart);
201  }
202 
204 
205  void CheckFracture(
206  const unsigned int& itFracture,
207  PropagationGlobalVariables& rPropagationData,
208  const UtilityVariables& AuxVariables,
209  Parameters& rParameters)
210  {
211  PropagationLocalVariables AuxPropagationVariables;
212 
213  // Original tip Coordinates
214  array_1d<double,3> LeftPointCoordinates;
215  array_1d<double,3> RightPointCoordinates;
216  for(unsigned int i = 0; i < 3; i++)
217  {
218  AuxPropagationVariables.TipCoordinates[i] = rParameters["fractures_list"][itFracture]["tip_point"]["coordinates"][i].GetDouble();
219  LeftPointCoordinates[i] = rParameters["fractures_list"][itFracture]["left_point"]["coordinates"][i].GetDouble();
220  RightPointCoordinates[i] = rParameters["fractures_list"][itFracture]["right_point"]["coordinates"][i].GetDouble();
221  }
222 
223  // Tip Rotation Matrix
224  this->CalculateTipRotationMatrix(AuxPropagationVariables.RotationMatrix,AuxPropagationVariables.TipCoordinates,LeftPointCoordinates,RightPointCoordinates);
225 
226  // Tip search area
227  const double PropagationLength = rParameters["fracture_data"]["propagation_length"].GetDouble();
228 
229  double X_left = AuxPropagationVariables.TipCoordinates[0] - PropagationLength;
230  double X_right = AuxPropagationVariables.TipCoordinates[0] + PropagationLength;
231  double Y_top = AuxPropagationVariables.TipCoordinates[1] + PropagationLength;
232  double Y_bot = AuxPropagationVariables.TipCoordinates[1] - PropagationLength;
233  double Z_back = AuxPropagationVariables.TipCoordinates[2] - PropagationLength;
234  double Z_front = AuxPropagationVariables.TipCoordinates[2] + PropagationLength;
235 
236  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
237  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
238  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
239  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
240  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
241  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
242 
243  if(Column_left < 0) Column_left = 0;
244  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
245  if(Row_top < 0) Row_top = 0;
246  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
247  if(Section_back < 0) Section_back = 0;
248  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
249 
250  // Search FracturePoints neighbours around the tip
251  std::vector<FracturePoint*> TipNeighbours;
252  noalias(AuxPropagationVariables.TipLocalCoordinates) = prod(AuxPropagationVariables.RotationMatrix,AuxPropagationVariables.TipCoordinates);
253  array_1d<double,3> OtherLocalCoordinates;
254  for(int s = Section_back; s <= Section_front; s++)
255  {
256  for(int i = Row_top; i <= Row_bot; i++)
257  {
258  for(int j = Column_left; j<= Column_right; j++)
259  {
260  for(unsigned int k = 0; k < rPropagationData.FracturePointsCellMatrix[i][j][s].size(); k++)
261  {
262  FracturePoint& rOtherPoint = rPropagationData.FracturePointsCellMatrix[i][j][s][k];
263 
264  rOtherPoint.TipDistance = sqrt((rOtherPoint.Coordinates[0]-AuxPropagationVariables.TipCoordinates[0])*(rOtherPoint.Coordinates[0]-AuxPropagationVariables.TipCoordinates[0]) +
265  (rOtherPoint.Coordinates[1]-AuxPropagationVariables.TipCoordinates[1])*(rOtherPoint.Coordinates[1]-AuxPropagationVariables.TipCoordinates[1]) +
266  (rOtherPoint.Coordinates[2]-AuxPropagationVariables.TipCoordinates[2])*(rOtherPoint.Coordinates[2]-AuxPropagationVariables.TipCoordinates[2]));
267 
268  if(rOtherPoint.TipDistance <= PropagationLength)
269  {
270  TipNeighbours.push_back(&rOtherPoint);
271 
272  noalias(OtherLocalCoordinates) = prod(AuxPropagationVariables.RotationMatrix,rOtherPoint.Coordinates);
273 
274  // FrontFracturePoints
275  if(OtherLocalCoordinates[0] >= AuxPropagationVariables.TipLocalCoordinates[0])
276  {
277  // TopFrontFracturePoints
278  if(OtherLocalCoordinates[2] >= AuxPropagationVariables.TipLocalCoordinates[2])
279  {
280  AuxPropagationVariables.TopFrontFracturePoints.push_back(&rOtherPoint);
281  }
282  // BotFrontFracturePoints
283  else
284  {
285  AuxPropagationVariables.BotFrontFracturePoints.push_back(&rOtherPoint);
286  }
287  }
288  }
289  }
290  }
291  }
292  }
293 
294  // Calculate Non-local damage around the tip
295  double NonlocalDamage = 0.0;
296  double WeightingFunctionDenominator = 0.0;
297  int NNeighbours = TipNeighbours.size();
298 
299  #pragma omp parallel for reduction(+:NonlocalDamage,WeightingFunctionDenominator)
300  for(int i = 0; i < NNeighbours; i++)
301  {
302  const FracturePoint& MyPoint = *(TipNeighbours[i]);
303  NonlocalDamage += (MyPoint.Weight)*
304  exp(-4.0*(MyPoint.TipDistance)*(MyPoint.TipDistance)/(PropagationLength*PropagationLength))*
305  (MyPoint.Damage);
306  WeightingFunctionDenominator += (MyPoint.Weight)*exp(-4.0*(MyPoint.TipDistance)*(MyPoint.TipDistance)/(PropagationLength*PropagationLength));
307  }
308 
309  if(WeightingFunctionDenominator > 1.0e-20)
310  NonlocalDamage = NonlocalDamage/WeightingFunctionDenominator;
311  else
312  NonlocalDamage = 0.0;
313 
314  // Check fracture propagation
315  if(NonlocalDamage >= rParameters["fracture_data"]["propagation_damage"].GetDouble())
316  this->PropagateFracture(itFracture,rPropagationData,AuxPropagationVariables,rParameters);
317  }
318 
320 
321  void FinalizeCheckFracture(
322  const PropagationGlobalVariables& PropagationData,
323  Parameters& rParameters,
324  ModelPart& rModelPart,
325  const bool& move_mesh_flag)
326  {
327  if(PropagationData.PropagateFractures==false)
328  {
329  // Move mesh to the current position
330  if(move_mesh_flag==true)
331  this->UpdateMeshPosition(rModelPart);
332  }
333  else
334  {
335  // Write PropagationData.tcl file
336  this->WritePropagationData(PropagationData,rParameters);
337  }
338  }
339 
341 
342  void WritePropagationData(
343  const PropagationGlobalVariables& PropagationData,
344  Parameters& rParameters)
345  {
346 
347  std::fstream PropDataFile;
348  PropDataFile.open ("PropagationData.tcl", std::fstream::out | std::fstream::trunc);
349  PropDataFile.precision(12);
350 
351  PropDataFile << "set PropagationData [list]" << std::endl;
352 
353  // GiDPath
354  PropDataFile << "lappend PropagationData \"" << rParameters["fracture_data"]["gid_path"].GetString() << "\"" << std::endl;
355 
356  int Id;
357 
358  // BodyVolumesDict
359  for(unsigned int i = 0; i < rParameters["body_volumes_list"].size(); i++)
360  {
361  Id = rParameters["body_volumes_list"][i]["id"].GetInt();
362 
363  PropDataFile << "set Groups [list]" << std::endl;
364  for(unsigned int j = 0; j < rParameters["body_volumes_list"][i]["groups"].size(); j++)
365  {
366  PropDataFile << "lappend Groups \"" << rParameters["body_volumes_list"][i]["groups"][j].GetString() << "\"" << std::endl;
367  }
368  PropDataFile << "dict set BodyVolumesDict " << Id << " Groups $Groups" << std::endl;
369  PropDataFile << "set Surfaces [list]" << std::endl;
370  for(unsigned int j = 0; j < rParameters["body_volumes_list"][i]["surfaces"].size(); j++)
371  {
372  PropDataFile << "lappend Surfaces " << rParameters["body_volumes_list"][i]["surfaces"][j].GetInt() << std::endl;
373  }
374  PropDataFile << "dict set BodyVolumesDict " << Id << " Surfaces $Surfaces" << std::endl;
375  PropDataFile << "dict set BodyVolumesDict " << Id << " MeshSize " << rParameters["body_volumes_list"][i]["mesh_size"].GetDouble() << std::endl;
376  }
377  PropDataFile << "lappend PropagationData $BodyVolumesDict" << std::endl;
378 
379  // FracturesDict
380  for(unsigned int i = 0; i < rParameters["fractures_list"].size(); i++)
381  {
382  Id = rParameters["fractures_list"][i]["id"].GetInt();
383 
384  // TipPoint
385  PropDataFile << "dict set FracturesDict " << Id << " TipPoint Id "
386  << rParameters["fractures_list"][i]["tip_point"]["id"].GetInt() << std::endl;
387  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["tip_point"]["coordinates"][0].GetDouble()
388  << " " << rParameters["fractures_list"][i]["tip_point"]["coordinates"][1].GetDouble() << " "
389  << rParameters["fractures_list"][i]["tip_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
390  PropDataFile << "dict set FracturesDict " << Id << " TipPoint Coordinates $Coordinates" << std::endl;
391  // TopPoint
392  PropDataFile << "dict set FracturesDict " << Id << " TopPoint Id "
393  << rParameters["fractures_list"][i]["top_point"]["id"].GetInt() << std::endl;
394  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["top_point"]["coordinates"][0].GetDouble()
395  << " " << rParameters["fractures_list"][i]["top_point"]["coordinates"][1].GetDouble() << " "
396  << rParameters["fractures_list"][i]["top_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
397  PropDataFile << "dict set FracturesDict " << Id << " TopPoint Coordinates $Coordinates" << std::endl;
398  // BotPoint
399  PropDataFile << "dict set FracturesDict " << Id << " BotPoint Id "
400  << rParameters["fractures_list"][i]["bot_point"]["id"].GetInt() << std::endl;
401  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["bot_point"]["coordinates"][0].GetDouble()
402  << " " << rParameters["fractures_list"][i]["bot_point"]["coordinates"][1].GetDouble() << " "
403  << rParameters["fractures_list"][i]["bot_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
404  PropDataFile << "dict set FracturesDict " << Id << " BotPoint Coordinates $Coordinates" << std::endl;
405  // LeftPoint
406  PropDataFile << "dict set FracturesDict " << Id << " LeftPoint Id "
407  << rParameters["fractures_list"][i]["left_point"]["id"].GetInt() << std::endl;
408  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["left_point"]["coordinates"][0].GetDouble()
409  << " " << rParameters["fractures_list"][i]["left_point"]["coordinates"][1].GetDouble() << " "
410  << rParameters["fractures_list"][i]["left_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
411  PropDataFile << "dict set FracturesDict " << Id << " LeftPoint Coordinates $Coordinates" << std::endl;
412  // RightPoint
413  PropDataFile << "dict set FracturesDict " << Id << " RightPoint Id "
414  << rParameters["fractures_list"][i]["right_point"]["id"].GetInt() << std::endl;
415  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["right_point"]["coordinates"][0].GetDouble()
416  << " " << rParameters["fractures_list"][i]["right_point"]["coordinates"][1].GetDouble() << " "
417  << rParameters["fractures_list"][i]["right_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
418  PropDataFile << "dict set FracturesDict " << Id << " RightPoint Coordinates $Coordinates" << std::endl;
419  // TopLine
420  PropDataFile << "dict set FracturesDict " << Id << " TopLine Id "
421  << rParameters["fractures_list"][i]["top_line"]["id"].GetInt() << std::endl;
422  // BotLine
423  PropDataFile << "dict set FracturesDict " << Id << " BotLine Id "
424  << rParameters["fractures_list"][i]["bot_line"]["id"].GetInt() << std::endl;
425  // LeftLine
426  PropDataFile << "dict set FracturesDict " << Id << " LeftLine Id "
427  << rParameters["fractures_list"][i]["left_line"]["id"].GetInt() << std::endl;
428  // RightLine
429  PropDataFile << "dict set FracturesDict " << Id << " RightLine Id "
430  << rParameters["fractures_list"][i]["right_line"]["id"].GetInt() << std::endl;
431  // TopLeftSurface
432  PropDataFile << "dict set FracturesDict " << Id << " TopLeftSurface Id "
433  << rParameters["fractures_list"][i]["top_left_surface"]["id"].GetInt() << std::endl;
434  PropDataFile << "set Lines [list]" << std::endl;
435  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["top_left_surface"]["lines"].size(); j++)
436  {
437  PropDataFile << "lappend Lines " << rParameters["fractures_list"][i]["top_left_surface"]["lines"][j].GetInt() << std::endl;
438  }
439  PropDataFile << "dict set FracturesDict " << Id << " TopLeftSurface Lines $Lines" << std::endl;
440  // TopRightSurface
441  PropDataFile << "dict set FracturesDict " << Id << " TopRightSurface Id "
442  << rParameters["fractures_list"][i]["top_right_surface"]["id"].GetInt() << std::endl;
443  PropDataFile << "set Lines [list]" << std::endl;
444  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["top_right_surface"]["lines"].size(); j++)
445  {
446  PropDataFile << "lappend Lines " << rParameters["fractures_list"][i]["top_right_surface"]["lines"][j].GetInt() << std::endl;
447  }
448  PropDataFile << "dict set FracturesDict " << Id << " TopRightSurface Lines $Lines" << std::endl;
449  // BotLeftSurface
450  PropDataFile << "dict set FracturesDict " << Id << " BotLeftSurface Id "
451  << rParameters["fractures_list"][i]["bot_left_surface"]["id"].GetInt() << std::endl;
452  PropDataFile << "set Lines [list]" << std::endl;
453  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["bot_left_surface"]["lines"].size(); j++)
454  {
455  PropDataFile << "lappend Lines " << rParameters["fractures_list"][i]["bot_left_surface"]["lines"][j].GetInt() << std::endl;
456  }
457  PropDataFile << "dict set FracturesDict " << Id << " BotLeftSurface Lines $Lines" << std::endl;
458  // BotRightSurface
459  PropDataFile << "dict set FracturesDict " << Id << " BotRightSurface Id "
460  << rParameters["fractures_list"][i]["bot_right_surface"]["id"].GetInt() << std::endl;
461  PropDataFile << "set Lines [list]" << std::endl;
462  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["bot_right_surface"]["lines"].size(); j++)
463  {
464  PropDataFile << "lappend Lines " << rParameters["fractures_list"][i]["bot_right_surface"]["lines"][j].GetInt() << std::endl;
465  }
466  PropDataFile << "dict set FracturesDict " << Id << " BotRightSurface Lines $Lines" << std::endl;
467  // LeftInterfaceVolume
468  PropDataFile << "dict set FracturesDict " << Id << " LeftInterfaceVolume Id "
469  << rParameters["fractures_list"][i]["left_interface_volume"]["id"].GetInt() << std::endl;
470  PropDataFile << "dict set FracturesDict " << Id << " LeftInterfaceVolume Layer \""
471  << rParameters["fractures_list"][i]["left_interface_volume"]["layer"].GetString() << "\"" << std::endl;
472  PropDataFile << "set Groups [list]" << std::endl;
473  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["left_interface_volume"]["groups"].size(); j++)
474  {
475  PropDataFile << "lappend Groups \"" << rParameters["fractures_list"][i]["left_interface_volume"]["groups"][j].GetString() << "\"" << std::endl;
476  }
477  PropDataFile << "dict set FracturesDict " << Id << " LeftInterfaceVolume Groups $Groups" << std::endl;
478  // RightInterfaceVolume
479  PropDataFile << "dict set FracturesDict " << Id << " RightInterfaceVolume Id "
480  << rParameters["fractures_list"][i]["right_interface_volume"]["id"].GetInt() << std::endl;
481  PropDataFile << "dict set FracturesDict " << Id << " RightInterfaceVolume Layer \""
482  << rParameters["fractures_list"][i]["right_interface_volume"]["layer"].GetString() << "\"" << std::endl;
483  PropDataFile << "set Groups [list]" << std::endl;
484  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["right_interface_volume"]["groups"].size(); j++)
485  {
486  PropDataFile << "lappend Groups \"" << rParameters["fractures_list"][i]["right_interface_volume"]["groups"][j].GetString() << "\"" << std::endl;
487  }
488  PropDataFile << "dict set FracturesDict " << Id << " RightInterfaceVolume Groups $Groups" << std::endl;
489  // BodyVolumes
490  PropDataFile << "set BodyVolumes [list]" << std::endl;
491  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["body_volumes"].size(); j++)
492  {
493  PropDataFile << "lappend BodyVolumes " << rParameters["fractures_list"][i]["body_volumes"][j].GetInt() << std::endl;
494  }
495  PropDataFile << "dict set FracturesDict " << Id << " BodyVolumes $BodyVolumes" << std::endl;
496  }
497  PropDataFile << "lappend PropagationData $FracturesDict" << std::endl;
498 
499  // PropagationDict
500  PropDataFile << "set PropagationDict [dict create]" << std::endl;
501  for(unsigned int i = 0; i < PropagationData.PropagationVector.size(); i++)
502  {
503  PropDataFile << "dict set PropagationDict " << i << " MotherFractureId "
504  << PropagationData.PropagationVector[i].MotherFractureId << std::endl;
505 
506  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].TopInitCoordinates[0]
507  << " " << PropagationData.PropagationVector[i].TopInitCoordinates[1] << " "
508  << PropagationData.PropagationVector[i].TopInitCoordinates[2] << "\"" << std::endl;
509  PropDataFile << "dict set PropagationDict " << i << " TopInitCoordinates $Coordinates" << std::endl;
510  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].BotInitCoordinates[0]
511  << " " << PropagationData.PropagationVector[i].BotInitCoordinates[1] << " "
512  << PropagationData.PropagationVector[i].BotInitCoordinates[2] << "\"" << std::endl;
513  PropDataFile << "dict set PropagationDict " << i << " BotInitCoordinates $Coordinates" << std::endl;
514  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].LeftInitCoordinates[0]
515  << " " << PropagationData.PropagationVector[i].LeftInitCoordinates[1] << " "
516  << PropagationData.PropagationVector[i].LeftInitCoordinates[2] << "\"" << std::endl;
517  PropDataFile << "dict set PropagationDict " << i << " LeftInitCoordinates $Coordinates" << std::endl;
518  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].RightInitCoordinates[0]
519  << " " << PropagationData.PropagationVector[i].RightInitCoordinates[1] << " "
520  << PropagationData.PropagationVector[i].RightInitCoordinates[2] << "\"" << std::endl;
521  PropDataFile << "dict set PropagationDict " << i << " RightInitCoordinates $Coordinates" << std::endl;
522  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].TopEndCoordinates[0]
523  << " " << PropagationData.PropagationVector[i].TopEndCoordinates[1] << " "
524  << PropagationData.PropagationVector[i].TopEndCoordinates[2] << "\"" << std::endl;
525  PropDataFile << "dict set PropagationDict " << i << " TopEndCoordinates $Coordinates" << std::endl;
526  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].BotEndCoordinates[0]
527  << " " << PropagationData.PropagationVector[i].BotEndCoordinates[1] << " "
528  << PropagationData.PropagationVector[i].BotEndCoordinates[2] << "\"" << std::endl;
529  PropDataFile << "dict set PropagationDict " << i << " BotEndCoordinates $Coordinates" << std::endl;
530  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].LeftEndCoordinates[0]
531  << " " << PropagationData.PropagationVector[i].LeftEndCoordinates[1] << " "
532  << PropagationData.PropagationVector[i].LeftEndCoordinates[2] << "\"" << std::endl;
533  PropDataFile << "dict set PropagationDict " << i << " LeftEndCoordinates $Coordinates" << std::endl;
534  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].RightEndCoordinates[0]
535  << " " << PropagationData.PropagationVector[i].RightEndCoordinates[1] << " "
536  << PropagationData.PropagationVector[i].RightEndCoordinates[2] << "\"" << std::endl;
537  PropDataFile << "dict set PropagationDict " << i << " RightEndCoordinates $Coordinates" << std::endl;
538  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].TipCoordinates[0]
539  << " " << PropagationData.PropagationVector[i].TipCoordinates[1] << " "
540  << PropagationData.PropagationVector[i].TipCoordinates[2] << "\"" << std::endl;
541  PropDataFile << "dict set PropagationDict " << i << " TipCoordinates $Coordinates" << std::endl;
542  }
543  PropDataFile << "lappend PropagationData $PropagationDict" << std::endl;
544 
545  // BifurcationDict
546  PropDataFile << "set BifurcationDict [dict create]" << std::endl;
547  for(unsigned int i = 0; i < PropagationData.BifurcationVector.size(); i++)
548  {
549  PropDataFile << "dict set BifurcationDict " << i << " MotherFractureId "
550  << PropagationData.BifurcationVector[i].MotherFractureId << std::endl;
551 
552  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopInitCoordinates[0]
553  << " " << PropagationData.BifurcationVector[i].TopInitCoordinates[1] << " "
554  << PropagationData.BifurcationVector[i].TopInitCoordinates[2] << "\"" << std::endl;
555  PropDataFile << "dict set BifurcationDict " << i << " TopInitCoordinates $Coordinates" << std::endl;
556  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotInitCoordinates[0]
557  << " " << PropagationData.BifurcationVector[i].BotInitCoordinates[1] << " "
558  << PropagationData.BifurcationVector[i].BotInitCoordinates[2] << "\"" << std::endl;
559  PropDataFile << "dict set BifurcationDict " << i << " BotInitCoordinates $Coordinates" << std::endl;
560  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].LeftInitCoordinates[0]
561  << " " << PropagationData.BifurcationVector[i].LeftInitCoordinates[1] << " "
562  << PropagationData.BifurcationVector[i].LeftInitCoordinates[2] << "\"" << std::endl;
563  PropDataFile << "dict set BifurcationDict " << i << " LeftInitCoordinates $Coordinates" << std::endl;
564  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].RightInitCoordinates[0]
565  << " " << PropagationData.BifurcationVector[i].RightInitCoordinates[1] << " "
566  << PropagationData.BifurcationVector[i].RightInitCoordinates[2] << "\"" << std::endl;
567  PropDataFile << "dict set BifurcationDict " << i << " RightInitCoordinates $Coordinates" << std::endl;
568  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopTopEndCoordinates[0]
569  << " " << PropagationData.BifurcationVector[i].TopTopEndCoordinates[1] << " "
570  << PropagationData.BifurcationVector[i].TopTopEndCoordinates[2] << "\"" << std::endl;
571  PropDataFile << "dict set BifurcationDict " << i << " TopTopEndCoordinates $Coordinates" << std::endl;
572  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopBotEndCoordinates[0]
573  << " " << PropagationData.BifurcationVector[i].TopBotEndCoordinates[1] << " "
574  << PropagationData.BifurcationVector[i].TopBotEndCoordinates[2] << "\"" << std::endl;
575  PropDataFile << "dict set BifurcationDict " << i << " TopBotEndCoordinates $Coordinates" << std::endl;
576  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopLeftEndCoordinates[0]
577  << " " << PropagationData.BifurcationVector[i].TopLeftEndCoordinates[1] << " "
578  << PropagationData.BifurcationVector[i].TopLeftEndCoordinates[2] << "\"" << std::endl;
579  PropDataFile << "dict set BifurcationDict " << i << " TopLeftEndCoordinates $Coordinates" << std::endl;
580  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopRightEndCoordinates[0]
581  << " " << PropagationData.BifurcationVector[i].TopRightEndCoordinates[1] << " "
582  << PropagationData.BifurcationVector[i].TopRightEndCoordinates[2] << "\"" << std::endl;
583  PropDataFile << "dict set BifurcationDict " << i << " TopRightEndCoordinates $Coordinates" << std::endl;
584  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopTipCoordinates[0]
585  << " " << PropagationData.BifurcationVector[i].TopTipCoordinates[1] << " "
586  << PropagationData.BifurcationVector[i].TopTipCoordinates[2] << "\"" << std::endl;
587  PropDataFile << "dict set BifurcationDict " << i << " TopTipCoordinates $Coordinates" << std::endl;
588  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotTopEndCoordinates[0]
589  << " " << PropagationData.BifurcationVector[i].BotTopEndCoordinates[1] << " "
590  << PropagationData.BifurcationVector[i].BotTopEndCoordinates[2] << "\"" << std::endl;
591  PropDataFile << "dict set BifurcationDict " << i << " BotTopEndCoordinates $Coordinates" << std::endl;
592  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotBotEndCoordinates[0]
593  << " " << PropagationData.BifurcationVector[i].BotBotEndCoordinates[1] << " "
594  << PropagationData.BifurcationVector[i].BotBotEndCoordinates[2] << "\"" << std::endl;
595  PropDataFile << "dict set BifurcationDict " << i << " BotBotEndCoordinates $Coordinates" << std::endl;
596  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotLeftEndCoordinates[0]
597  << " " << PropagationData.BifurcationVector[i].BotLeftEndCoordinates[1] << " "
598  << PropagationData.BifurcationVector[i].BotLeftEndCoordinates[2] << "\"" << std::endl;
599  PropDataFile << "dict set BifurcationDict " << i << " BotLeftEndCoordinates $Coordinates" << std::endl;
600  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotRightEndCoordinates[0]
601  << " " << PropagationData.BifurcationVector[i].BotRightEndCoordinates[1] << " "
602  << PropagationData.BifurcationVector[i].BotRightEndCoordinates[2] << "\"" << std::endl;
603  PropDataFile << "dict set BifurcationDict " << i << " BotRightEndCoordinates $Coordinates" << std::endl;
604  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotTipCoordinates[0]
605  << " " << PropagationData.BifurcationVector[i].BotTipCoordinates[1] << " "
606  << PropagationData.BifurcationVector[i].BotTipCoordinates[2] << "\"" << std::endl;
607  PropDataFile << "dict set BifurcationDict " << i << " BotTipCoordinates $Coordinates" << std::endl;
608  }
609  PropDataFile << "lappend PropagationData $BifurcationDict" << std::endl;
610 
611  PropDataFile << "return $PropagationData" << std::endl;
612 
613  PropDataFile.close();
614  }
615 
617 
618  void InitializeMapping(
619  UtilityVariables& rAuxVariables,
620  ModelPart& rModelPartNew,
621  const bool& move_mesh_flag)
622  {
623  // Move mesh to the original position to work in the reference state
624  if(move_mesh_flag==true)
625  this->ResetMeshPosition(rModelPartNew);
626 
627  this->ComputeCellMatrixDimensions(rAuxVariables,rModelPartNew);
628  }
629 
631 
632  void NodalVariablesMapping(
633  const UtilityVariables& AuxVariables,
634  Parameters& rParameters,
635  ModelPart& rModelPartOld,
636  ModelPart& rModelPartNew)
637  {
638  // Define ElementOld Cell matrix
639  std::vector< std::vector< std::vector< std::vector<Element::Pointer> > > > ElementOldCellMatrix;
640  ElementOldCellMatrix.resize(AuxVariables.NRows);
641  for(int i = 0; i < AuxVariables.NRows; i++)
642  ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
643  for(int i = 0; i < AuxVariables.NRows; i++)
644  {
645  for(int j = 0; j < AuxVariables.NColumns; j++)
646  {
647  ElementOldCellMatrix[i][j].resize(AuxVariables.NSections);
648  }
649  }
650 
651  // Locate Old Elments in cells
652  double X_me;
653  double Y_me;
654  double Z_me;
655  int PointsNumber;
656 
657  unsigned int NumBodySubModelParts = rParameters["fracture_data"]["body_domain_sub_model_part_list"].size();
658 
659  // Loop through all BodySubModelParts
660  for(unsigned int m = 0; m < NumBodySubModelParts; m++)
661  {
662  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][m].GetString());
663 
664  int NElems = static_cast<int>(SubModelPart.Elements().size());
665  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
666 
667  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber)
668  for(int i = 0; i < NElems; i++)
669  {
670  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
671 
672  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
673  double X_right = X_left;
674  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
675  double Y_bot = Y_top;
676  double Z_back = itElemOld->GetGeometry().GetPoint(0).Z0();
677  double Z_front = Z_back;
678  PointsNumber = itElemOld->GetGeometry().PointsNumber();
679 
680  for(int j = 1; j < PointsNumber; j++)
681  {
682  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
683  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
684  Z_me = itElemOld->GetGeometry().GetPoint(j).Z0();
685 
686  if(X_me > X_right) X_right = X_me;
687  else if(X_me < X_left) X_left = X_me;
688  if(Y_me > Y_top) Y_top = Y_me;
689  else if(Y_me < Y_bot) Y_bot = Y_me;
690  if(Z_me > Z_front) Z_front = Z_me;
691  else if(Z_me < Z_back) Z_back = Z_me;
692  }
693 
694  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
695  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
696  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
697  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
698  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
699  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
700 
701  if(Column_left < 0) Column_left = 0;
702  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
703  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
704  else if(Column_right < 0) Column_right = 0;
705  if(Row_top < 0) Row_top = 0;
706  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
707  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
708  else if(Row_bot < 0) Row_bot = 0;
709  if(Section_back < 0) Section_back = 0;
710  else if(Section_back >= AuxVariables.NSections) Section_back = AuxVariables.NSections-1;
711  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
712  else if(Section_front < 0) Section_front = 0;
713 
714  for(int s = Section_back; s <= Section_front; s++)
715  {
716  for(int k = Row_top; k <= Row_bot; k++)
717  {
718  for(int l = Column_left; l <= Column_right; l++)
719  {
720  #pragma omp critical
721  {
722  ElementOldCellMatrix[k][l][s].push_back((*(itElemOld.base())));
723  }
724  }
725  }
726  }
727  }
728  }
729 
730  unsigned int NumInterfaceSubModelPartsOld = rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"].size();
731 
732  // Loop through all InterfaceSubModelParts
733  for(unsigned int m = 0; m < NumInterfaceSubModelPartsOld; m++)
734  {
735  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"][m].GetString());
736 
737  int NElems = static_cast<int>(SubModelPart.Elements().size());
738  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
739 
740  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber)
741  for(int i = 0; i < NElems; i++)
742  {
743  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
744 
745  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
746  double X_right = X_left;
747  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
748  double Y_bot = Y_top;
749  double Z_back = itElemOld->GetGeometry().GetPoint(0).Z0();
750  double Z_front = Z_back;
751  PointsNumber = itElemOld->GetGeometry().PointsNumber();
752 
753  for(int j = 1; j < PointsNumber; j++)
754  {
755  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
756  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
757  Z_me = itElemOld->GetGeometry().GetPoint(j).Z0();
758 
759  if(X_me > X_right) X_right = X_me;
760  else if(X_me < X_left) X_left = X_me;
761  if(Y_me > Y_top) Y_top = Y_me;
762  else if(Y_me < Y_bot) Y_bot = Y_me;
763  if(Z_me > Z_front) Z_front = Z_me;
764  else if(Z_me < Z_back) Z_back = Z_me;
765  }
766 
767  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
768  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
769  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
770  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
771  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
772  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
773 
774  if(Column_left < 0) Column_left = 0;
775  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
776  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
777  else if(Column_right < 0) Column_right = 0;
778  if(Row_top < 0) Row_top = 0;
779  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
780  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
781  else if(Row_bot < 0) Row_bot = 0;
782  if(Section_back < 0) Section_back = 0;
783  else if(Section_back >= AuxVariables.NSections) Section_back = AuxVariables.NSections-1;
784  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
785  else if(Section_front < 0) Section_front = 0;
786 
787  for(int s = Section_back; s <= Section_front; s++)
788  {
789  for(int k = Row_top; k <= Row_bot; k++)
790  {
791  for(int l = Column_left; l <= Column_right; l++)
792  {
793  #pragma omp critical
794  {
795  ElementOldCellMatrix[k][l][s].push_back((*(itElemOld.base())));
796  }
797  }
798  }
799  }
800  }
801  }
802 
803  // Locate new nodes inside old elements and interpolate nodal variables
804  const int NNodes = static_cast<int>(rModelPartNew.Nodes().size());
805  ModelPart::NodesContainerType::iterator node_begin = rModelPartNew.NodesBegin();
806 
807  array_1d<double,3> GlobalCoordinates;
808  array_1d<double,3> LocalCoordinates;
809 
810  #pragma omp parallel for private(X_me,Y_me,Z_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
811  for(int i = 0; i < NNodes; i++)
812  {
813  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
814 
815  X_me = itNodeNew->X0();
816  Y_me = itNodeNew->Y0();
817  Z_me = itNodeNew->Z0();
818 
819  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
820  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
821  int Section = int((Z_me-AuxVariables.Z_min)/AuxVariables.SectionSize);
822 
823  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
824  else if(Column < 0) Column = 0;
825  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
826  else if(Row < 0) Row = 0;
827  if(Section >= AuxVariables.NSections) Section = AuxVariables.NSections-1;
828  else if(Section < 0) Section = 0;
829 
830  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
831  bool IsInside = false;
832  Element::Pointer pElementOld;
833 
834  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column][Section]).size(); m++)
835  {
836  pElementOld = ElementOldCellMatrix[Row][Column][Section][m];
837  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
838  if(IsInside) break;
839  }
840  if(IsInside==false)
841  {
842  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column][Section]).size(); m++)
843  {
844  pElementOld = ElementOldCellMatrix[Row][Column][Section][m];
845  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,1.0e-5);
846  if(IsInside) break;
847  }
848  }
849  if(IsInside == false)
850  std::cout << "ERROR!!, NONE OF THE OLD ELEMENTS CONTAINS NODE: " << itNodeNew->Id() << std::endl;
851 
852  PointsNumber = pElementOld->GetGeometry().PointsNumber();
853  Vector ShapeFunctionsValuesVector(PointsNumber);
854  Vector NodalVariableVector(PointsNumber);
855 
856  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
857 
858  // Interpolation of nodal variables
859  if( itNodeNew->IsFixed(DISPLACEMENT_X)==false )
860  {
861  for(int j = 0; j < PointsNumber; j++)
862  {
863  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_X);
864  }
865  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
866  }
867  if( itNodeNew->IsFixed(VELOCITY_X)==false )
868  {
869  for(int j = 0; j < PointsNumber; j++)
870  {
871  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(VELOCITY_X);
872  }
873  itNodeNew->FastGetSolutionStepValue(VELOCITY_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
874  }
875  if( itNodeNew->IsFixed(ACCELERATION_X)==false )
876  {
877  for(int j = 0; j < PointsNumber; j++)
878  {
879  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(ACCELERATION_X);
880  }
881  itNodeNew->FastGetSolutionStepValue(ACCELERATION_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
882  }
883  if( itNodeNew->IsFixed(DISPLACEMENT_Y)==false )
884  {
885  for(int j = 0; j < PointsNumber; j++)
886  {
887  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_Y);
888  }
889  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
890  }
891  if( itNodeNew->IsFixed(VELOCITY_Y)==false )
892  {
893  for(int j = 0; j < PointsNumber; j++)
894  {
895  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(VELOCITY_Y);
896  }
897  itNodeNew->FastGetSolutionStepValue(VELOCITY_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
898  }
899  if( itNodeNew->IsFixed(ACCELERATION_Y)==false )
900  {
901  for(int j = 0; j < PointsNumber; j++)
902  {
903  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(ACCELERATION_Y);
904  }
905  itNodeNew->FastGetSolutionStepValue(ACCELERATION_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
906  }
907  if( itNodeNew->IsFixed(DISPLACEMENT_Z)==false )
908  {
909  for(int j = 0; j < PointsNumber; j++)
910  {
911  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_Z);
912  }
913  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Z) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
914  }
915  if( itNodeNew->IsFixed(VELOCITY_Z)==false )
916  {
917  for(int j = 0; j < PointsNumber; j++)
918  {
919  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(VELOCITY_Z);
920  }
921  itNodeNew->FastGetSolutionStepValue(VELOCITY_Z) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
922  }
923  if( itNodeNew->IsFixed(ACCELERATION_Z)==false )
924  {
925  for(int j = 0; j < PointsNumber; j++)
926  {
927  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(ACCELERATION_Z);
928  }
929  itNodeNew->FastGetSolutionStepValue(ACCELERATION_Z) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
930  }
931  if( itNodeNew->IsFixed(WATER_PRESSURE)==false )
932  {
933  for(int j = 0; j < PointsNumber; j++)
934  {
935  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(WATER_PRESSURE);
936  }
937  itNodeNew->FastGetSolutionStepValue(WATER_PRESSURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
938  for(int j = 0; j < PointsNumber; j++)
939  {
940  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DT_WATER_PRESSURE);
941  }
942  itNodeNew->FastGetSolutionStepValue(DT_WATER_PRESSURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
943  }
944  }
945  }
946 
948 
949  void GaussPointStateVariableMapping(
950  const UtilityVariables& AuxVariables,
951  Parameters& rParameters,
952  ModelPart& rModelPartOld,
953  ModelPart& rModelPartNew)
954  {
955  // Define GaussPointOld Cell matrix
956  std::vector< std::vector< std::vector< std::vector<GaussPointOld> > > > BodyGaussPointOldCellMatrix;
957  BodyGaussPointOldCellMatrix.resize(AuxVariables.NRows);
958  for(int i = 0; i < AuxVariables.NRows; i++)
959  BodyGaussPointOldCellMatrix[i].resize(AuxVariables.NColumns);
960  for(int i = 0; i < AuxVariables.NRows; i++)
961  {
962  for(int j = 0; j < AuxVariables.NColumns; j++)
963  {
964  BodyGaussPointOldCellMatrix[i][j].resize(AuxVariables.NSections);
965  }
966  }
967 
968  std::vector< std::vector< std::vector< std::vector<GaussPointOld> > > > InterfaceGaussPointOldCellMatrix;
969  InterfaceGaussPointOldCellMatrix.resize(AuxVariables.NRows);
970  for(int i = 0; i < AuxVariables.NRows; i++)
971  InterfaceGaussPointOldCellMatrix[i].resize(AuxVariables.NColumns);
972  for(int i = 0; i < AuxVariables.NRows; i++)
973  {
974  for(int j = 0; j < AuxVariables.NColumns; j++)
975  {
976  InterfaceGaussPointOldCellMatrix[i][j].resize(AuxVariables.NSections);
977  }
978  }
979 
980  // Locate Old Gauss Points in cells
981  GaussPointOld MyGaussPointOld;
982  GeometryData::IntegrationMethod MyIntegrationMethod;
983  const ProcessInfo& CurrentProcessInfoOld = rModelPartOld.GetProcessInfo();
984  array_1d<double,3> AuxLocalCoordinates;
985 
986  unsigned int NumBodySubModelParts = rParameters["fracture_data"]["body_domain_sub_model_part_list"].size();
987 
988  // Loop through all OLD BodySubModelParts
989  for(unsigned int i = 0; i < NumBodySubModelParts; i++)
990  {
991  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][i].GetString());
992 
993  int NElems = static_cast<int>(SubModelPart.Elements().size());
994  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
995 
996  #pragma omp parallel for private(MyGaussPointOld,MyIntegrationMethod,AuxLocalCoordinates)
997  for(int j = 0; j < NElems; j++)
998  {
999  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
1000 
1001  Element::GeometryType& rGeom = itElem->GetGeometry();
1002  MyIntegrationMethod = itElem->GetIntegrationMethod();
1003  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
1004  unsigned int NumGPoints = IntegrationPoints.size();
1005  Vector detJContainer(NumGPoints);
1006  rGeom.DeterminantOfJacobian(detJContainer,MyIntegrationMethod);
1007  std::vector<double> StateVariableVector(NumGPoints);
1008  itElem->CalculateOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoOld);
1009  int Row;
1010  int Column;
1011  int Section;
1012 
1013  // Loop through GaussPoints
1014  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
1015  {
1016  // GaussPointOld Coordinates
1017  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
1018  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
1019  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
1020  rGeom.GlobalCoordinates(MyGaussPointOld.Coordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
1021 
1022  // GaussPointOld Weight
1023  MyGaussPointOld.Weight = detJContainer[GPoint]*IntegrationPoints[GPoint].Weight();
1024 
1025  // GaussPointOld StateVariable and Damage
1026  MyGaussPointOld.StateVariable = StateVariableVector[GPoint];
1027 
1028  // GaussPointOld Row and Column
1029  Row = int((AuxVariables.Y_max-MyGaussPointOld.Coordinates[1])/AuxVariables.RowSize);
1030  Column = int((MyGaussPointOld.Coordinates[0]-AuxVariables.X_min)/AuxVariables.ColumnSize);
1031  Section = int((MyGaussPointOld.Coordinates[2]-AuxVariables.Z_min)/AuxVariables.SectionSize);
1032  #pragma omp critical
1033  {
1034  BodyGaussPointOldCellMatrix[Row][Column][Section].push_back(MyGaussPointOld);
1035  }
1036  }
1037  }
1038  }
1039 
1040  unsigned int NumInterfaceSubModelPartsOld = rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"].size();
1041 
1042  // Loop through all OLD InterfaceSubModelParts
1043  for(unsigned int i = 0; i < NumInterfaceSubModelPartsOld; i++)
1044  {
1045  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"][i].GetString());
1046 
1047  int NElems = static_cast<int>(SubModelPart.Elements().size());
1048  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
1049 
1050  #pragma omp parallel for private(MyGaussPointOld,MyIntegrationMethod,AuxLocalCoordinates)
1051  for(int j = 0; j < NElems; j++)
1052  {
1053  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
1054 
1055  Element::GeometryType& rGeom = itElem->GetGeometry();
1056  MyIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1;
1057  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
1058  unsigned int NumGPoints = IntegrationPoints.size();
1059  Vector detJContainer(NumGPoints);
1060  rGeom.DeterminantOfJacobian(detJContainer,MyIntegrationMethod);
1061  std::vector<double> StateVariableVector(NumGPoints);
1062  itElem->CalculateOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoOld);
1063  int Row;
1064  int Column;
1065  int Section;
1066 
1067  // Loop through GaussPoints
1068  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
1069  {
1070  // GaussPointOld Coordinates
1071  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
1072  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
1073  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
1074  rGeom.GlobalCoordinates(MyGaussPointOld.Coordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
1075 
1076  // GaussPointOld Weight
1077  MyGaussPointOld.Weight = detJContainer[GPoint]*IntegrationPoints[GPoint].Weight();
1078 
1079  // GaussPointOld StateVariable
1080  MyGaussPointOld.StateVariable = StateVariableVector[GPoint];
1081 
1082  // GaussPointOld Row and Column
1083  Row = int((AuxVariables.Y_max-MyGaussPointOld.Coordinates[1])/AuxVariables.RowSize);
1084  Column = int((MyGaussPointOld.Coordinates[0]-AuxVariables.X_min)/AuxVariables.ColumnSize);
1085  Section = int((MyGaussPointOld.Coordinates[2]-AuxVariables.Z_min)/AuxVariables.SectionSize);
1086  #pragma omp critical
1087  {
1088  InterfaceGaussPointOldCellMatrix[Row][Column][Section].push_back(MyGaussPointOld);
1089  }
1090  }
1091  }
1092  }
1093 
1094  // Transfer state variables from old Gauss points to new Gauss Points (nonlocal average)
1095  const ProcessInfo& CurrentProcessInfoNew = rModelPartNew.GetProcessInfo();
1096  const double PropagationLength = rParameters["fracture_data"]["propagation_length"].GetDouble();
1097  array_1d<double,3> AuxGlobalCoordinates;
1098 
1099  // Loop through all NEW BodySubModelParts
1100  for(unsigned int i = 0; i < NumBodySubModelParts; i++)
1101  {
1102  ModelPart& SubModelPart = rModelPartNew.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][i].GetString());
1103 
1104  int NElems = static_cast<int>(SubModelPart.Elements().size());
1105  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
1106 
1107  double DamageThreshold = el_begin->GetProperties()[DAMAGE_THRESHOLD];
1108 
1109  #pragma omp parallel for private(MyIntegrationMethod,AuxLocalCoordinates,AuxGlobalCoordinates)
1110  for(int j = 0; j < NElems; j++)
1111  {
1112  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
1113 
1114  Element::GeometryType& rGeom = itElem->GetGeometry();
1115  MyIntegrationMethod = itElem->GetIntegrationMethod();
1116  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
1117  unsigned int NumGPoints = IntegrationPoints.size();
1118  std::vector<double> StateVariableVector(NumGPoints);
1119 
1120  // Loop through NEW GaussPoints
1121  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
1122  {
1123  // GaussPointNew Coordinates
1124  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
1125  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
1126  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
1127  rGeom.GlobalCoordinates(AuxGlobalCoordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
1128  double X_me = AuxGlobalCoordinates[0];
1129  double Y_me = AuxGlobalCoordinates[1];
1130  double Z_me = AuxGlobalCoordinates[2];
1131 
1132  // GaussPointNew search area
1133  double X_left = X_me - PropagationLength;
1134  double X_right = X_me + PropagationLength;
1135  double Y_top = Y_me + PropagationLength;
1136  double Y_bot = Y_me - PropagationLength;
1137  double Z_back = Z_me - PropagationLength;
1138  double Z_front = Z_me + PropagationLength;
1139 
1140  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
1141  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
1142  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
1143  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
1144  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
1145  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
1146 
1147  if(Column_left < 0) Column_left = 0;
1148  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
1149  if(Row_top < 0) Row_top = 0;
1150  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
1151  if(Section_back < 0) Section_back = 0;
1152  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
1153 
1154  // Search GaussPointOld neighbours around the GaussPointNew and compute nonlocal state variable
1155  double Numerator = 0.0;
1156  double WeightingFunctionDenominator = 0.0;
1157  double Distance;
1158  for(int s = Section_back; s <= Section_front; s++)
1159  {
1160  for(int k = Row_top; k <= Row_bot; k++)
1161  {
1162  for(int l = Column_left; l<= Column_right; l++)
1163  {
1164  for(unsigned int m = 0; m < BodyGaussPointOldCellMatrix[k][l][s].size(); m++)
1165  {
1166  GaussPointOld& rOtherGaussPointOld = BodyGaussPointOldCellMatrix[k][l][s][m];
1167 
1168  Distance = sqrt((rOtherGaussPointOld.Coordinates[0]-X_me)*(rOtherGaussPointOld.Coordinates[0]-X_me) +
1169  (rOtherGaussPointOld.Coordinates[1]-Y_me)*(rOtherGaussPointOld.Coordinates[1]-Y_me) +
1170  (rOtherGaussPointOld.Coordinates[2]-Z_me)*(rOtherGaussPointOld.Coordinates[2]-Z_me));
1171 
1172  if(Distance <= PropagationLength)
1173  {
1174  Numerator += rOtherGaussPointOld.Weight
1175  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength))
1176  *rOtherGaussPointOld.StateVariable;
1177  WeightingFunctionDenominator += rOtherGaussPointOld.Weight
1178  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength));
1179  }
1180  }
1181  }
1182  }
1183  }
1184  // Save computed stateVariable
1185  if(WeightingFunctionDenominator > 0.0)
1186  StateVariableVector[GPoint] = Numerator/WeightingFunctionDenominator;
1187  else
1188  StateVariableVector[GPoint] = DamageThreshold;
1189  }
1190  // Set stateVariable of new GaussPoints
1191  itElem->SetValuesOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoNew);
1192  }
1193  }
1194 
1195  unsigned int NumInterfaceSubModelParts = rParameters["fracture_data"]["interface_domain_sub_model_part_list"].size();
1196  const double PropagationDamage = rParameters["fracture_data"]["propagation_damage"].GetDouble();
1197 
1198  // Loop through all NEW InterfaceSubModelParts
1199  for(unsigned int i = 0; i < NumInterfaceSubModelParts; i++)
1200  {
1201  ModelPart& SubModelPart = rModelPartNew.GetSubModelPart(rParameters["fracture_data"]["interface_domain_sub_model_part_list"][i].GetString());
1202 
1203  int NElems = static_cast<int>(SubModelPart.Elements().size());
1204  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
1205 
1206  //double DamageThreshold = el_begin->GetProperties()[DAMAGE_THRESHOLD];
1207 
1208  #pragma omp parallel for private(MyIntegrationMethod,AuxLocalCoordinates,AuxGlobalCoordinates)
1209  for(int j = 0; j < NElems; j++)
1210  {
1211  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
1212 
1213  Element::GeometryType& rGeom = itElem->GetGeometry();
1214  MyIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1;
1215  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
1216  unsigned int NumGPoints = IntegrationPoints.size();
1217  std::vector<double> StateVariableVector(NumGPoints);
1218 
1219  // Loop through NEW GaussPoints
1220  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
1221  {
1222  // GaussPointNew Coordinates
1223  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
1224  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
1225  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
1226  rGeom.GlobalCoordinates(AuxGlobalCoordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
1227  double X_me = AuxGlobalCoordinates[0];
1228  double Y_me = AuxGlobalCoordinates[1];
1229  double Z_me = AuxGlobalCoordinates[2];
1230 
1231  // GaussPointNew search area
1232  double X_left = X_me - PropagationLength;
1233  double X_right = X_me + PropagationLength;
1234  double Y_top = Y_me + PropagationLength;
1235  double Y_bot = Y_me - PropagationLength;
1236  double Z_back = Z_me - PropagationLength;
1237  double Z_front = Z_me + PropagationLength;
1238 
1239  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
1240  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
1241  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
1242  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
1243  int Section_back = int((Z_back-AuxVariables.Z_min)/AuxVariables.SectionSize);
1244  int Section_front = int((Z_front-AuxVariables.Z_min)/AuxVariables.SectionSize);
1245 
1246  if(Column_left < 0) Column_left = 0;
1247  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
1248  if(Row_top < 0) Row_top = 0;
1249  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
1250  if(Section_back < 0) Section_back = 0;
1251  if(Section_front >= AuxVariables.NSections) Section_front = AuxVariables.NSections-1;
1252 
1253  // Search GaussPointOld neighbours around the GaussPointNew and compute nonlocal state variable
1254  double Numerator = 0.0;
1255  double WeightingFunctionDenominator = 0.0;
1256  double Distance;
1257  for(int s = Section_back; s <= Section_front; s++)
1258  {
1259  for(int k = Row_top; k <= Row_bot; k++)
1260  {
1261  for(int l = Column_left; l<= Column_right; l++)
1262  {
1263  for(unsigned int m = 0; m < InterfaceGaussPointOldCellMatrix[k][l][s].size(); m++)
1264  {
1265  GaussPointOld& rOtherGaussPointOld = InterfaceGaussPointOldCellMatrix[k][l][s][m];
1266 
1267  Distance = sqrt((rOtherGaussPointOld.Coordinates[0]-X_me)*(rOtherGaussPointOld.Coordinates[0]-X_me) +
1268  (rOtherGaussPointOld.Coordinates[1]-Y_me)*(rOtherGaussPointOld.Coordinates[1]-Y_me) +
1269  (rOtherGaussPointOld.Coordinates[2]-Z_me)*(rOtherGaussPointOld.Coordinates[2]-Z_me));
1270 
1271  if(Distance <= PropagationLength)
1272  {
1273  Numerator += rOtherGaussPointOld.Weight
1274  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength))
1275  *rOtherGaussPointOld.StateVariable;
1276  WeightingFunctionDenominator += rOtherGaussPointOld.Weight
1277  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength));
1278  }
1279  }
1280  }
1281  }
1282  }
1283 
1284  // Save computed stateVariable
1285  if(WeightingFunctionDenominator > 0.0)
1286  {
1287  StateVariableVector[GPoint] = Numerator/WeightingFunctionDenominator;
1288  if(StateVariableVector[GPoint] < PropagationDamage)
1289  {
1290  StateVariableVector[GPoint] = PropagationDamage;
1291  }
1292  }
1293  else
1294  StateVariableVector[GPoint] = PropagationDamage;
1295  }
1296  // Set stateVariable of new GaussPoints
1297  itElem->SetValuesOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoNew);
1298  }
1299  }
1300  }
1301 
1303 
1304  void ResetMeshPosition(ModelPart& rModelPart)
1305  {
1306  // Move mesh to the original position
1307  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
1308  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
1309 
1310  #pragma omp parallel for
1311  for(int i = 0; i < NNodes; i++)
1312  {
1313  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
1314 
1315  itNode->X() = itNode->X0();
1316  itNode->Y() = itNode->Y0();
1317  itNode->Z() = itNode->Z0();
1318  }
1319  }
1320 
1322 
1323  void UpdateMeshPosition(ModelPart& rModelPart)
1324  {
1325  // Move mesh to the current position
1326  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
1327  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
1328 
1329  #pragma omp parallel for
1330  for(int i = 0; i < NNodes; i++)
1331  {
1332  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
1333 
1334  itNode->X() = itNode->X0() + itNode->FastGetSolutionStepValue(DISPLACEMENT_X);
1335  itNode->Y() = itNode->Y0() + itNode->FastGetSolutionStepValue(DISPLACEMENT_Y);
1336  itNode->Z() = itNode->Z0() + itNode->FastGetSolutionStepValue(DISPLACEMENT_Z);
1337  }
1338  }
1339 
1340 private:
1341 
1343 
1344  void ComputeCellMatrixDimensions(
1345  UtilityVariables& rAuxVariables,
1346  ModelPart& rModelPart)
1347  {
1348  // Compute X, Y and Z limits of the current geometry
1349  unsigned int NumThreads = ParallelUtilities::GetNumThreads();
1350  std::vector<double> X_max_partition(NumThreads);
1351  std::vector<double> X_min_partition(NumThreads);
1352  std::vector<double> Y_max_partition(NumThreads);
1353  std::vector<double> Y_min_partition(NumThreads);
1354  std::vector<double> Z_max_partition(NumThreads);
1355  std::vector<double> Z_min_partition(NumThreads);
1356 
1357  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
1358  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
1359 
1360  #pragma omp parallel
1361  {
1362  int k = OpenMPUtils::ThisThread();
1363 
1364  X_max_partition[k] = node_begin->X();
1365  X_min_partition[k] = X_max_partition[k];
1366  Y_max_partition[k] = node_begin->Y();
1367  Y_min_partition[k] = Y_max_partition[k];
1368  Z_max_partition[k] = node_begin->Z();
1369  Z_min_partition[k] = Z_max_partition[k];
1370 
1371  double X_me, Y_me, Z_me;
1372 
1373  #pragma omp for
1374  for(int i = 0; i < NNodes; i++)
1375  {
1376  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
1377 
1378  X_me = itNode->X();
1379  Y_me = itNode->Y();
1380  Z_me = itNode->Z();
1381 
1382  if( X_me > X_max_partition[k] ) X_max_partition[k] = X_me;
1383  else if( X_me < X_min_partition[k] ) X_min_partition[k] = X_me;
1384 
1385  if( Y_me > Y_max_partition[k] ) Y_max_partition[k] = Y_me;
1386  else if( Y_me < Y_min_partition[k] ) Y_min_partition[k] = Y_me;
1387 
1388  if( Z_me > Z_max_partition[k] ) Z_max_partition[k] = Z_me;
1389  else if( Z_me < Z_min_partition[k] ) Z_min_partition[k] = Z_me;
1390  }
1391  }
1392 
1393  rAuxVariables.X_max = X_max_partition[0];
1394  rAuxVariables.X_min = X_min_partition[0];
1395  rAuxVariables.Y_max = Y_max_partition[0];
1396  rAuxVariables.Y_min = Y_min_partition[0];
1397  rAuxVariables.Z_max = Z_max_partition[0];
1398  rAuxVariables.Z_min = Z_min_partition[0];
1399 
1400  for(unsigned int i=1; i < NumThreads; i++)
1401  {
1402  if(X_max_partition[i] > rAuxVariables.X_max) rAuxVariables.X_max = X_max_partition[i];
1403  if(X_min_partition[i] < rAuxVariables.X_min) rAuxVariables.X_min = X_min_partition[i];
1404  if(Y_max_partition[i] > rAuxVariables.Y_max) rAuxVariables.Y_max = Y_max_partition[i];
1405  if(Y_min_partition[i] < rAuxVariables.Y_min) rAuxVariables.Y_min = Y_min_partition[i];
1406  if(Z_max_partition[i] > rAuxVariables.Z_max) rAuxVariables.Z_max = Z_max_partition[i];
1407  if(Z_min_partition[i] < rAuxVariables.Z_min) rAuxVariables.Z_min = Z_min_partition[i];
1408  }
1409 
1410  // Calculate Average Element Length
1411  double AverageElementLength = 0.0;
1412 
1413  int NElems = static_cast<int>(rModelPart.Elements().size());
1414  ModelPart::ElementsContainerType::iterator el_begin = rModelPart.ElementsBegin();
1415 
1416  #pragma omp parallel for reduction(+:AverageElementLength)
1417  for(int i = 0; i < NElems; i++)
1418  {
1419  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
1420 
1421  AverageElementLength += itElem->GetGeometry().Length();
1422  }
1423 
1424  AverageElementLength = AverageElementLength/NElems;
1425 
1426  // Compute FracturePoints CellMatrix dimensions
1427  rAuxVariables.NRows = int((rAuxVariables.Y_max-rAuxVariables.Y_min)/AverageElementLength);
1428  rAuxVariables.NColumns = int((rAuxVariables.X_max-rAuxVariables.X_min)/AverageElementLength);
1429  rAuxVariables.NSections = int((rAuxVariables.Z_max-rAuxVariables.Z_min)/AverageElementLength);
1430  if(rAuxVariables.NRows <= 0) rAuxVariables.NRows = 1;
1431  if(rAuxVariables.NColumns <= 0) rAuxVariables.NColumns = 1;
1432  if(rAuxVariables.NSections <= 0) rAuxVariables.NSections = 1;
1433  rAuxVariables.RowSize = (rAuxVariables.Y_max-rAuxVariables.Y_min)/rAuxVariables.NRows;
1434  rAuxVariables.ColumnSize = (rAuxVariables.X_max-rAuxVariables.X_min)/rAuxVariables.NColumns;
1435  rAuxVariables.SectionSize = (rAuxVariables.Z_max-rAuxVariables.Z_min)/rAuxVariables.NSections;
1436  }
1437 
1439 
1440  void SetFracturePoints(
1441  PropagationGlobalVariables& rPropagationData,
1442  UtilityVariables& rAuxVariables,
1443  Parameters& rParameters,
1444  ModelPart& rModelPart)
1445  {
1446  // Compute FracturePointsCellMatrix dimensions
1447  this->ComputeCellMatrixDimensions(rAuxVariables,rModelPart);
1448 
1449  rPropagationData.FracturePointsCellMatrix.resize(rAuxVariables.NRows);
1450  for(int i = 0; i < rAuxVariables.NRows; i++)
1451  rPropagationData.FracturePointsCellMatrix[i].resize(rAuxVariables.NColumns);
1452  for(int i = 0; i < rAuxVariables.NRows; i++)
1453  {
1454  for(int j = 0; j < rAuxVariables.NColumns; j++)
1455  {
1456  rPropagationData.FracturePointsCellMatrix[i][j].resize(rAuxVariables.NSections);
1457  }
1458  }
1459 
1460  // Locate FracturePoints inside CellMatrix
1461  FracturePoint MyFracturePoint;
1462  GeometryData::IntegrationMethod MyIntegrationMethod;
1463  const ProcessInfo& CurrentProcessInfo = rModelPart.GetProcessInfo();
1464  rPropagationData.pProcessInfo = rModelPart.pGetProcessInfo();
1465  array_1d<double,3> AuxLocalCoordinates;
1466 
1467  unsigned int NumBodySubModelParts = rParameters["fracture_data"]["body_domain_sub_model_part_list"].size();
1468 
1469  // Loop through all BodySubModelParts
1470  for(unsigned int i = 0; i < NumBodySubModelParts; i++)
1471  {
1472  ModelPart& BodySubModelPart = rModelPart.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][i].GetString());
1473 
1474  int NElems = static_cast<int>(BodySubModelPart.Elements().size());
1475  ModelPart::ElementsContainerType::iterator el_begin = BodySubModelPart.ElementsBegin();
1476 
1477  // Loop through all body elements
1478  #pragma omp parallel for private(MyFracturePoint,MyIntegrationMethod,AuxLocalCoordinates)
1479  for(int j = 0; j < NElems; j++)
1480  {
1481  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
1482 
1483  Element::GeometryType& rGeom = itElem->GetGeometry();
1484  MyIntegrationMethod = itElem->GetIntegrationMethod();
1485  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
1486  unsigned int NumGPoints = IntegrationPoints.size();
1487  Vector detJContainer(NumGPoints);
1488  rGeom.DeterminantOfJacobian(detJContainer,MyIntegrationMethod);
1489  std::vector<double> DamageVector(NumGPoints);
1490  itElem->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1491  int Row;
1492  int Column;
1493  int Section;
1494 
1495  // Loop through GaussPoints
1496  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
1497  {
1498  // FracturePoint Coordinates
1499  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
1500  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
1501  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
1502  rGeom.GlobalCoordinates(MyFracturePoint.Coordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
1503 
1504  // FracturePoint Weight
1505  MyFracturePoint.Weight = detJContainer[GPoint]*IntegrationPoints[GPoint].Weight();
1506 
1507  // FracturePoint Damage
1508  MyFracturePoint.Damage = DamageVector[GPoint];
1509 
1510  // FracturePoint Row, Column and Section
1511  Row = int((rAuxVariables.Y_max-MyFracturePoint.Coordinates[1])/rAuxVariables.RowSize);
1512  Column = int((MyFracturePoint.Coordinates[0]-rAuxVariables.X_min)/rAuxVariables.ColumnSize);
1513  Section = int((MyFracturePoint.Coordinates[2]-rAuxVariables.Z_min)/rAuxVariables.SectionSize);
1514 
1515  // Element containing the FracturePoint
1516  MyFracturePoint.pElement = (*(itElem.base()));
1517 
1518  #pragma omp critical
1519  {
1520  rPropagationData.FracturePointsCellMatrix[Row][Column][Section].push_back(MyFracturePoint);
1521  }
1522  }
1523  }
1524  }
1525  }
1526 
1528 
1529  void PropagateFracture(
1530  const unsigned int& itFracture,
1531  PropagationGlobalVariables& rPropagationData,
1532  PropagationLocalVariables& rAuxPropagationVariables,
1533  Parameters& rParameters)
1534  {
1535  // Compute Propagation Coordinates
1536  double TipX = 0.0;
1537  double TipY = 0.0;
1538  double TipZ = 0.0;
1539  double TipDen = 0.0;
1540 
1541  #pragma omp parallel sections reduction(+:TipX,TipY,TipZ,TipDen)
1542  {
1543  #pragma omp section
1544  {
1545  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1546  {
1547  this->AverageTipCoordinates(TipX,TipY,TipZ,TipDen,*(rAuxPropagationVariables.TopFrontFracturePoints[i]));
1548  }
1549  }
1550  #pragma omp section
1551  {
1552  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1553  {
1554  this->AverageTipCoordinates(TipX,TipY,TipZ,TipDen,*(rAuxPropagationVariables.BotFrontFracturePoints[i]));
1555  }
1556  }
1557  }
1558 
1559  array_1d<double,3> AuxArray1;
1560  array_1d<double,3> AuxArray2;
1561  const double PropagationLength = rParameters["fracture_data"]["propagation_length"].GetDouble();
1562  const double PropagationWidth = rParameters["fracture_data"]["propagation_width"].GetDouble();
1563  const double PropagationHeight = rParameters["fracture_data"]["propagation_height"].GetDouble();
1564  int MotherFractureId = rParameters["fractures_list"][itFracture]["id"].GetInt();
1565  Propagation MyPropagation;
1566 
1567  MyPropagation.MotherFractureId = MotherFractureId;
1568 
1569  MyPropagation.TipCoordinates[0] = TipX/TipDen;
1570  MyPropagation.TipCoordinates[1] = TipY/TipDen;
1571  MyPropagation.TipCoordinates[2] = TipZ/TipDen;
1572 
1573  // LeftInitCoordinates
1574  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1575  AuxArray1[1] += 0.5*PropagationHeight;
1576  noalias(MyPropagation.LeftInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1577  // RightInitCoordinates
1578  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1579  AuxArray1[1] -= 0.5*PropagationHeight;
1580  noalias(MyPropagation.RightInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1581  // TopInitCoordinates
1582  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1583  AuxArray1[2] += 0.5*PropagationWidth;
1584  noalias(MyPropagation.TopInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1585  // BotInitCoordinates
1586  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1587  AuxArray1[2] -= 0.5*PropagationWidth;
1588  noalias(MyPropagation.BotInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1589 
1590  // Check straight propagation
1591  const double CorrectionTol = rParameters["fracture_data"]["correction_tolerance"].GetDouble();
1592  noalias(AuxArray2) = MyPropagation.TipCoordinates;
1593  noalias(AuxArray1) = prod(rAuxPropagationVariables.RotationMatrix,AuxArray2);
1594  AuxArray1[1] = rAuxPropagationVariables.TipLocalCoordinates[1];
1595  AuxArray1[2] = rAuxPropagationVariables.TipLocalCoordinates[2];
1596  noalias(AuxArray2) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1597 
1598  double Distance = sqrt((MyPropagation.TipCoordinates[0]-AuxArray2[0])*(MyPropagation.TipCoordinates[0]-AuxArray2[0])+
1599  (MyPropagation.TipCoordinates[1]-AuxArray2[1])*(MyPropagation.TipCoordinates[1]-AuxArray2[1])+
1600  (MyPropagation.TipCoordinates[2]-AuxArray2[2])*(MyPropagation.TipCoordinates[2]-AuxArray2[2]));
1601  if (Distance <= PropagationLength*CorrectionTol)
1602  noalias(MyPropagation.TipCoordinates) = AuxArray2;
1603 
1604  // Check whether new tip falls inside a valid element
1605  array_1d<double,3> LocalCoordinates;
1606  Element::Pointer pElement;
1607  bool IsInside = false;
1608 
1609  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1610  {
1611  pElement = rAuxPropagationVariables.TopFrontFracturePoints[i]->pElement;
1612  IsInside = pElement->GetGeometry().IsInside(MyPropagation.TipCoordinates,LocalCoordinates);
1613  if(IsInside) break;
1614  }
1615  if(IsInside == false)
1616  {
1617  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1618  {
1619  pElement = rAuxPropagationVariables.BotFrontFracturePoints[i]->pElement;
1620  IsInside = pElement->GetGeometry().IsInside(MyPropagation.TipCoordinates,LocalCoordinates);
1621  if(IsInside) break;
1622  }
1623  }
1624 
1625  double PropagationDistance = 0.1*PropagationLength;
1626  const double PropagationDamage = rParameters["fracture_data"]["propagation_damage"].GetDouble();
1627  const ProcessInfo& CurrentProcessInfo = *(rPropagationData.pProcessInfo);
1628 
1629  if (IsInside == true)
1630  {
1631  std::vector<double> DamageVector;
1632  pElement->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1633  unsigned int NumGPoints = DamageVector.size();
1634  double InvNumGPoints = 1.0/static_cast<double>(NumGPoints);
1635  double ElementDamage = 0.0;
1636  for (unsigned int i = 0; i < NumGPoints; i++)
1637  {
1638  ElementDamage += DamageVector[i];
1639  }
1640  ElementDamage *= InvNumGPoints;
1641  if (ElementDamage >= 0.5*PropagationDamage)
1642  {
1643  // Compute new Tip RotationMatrix
1644  this->CalculateNewTipRotationMatrix(rAuxPropagationVariables,MyPropagation.TipCoordinates,
1645  MyPropagation.LeftInitCoordinates,PropagationDistance);
1646  // LeftEndCoordinates
1647  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1648  AuxArray1[1] += 0.5*PropagationHeight;
1649  noalias(MyPropagation.LeftEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1650  // RightEndCoordinates
1651  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1652  AuxArray1[1] -= 0.5*PropagationHeight;
1653  noalias(MyPropagation.RightEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1654  // TopEndCoordinates
1655  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1656  AuxArray1[2] += 0.5*PropagationWidth;
1657  noalias(MyPropagation.TopEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1658  // BotEndCoordinates
1659  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1660  AuxArray1[2] -= 0.5*PropagationWidth;
1661  noalias(MyPropagation.BotEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1662 
1663  rPropagationData.PropagationVector.push_back(MyPropagation);
1664  rPropagationData.PropagateFractures = true;
1665  return;
1666  }
1667  }
1668 
1669  // Being here means that the new tip does not fall inside a valid element. We need to check top and bot fracture points
1670  bool PropagateTop = false;
1671  bool PropagateBot = false;
1672  double TopEndX = 0.0;
1673  double TopEndY = 0.0;
1674  double TopEndZ = 0.0;
1675  double TopEndDen = 0.0;
1676  double BotEndX = 0.0;
1677  double BotEndY = 0.0;
1678  double BotEndZ = 0.0;
1679  double BotEndDen = 0.0;
1680  array_1d<double,3> GlobalCoordinates;
1681 
1682  #pragma omp parallel sections private(GlobalCoordinates,LocalCoordinates,pElement,IsInside)
1683  {
1684  #pragma omp section
1685  {
1686  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1687  {
1688  this->AverageTipCoordinates(TopEndX,TopEndY,TopEndZ,TopEndDen,*(rAuxPropagationVariables.TopFrontFracturePoints[i]));
1689  }
1690  GlobalCoordinates[0] = TopEndX/TopEndDen;
1691  GlobalCoordinates[1] = TopEndY/TopEndDen;
1692  GlobalCoordinates[2] = TopEndZ/TopEndDen;
1693 
1694  // Check whether new tip falls inside a valid element
1695  IsInside = false;
1696  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1697  {
1698  pElement = rAuxPropagationVariables.TopFrontFracturePoints[i]->pElement;
1699  IsInside = pElement->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates);
1700  if(IsInside) break;
1701  }
1702 
1703  if (IsInside == true)
1704  {
1705  std::vector<double> DamageVector;
1706  pElement->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1707  unsigned int NumGPoints = DamageVector.size();
1708  double InvNumGPoints = 1.0/static_cast<double>(NumGPoints);
1709  double ElementDamage = 0.0;
1710  for (unsigned int i = 0; i < NumGPoints; i++)
1711  {
1712  ElementDamage += DamageVector[i];
1713  }
1714  ElementDamage *= InvNumGPoints;
1715  if (ElementDamage >= PropagationDamage)
1716  PropagateTop = true;
1717  }
1718  }
1719  #pragma omp section
1720  {
1721  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1722  {
1723  this->AverageTipCoordinates(BotEndX,BotEndY,BotEndZ,BotEndDen,*(rAuxPropagationVariables.BotFrontFracturePoints[i]));
1724  }
1725  GlobalCoordinates[0] = BotEndX/BotEndDen;
1726  GlobalCoordinates[1] = BotEndY/BotEndDen;
1727  GlobalCoordinates[2] = BotEndZ/BotEndDen;
1728 
1729  // Check whether new tip falls inside a valid element
1730  IsInside = false;
1731  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1732  {
1733  pElement = rAuxPropagationVariables.BotFrontFracturePoints[i]->pElement;
1734  IsInside = pElement->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates);
1735  if(IsInside) break;
1736  }
1737 
1738  if (IsInside == true)
1739  {
1740  std::vector<double> DamageVector;
1741  pElement->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1742  unsigned int NumGPoints = DamageVector.size();
1743  double InvNumGPoints = 1.0/static_cast<double>(NumGPoints);
1744  double ElementDamage = 0.0;
1745  for (unsigned int i = 0; i < NumGPoints; i++)
1746  {
1747  ElementDamage += DamageVector[i];
1748  }
1749  ElementDamage *= InvNumGPoints;
1750  if (ElementDamage >= PropagationDamage)
1751  PropagateBot = true;
1752  }
1753  }
1754  }
1755 
1756  if (PropagateTop == true && PropagateBot == true) // Bifurcation
1757  {
1758  Bifurcation MyBifurcation;
1759  MyBifurcation.MotherFractureId = MotherFractureId;
1760 
1761  MyBifurcation.TopTipCoordinates[0] = TopEndX/TopEndDen;
1762  MyBifurcation.TopTipCoordinates[1] = TopEndY/TopEndDen;
1763  MyBifurcation.TopTipCoordinates[2] = TopEndZ/TopEndDen;
1764 
1765  MyBifurcation.BotTipCoordinates[0] = BotEndX/BotEndDen;
1766  MyBifurcation.BotTipCoordinates[1] = BotEndY/BotEndDen;
1767  MyBifurcation.BotTipCoordinates[2] = BotEndZ/BotEndDen;
1768 
1769  // LeftInitCoordinates
1770  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1771  AuxArray1[0] -= PropagationDistance;
1772  AuxArray1[1] += 0.5*PropagationHeight;
1773  noalias(MyBifurcation.LeftInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1774  // RightInitCoordinates
1775  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1776  AuxArray1[0] -= PropagationDistance;
1777  AuxArray1[1] -= 0.5*PropagationHeight;
1778  noalias(MyBifurcation.RightInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1779  // TopInitCoordinates
1780  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1781  AuxArray1[0] -= PropagationDistance;
1782  AuxArray1[2] += 0.5*PropagationWidth;
1783  noalias(MyBifurcation.TopInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1784  // BotInitCoordinates
1785  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1786  AuxArray1[0] -= PropagationDistance;
1787  AuxArray1[2] -= 0.5*PropagationWidth;
1788  noalias(MyBifurcation.BotInitCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1789 
1790  // Compute new TopTip RotationMatrix
1791  this->CalculateNewTipRotationMatrix(rAuxPropagationVariables,MyBifurcation.TopTipCoordinates,
1792  MyBifurcation.LeftInitCoordinates,PropagationDistance);
1793  // TopLeftEndCoordinates
1794  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1795  AuxArray1[1] += 0.5*PropagationHeight;
1796  noalias(MyBifurcation.TopLeftEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1797  // TopRightEndCoordinates
1798  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1799  AuxArray1[1] -= 0.5*PropagationHeight;
1800  noalias(MyBifurcation.TopRightEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1801  // TopTopEndCoordinates
1802  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1803  AuxArray1[2] += 0.5*PropagationWidth;
1804  noalias(MyBifurcation.TopTopEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1805  // TopBotEndCoordinates
1806  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1807  AuxArray1[2] -= 0.5*PropagationWidth;
1808  noalias(MyBifurcation.TopBotEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1809 
1810  // Compute new BotTip RotationMatrix
1811  this->CalculateNewTipRotationMatrix(rAuxPropagationVariables,MyBifurcation.BotTipCoordinates,
1812  MyBifurcation.LeftInitCoordinates,PropagationDistance);
1813  // BotLeftEndCoordinates
1814  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1815  AuxArray1[1] += 0.5*PropagationHeight;
1816  noalias(MyBifurcation.BotLeftEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1817  // BotRightEndCoordinates
1818  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1819  AuxArray1[1] -= 0.5*PropagationHeight;
1820  noalias(MyBifurcation.BotRightEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1821  // BotTopEndCoordinates
1822  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1823  AuxArray1[2] += 0.5*PropagationWidth;
1824  noalias(MyBifurcation.BotTopEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1825  // BotBotEndCoordinates
1826  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1827  AuxArray1[2] -= 0.5*PropagationWidth;
1828  noalias(MyBifurcation.BotBotEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1829 
1830  rPropagationData.BifurcationVector.push_back(MyBifurcation);
1831  rPropagationData.PropagateFractures = true;
1832  return;
1833  }
1834  else if (PropagateTop == true) // Top Propagation
1835  {
1836  MyPropagation.TipCoordinates[0] = TopEndX/TopEndDen;
1837  MyPropagation.TipCoordinates[1] = TopEndY/TopEndDen;
1838  MyPropagation.TipCoordinates[2] = TopEndZ/TopEndDen;
1839 
1840  // Compute new Tip RotationMatrix
1841  this->CalculateNewTipRotationMatrix(rAuxPropagationVariables,MyPropagation.TipCoordinates,
1842  MyPropagation.LeftInitCoordinates,PropagationDistance);
1843  // LeftEndCoordinates
1844  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1845  AuxArray1[1] += 0.5*PropagationHeight;
1846  noalias(MyPropagation.LeftEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1847  // RightEndCoordinates
1848  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1849  AuxArray1[1] -= 0.5*PropagationHeight;
1850  noalias(MyPropagation.RightEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1851  // TopEndCoordinates
1852  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1853  AuxArray1[2] += 0.5*PropagationWidth;
1854  noalias(MyPropagation.TopEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1855  // BotEndCoordinates
1856  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1857  AuxArray1[2] -= 0.5*PropagationWidth;
1858  noalias(MyPropagation.BotEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1859 
1860  rPropagationData.PropagationVector.push_back(MyPropagation);
1861  rPropagationData.PropagateFractures = true;
1862  return;
1863  }
1864  else if (PropagateBot == true) // Bot Propagation
1865  {
1866  MyPropagation.TipCoordinates[0] = BotEndX/BotEndDen;
1867  MyPropagation.TipCoordinates[1] = BotEndY/BotEndDen;
1868  MyPropagation.TipCoordinates[2] = BotEndZ/BotEndDen;
1869 
1870  // Compute new Tip RotationMatrix
1871  this->CalculateNewTipRotationMatrix(rAuxPropagationVariables,MyPropagation.TipCoordinates,
1872  MyPropagation.LeftInitCoordinates,PropagationDistance);
1873  // LeftEndCoordinates
1874  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1875  AuxArray1[1] += 0.5*PropagationHeight;
1876  noalias(MyPropagation.LeftEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1877  // RightEndCoordinates
1878  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1879  AuxArray1[1] -= 0.5*PropagationHeight;
1880  noalias(MyPropagation.RightEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1881  // TopEndCoordinates
1882  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1883  AuxArray1[2] += 0.5*PropagationWidth;
1884  noalias(MyPropagation.TopEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1885  // BotEndCoordinates
1886  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1887  AuxArray1[2] -= 0.5*PropagationWidth;
1888  noalias(MyPropagation.BotEndCoordinates) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1889 
1890  rPropagationData.PropagationVector.push_back(MyPropagation);
1891  rPropagationData.PropagateFractures = true;
1892  return;
1893  }
1894  }
1895 
1897 
1898  void CalculateTipRotationMatrix(
1899  BoundedMatrix<double,3,3>& rRotationMatrix,
1900  const array_1d<double,3>& TipPointCoordinates,
1901  const array_1d<double,3>& LeftPointCoordinates,
1902  const array_1d<double,3>& RightPointCoordinates)
1903  {
1904  array_1d<double, 3> CenterPoint;
1905  noalias(CenterPoint) = 0.5 * (LeftPointCoordinates + RightPointCoordinates);
1906 
1907  // Unitary vector in local x direction
1908  array_1d<double, 3> Vx;
1909  noalias(Vx) = TipPointCoordinates - CenterPoint;
1910  double inv_norm = 1.0/norm_2(Vx);
1911  Vx[0] *= inv_norm;
1912  Vx[1] *= inv_norm;
1913  Vx[2] *= inv_norm;
1914 
1915  // Vector in local y direction
1916  array_1d<double, 3> Vy;
1917  noalias(Vy) = LeftPointCoordinates - CenterPoint;
1918  // Unitary vector in local z direction
1919  array_1d<double, 3> Vz;
1920  MathUtils<double>::CrossProduct(Vz, Vx, Vy);
1921  inv_norm = 1.0/norm_2(Vz);
1922  Vz[0] *= inv_norm;
1923  Vz[1] *= inv_norm;
1924  Vz[2] *= inv_norm;
1925 
1926  // Unitary vector in local y direction
1927  MathUtils<double>::CrossProduct( Vy, Vz, Vx);
1928 
1929  // Rotation Matrix
1930  rRotationMatrix(0,0) = Vx[0];
1931  rRotationMatrix(0,1) = Vx[1];
1932  rRotationMatrix(0,2) = Vx[2];
1933 
1934  rRotationMatrix(1,0) = Vy[0];
1935  rRotationMatrix(1,1) = Vy[1];
1936  rRotationMatrix(1,2) = Vy[2];
1937 
1938  rRotationMatrix(2,0) = Vz[0];
1939  rRotationMatrix(2,1) = Vz[1];
1940  rRotationMatrix(2,2) = Vz[2];
1941  }
1942 
1944 
1945  void AverageTipCoordinates(
1946  double& rTipX,
1947  double& rTipY,
1948  double& rTipZ,
1949  double& rTipDenominator,
1950  const FracturePoint& MyFracturePoint)
1951  {
1952  rTipX += MyFracturePoint.Weight * MyFracturePoint.Damage * (MyFracturePoint.Coordinates[0]);
1953  rTipY += MyFracturePoint.Weight * MyFracturePoint.Damage * (MyFracturePoint.Coordinates[1]);
1954  rTipZ += MyFracturePoint.Weight * MyFracturePoint.Damage * (MyFracturePoint.Coordinates[2]);
1955  rTipDenominator += MyFracturePoint.Weight * MyFracturePoint.Damage;
1956  }
1957 
1959 
1960  void CalculateNewTipRotationMatrix(
1961  PropagationLocalVariables& rAuxPropagationVariables,
1962  const array_1d<double,3>& NewTipCoordinates,
1963  const array_1d<double,3>& LeftInitCoordinates,
1964  const double& PropagationDistance)
1965  {
1966  // Unitary vector in local x direction
1967  array_1d<double,3> Vx;
1968  noalias(Vx) = NewTipCoordinates - rAuxPropagationVariables.TipCoordinates;
1969  double inv_norm = 1.0/norm_2(Vx);
1970  Vx[0] *= inv_norm;
1971  Vx[1] *= inv_norm;
1972  Vx[2] *= inv_norm;
1973 
1974  array_1d<double,3> CenterPoint;
1975  noalias(CenterPoint) = rAuxPropagationVariables.TipCoordinates + PropagationDistance*(Vx);
1976 
1977  array_1d<double,3> LeftLine;
1978  noalias(LeftLine) = NewTipCoordinates - LeftInitCoordinates;
1979 
1980  const double Lambda = (Vx[0]*(CenterPoint[0]-LeftInitCoordinates[0])+Vx[1]*(CenterPoint[1]-LeftInitCoordinates[1])+
1981  Vx[2]*(CenterPoint[2]-LeftInitCoordinates[2]))/(LeftLine[0]*Vx[0]+LeftLine[1]*Vx[1]+LeftLine[2]*Vx[2]);
1982 
1983  array_1d<double,3> LeftPoint;
1984  noalias(LeftPoint) = LeftInitCoordinates + Lambda*LeftLine;
1985 
1986  // Vector in local y direction
1987  array_1d<double,3> Vy;
1988  noalias(Vy) = LeftPoint - CenterPoint;
1989  // Unitary vector in local z direction
1990  array_1d<double,3> Vz;
1991  MathUtils<double>::CrossProduct(Vz,Vx,Vy);
1992  inv_norm = 1.0/norm_2(Vz);
1993  Vz[0] *= inv_norm;
1994  Vz[1] *= inv_norm;
1995  Vz[2] *= inv_norm;
1996 
1997  // Unitary vector in local y direction
1998  MathUtils<double>::CrossProduct(Vy, Vz, Vx);
1999 
2000  // Rotation Matrix
2001  rAuxPropagationVariables.RotationMatrix(0,0) = Vx[0];
2002  rAuxPropagationVariables.RotationMatrix(0,1) = Vx[1];
2003  rAuxPropagationVariables.RotationMatrix(0,2) = Vx[2];
2004 
2005  rAuxPropagationVariables.RotationMatrix(1,0) = Vy[0];
2006  rAuxPropagationVariables.RotationMatrix(1,1) = Vy[1];
2007  rAuxPropagationVariables.RotationMatrix(1,2) = Vy[2];
2008 
2009  rAuxPropagationVariables.RotationMatrix(2,0) = Vz[0];
2010  rAuxPropagationVariables.RotationMatrix(2,1) = Vz[1];
2011  rAuxPropagationVariables.RotationMatrix(2,2) = Vz[2];
2012 
2013  noalias(rAuxPropagationVariables.TipLocalCoordinates) = prod(rAuxPropagationVariables.RotationMatrix,CenterPoint);
2014  }
2015 
2016 }; // Class FracturePropagation3DUtilities
2017 
2018 } // namespace Kratos.
2019 
2020 #endif /* KRATOS_PROPAGATE_FRACTURES_3D_UTILITIES defined */
Kratos::Element::GeometryType
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: element.h:83
Kratos::FracturePropagation3DUtilities
Definition: fracture_propagation_3D_utilities.hpp:38
Kratos::FracturePropagation3DUtilities::MappingModelParts
void MappingModelParts(Parameters &rParameters, ModelPart &rModelPartOld, ModelPart &rModelPartNew, const bool &move_mesh_flag)
Definition: fracture_propagation_3D_utilities.hpp:162
Kratos::FracturePropagation3DUtilities::CheckFracture
void CheckFracture(const unsigned int &itFracture, PropagationGlobalVariables &rPropagationData, const UtilityVariables &AuxVariables, Parameters &rParameters)
Definition: fracture_propagation_3D_utilities.hpp:205
Kratos::FracturePropagation3DUtilities::FinalizeCheckFracture
void FinalizeCheckFracture(const PropagationGlobalVariables &PropagationData, Parameters &rParameters, ModelPart &rModelPart, const bool &move_mesh_flag)
Definition: fracture_propagation_3D_utilities.hpp:321
Kratos::FracturePropagation3DUtilities::WritePropagationData
void WritePropagationData(const PropagationGlobalVariables &PropagationData, Parameters &rParameters)
Definition: fracture_propagation_3D_utilities.hpp:342
Kratos::FracturePropagation3DUtilities::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(FracturePropagation3DUtilities)
Kratos::FracturePropagation3DUtilities::FracturePropagation3DUtilities
FracturePropagation3DUtilities()
Constructor.
Definition: fracture_propagation_3D_utilities.hpp:132
Kratos::FracturePropagation3DUtilities::CheckFracturePropagation
bool CheckFracturePropagation(Parameters &rParameters, ModelPart &rModelPart, const bool &move_mesh_flag)
Definition: fracture_propagation_3D_utilities.hpp:141
Kratos::FracturePropagation3DUtilities::ResetMeshPosition
void ResetMeshPosition(ModelPart &rModelPart)
Common ----------------------------------------------------------------------------------------------...
Definition: fracture_propagation_3D_utilities.hpp:1304
Kratos::FracturePropagation3DUtilities::UpdateMeshPosition
void UpdateMeshPosition(ModelPart &rModelPart)
Definition: fracture_propagation_3D_utilities.hpp:1323
Kratos::FracturePropagation3DUtilities::InitializeCheckFracture
void InitializeCheckFracture(PropagationGlobalVariables &rPropagationData, UtilityVariables &rAuxVariables, Parameters &rParameters, ModelPart &rModelPart, const bool &move_mesh_flag)
Member Variables.
Definition: fracture_propagation_3D_utilities.hpp:185
Kratos::FracturePropagation3DUtilities::GaussPointStateVariableMapping
void GaussPointStateVariableMapping(const UtilityVariables &AuxVariables, Parameters &rParameters, ModelPart &rModelPartOld, ModelPart &rModelPartNew)
Definition: fracture_propagation_3D_utilities.hpp:949
Kratos::FracturePropagation3DUtilities::~FracturePropagation3DUtilities
virtual ~FracturePropagation3DUtilities()
Destructor.
Definition: fracture_propagation_3D_utilities.hpp:137
Kratos::FracturePropagation3DUtilities::NodalVariablesMapping
void NodalVariablesMapping(const UtilityVariables &AuxVariables, Parameters &rParameters, ModelPart &rModelPartOld, ModelPart &rModelPartNew)
Definition: fracture_propagation_3D_utilities.hpp:632
Kratos::FracturePropagation3DUtilities::InitializeMapping
void InitializeMapping(UtilityVariables &rAuxVariables, ModelPart &rModelPartNew, const bool &move_mesh_flag)
Mapping ---------------------------------------------------------------------------------------------...
Definition: fracture_propagation_3D_utilities.hpp:618
Kratos::GeometryData::IntegrationMethod
IntegrationMethod
Definition: geometry_data.h:76
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::GlobalCoordinates
virtual CoordinatesArrayType & GlobalCoordinates(CoordinatesArrayType &rResult, CoordinatesArrayType const &LocalCoordinates) const
Definition: geometry.h:2377
Kratos::Geometry::IntegrationPointsArrayType
std::vector< IntegrationPointType > IntegrationPointsArrayType
Definition: geometry.h:161
Kratos::Geometry::DeterminantOfJacobian
Vector & DeterminantOfJacobian(Vector &rResult) const
Definition: geometry.h:3154
Kratos::Geometry::IntegrationPoints
const IntegrationPointsArrayType & IntegrationPoints() const
Definition: geometry.h:2284
Kratos::Internals::Matrix
Definition: amatrix_interface.h:41
Kratos::MathUtils::CrossProduct
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
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::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::Elements
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
Kratos::ModelPart::GetSubModelPart
ModelPart & GetSubModelPart(std::string const &SubModelPartName)
Definition: model_part.cpp:2029
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::Parameters
This class provides to Kratos a data structure for I/O based on the standard of JSON.
Definition: kratos_parameters.h:59
Kratos::Parameters::GetDouble
double GetDouble() const
This method returns the double contained in the current Parameter.
Definition: kratos_parameters.cpp:657
Kratos::Parameters::GetInt
int GetInt() const
This method returns the integer contained in the current Parameter.
Definition: kratos_parameters.cpp:666
Kratos::Parameters::size
SizeType size() const
This method returns the total size of the current array parameter.
Definition: kratos_parameters.cpp:993
Kratos::Parameters::GetString
std::string GetString() const
This method returns the string contained in the current Parameter.
Definition: kratos_parameters.cpp:684
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::array_1d< double, 3 >
kratos_parameters.h
model_part.h
CSharpKratosWrapper::ModelPart
Kratos::ModelPart ModelPart
Definition: kratos_wrapper.h:31
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::norm_2
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
Kratos::Vector
Internals::Matrix< double, AMatrix::dynamic, 1 > Vector
Definition: amatrix_interface.h:472
Kratos::prod
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
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::trans
AMatrix::TransposeMatrix< const T > trans(const T &TheMatrix)
Definition: amatrix_interface.h:486
Kratos::int
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
edgebased_PureConvection.ProcessInfo
ProcessInfo
Definition: edgebased_PureConvection.py:116
isotropic_damage_automatic_differentiation.out
out
Definition: isotropic_damage_automatic_differentiation.py:200
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
parallel_utilities.h
poromechanics_application_variables.h
Kratos::FracturePropagation3DUtilities::Bifurcation
Definition: fracture_propagation_3D_utilities.hpp:79
Kratos::FracturePropagation3DUtilities::Bifurcation::BotBotEndCoordinates
array_1d< double, 3 > BotBotEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:91
Kratos::FracturePropagation3DUtilities::Bifurcation::BotLeftEndCoordinates
array_1d< double, 3 > BotLeftEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:92
Kratos::FracturePropagation3DUtilities::Bifurcation::TopLeftEndCoordinates
array_1d< double, 3 > TopLeftEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:87
Kratos::FracturePropagation3DUtilities::Bifurcation::BotTipCoordinates
array_1d< double, 3 > BotTipCoordinates
Definition: fracture_propagation_3D_utilities.hpp:94
Kratos::FracturePropagation3DUtilities::Bifurcation::LeftInitCoordinates
array_1d< double, 3 > LeftInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:83
Kratos::FracturePropagation3DUtilities::Bifurcation::MotherFractureId
int MotherFractureId
Definition: fracture_propagation_3D_utilities.hpp:80
Kratos::FracturePropagation3DUtilities::Bifurcation::TopTopEndCoordinates
array_1d< double, 3 > TopTopEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:85
Kratos::FracturePropagation3DUtilities::Bifurcation::TopBotEndCoordinates
array_1d< double, 3 > TopBotEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:86
Kratos::FracturePropagation3DUtilities::Bifurcation::TopTipCoordinates
array_1d< double, 3 > TopTipCoordinates
Definition: fracture_propagation_3D_utilities.hpp:89
Kratos::FracturePropagation3DUtilities::Bifurcation::BotInitCoordinates
array_1d< double, 3 > BotInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:82
Kratos::FracturePropagation3DUtilities::Bifurcation::TopInitCoordinates
array_1d< double, 3 > TopInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:81
Kratos::FracturePropagation3DUtilities::Bifurcation::TopRightEndCoordinates
array_1d< double, 3 > TopRightEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:88
Kratos::FracturePropagation3DUtilities::Bifurcation::BotRightEndCoordinates
array_1d< double, 3 > BotRightEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:93
Kratos::FracturePropagation3DUtilities::Bifurcation::RightInitCoordinates
array_1d< double, 3 > RightInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:84
Kratos::FracturePropagation3DUtilities::Bifurcation::BotTopEndCoordinates
array_1d< double, 3 > BotTopEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:90
Kratos::FracturePropagation3DUtilities::FracturePoint
Structs for fracture propagation check --------------------------------------------------------------...
Definition: fracture_propagation_3D_utilities.hpp:54
Kratos::FracturePropagation3DUtilities::FracturePoint::TipDistance
double TipDistance
Definition: fracture_propagation_3D_utilities.hpp:56
Kratos::FracturePropagation3DUtilities::FracturePoint::Weight
double Weight
Definition: fracture_propagation_3D_utilities.hpp:56
Kratos::FracturePropagation3DUtilities::FracturePoint::Damage
double Damage
Definition: fracture_propagation_3D_utilities.hpp:56
Kratos::FracturePropagation3DUtilities::FracturePoint::Coordinates
array_1d< double, 3 > Coordinates
Definition: fracture_propagation_3D_utilities.hpp:55
Kratos::FracturePropagation3DUtilities::FracturePoint::pElement
Element::Pointer pElement
Definition: fracture_propagation_3D_utilities.hpp:57
Kratos::FracturePropagation3DUtilities::GaussPointOld
Structs for mapping model parts ---------------------------------------------------------------------...
Definition: fracture_propagation_3D_utilities.hpp:122
Kratos::FracturePropagation3DUtilities::GaussPointOld::Coordinates
array_1d< double, 3 > Coordinates
Definition: fracture_propagation_3D_utilities.hpp:123
Kratos::FracturePropagation3DUtilities::GaussPointOld::Weight
double Weight
Definition: fracture_propagation_3D_utilities.hpp:124
Kratos::FracturePropagation3DUtilities::GaussPointOld::StateVariable
double StateVariable
Definition: fracture_propagation_3D_utilities.hpp:124
Kratos::FracturePropagation3DUtilities::PropagationGlobalVariables
Definition: fracture_propagation_3D_utilities.hpp:100
Kratos::FracturePropagation3DUtilities::PropagationGlobalVariables::PropagateFractures
bool PropagateFractures
Definition: fracture_propagation_3D_utilities.hpp:103
Kratos::FracturePropagation3DUtilities::PropagationGlobalVariables::PropagationVector
std::vector< Propagation > PropagationVector
Definition: fracture_propagation_3D_utilities.hpp:104
Kratos::FracturePropagation3DUtilities::PropagationGlobalVariables::BifurcationVector
std::vector< Bifurcation > BifurcationVector
Definition: fracture_propagation_3D_utilities.hpp:105
Kratos::FracturePropagation3DUtilities::PropagationGlobalVariables::FracturePointsCellMatrix
std::vector< std::vector< std::vector< std::vector< FracturePoint > > > > FracturePointsCellMatrix
Definition: fracture_propagation_3D_utilities.hpp:101
Kratos::FracturePropagation3DUtilities::PropagationGlobalVariables::pProcessInfo
ProcessInfo::Pointer pProcessInfo
Definition: fracture_propagation_3D_utilities.hpp:102
Kratos::FracturePropagation3DUtilities::Propagation
Definition: fracture_propagation_3D_utilities.hpp:63
Kratos::FracturePropagation3DUtilities::Propagation::RightInitCoordinates
array_1d< double, 3 > RightInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:68
Kratos::FracturePropagation3DUtilities::Propagation::BotInitCoordinates
array_1d< double, 3 > BotInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:66
Kratos::FracturePropagation3DUtilities::Propagation::TopInitCoordinates
array_1d< double, 3 > TopInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:65
Kratos::FracturePropagation3DUtilities::Propagation::BotEndCoordinates
array_1d< double, 3 > BotEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:70
Kratos::FracturePropagation3DUtilities::Propagation::RightEndCoordinates
array_1d< double, 3 > RightEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:72
Kratos::FracturePropagation3DUtilities::Propagation::TipCoordinates
array_1d< double, 3 > TipCoordinates
Definition: fracture_propagation_3D_utilities.hpp:73
Kratos::FracturePropagation3DUtilities::Propagation::TopEndCoordinates
array_1d< double, 3 > TopEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:69
Kratos::FracturePropagation3DUtilities::Propagation::LeftInitCoordinates
array_1d< double, 3 > LeftInitCoordinates
Definition: fracture_propagation_3D_utilities.hpp:67
Kratos::FracturePropagation3DUtilities::Propagation::MotherFractureId
int MotherFractureId
Definition: fracture_propagation_3D_utilities.hpp:64
Kratos::FracturePropagation3DUtilities::Propagation::LeftEndCoordinates
array_1d< double, 3 > LeftEndCoordinates
Definition: fracture_propagation_3D_utilities.hpp:71
Kratos::FracturePropagation3DUtilities::PropagationLocalVariables
Definition: fracture_propagation_3D_utilities.hpp:111
Kratos::FracturePropagation3DUtilities::PropagationLocalVariables::BotFrontFracturePoints
std::vector< FracturePoint * > BotFrontFracturePoints
Definition: fracture_propagation_3D_utilities.hpp:116
Kratos::FracturePropagation3DUtilities::PropagationLocalVariables::TipCoordinates
array_1d< double, 3 > TipCoordinates
Definition: fracture_propagation_3D_utilities.hpp:113
Kratos::FracturePropagation3DUtilities::PropagationLocalVariables::TipLocalCoordinates
array_1d< double, 3 > TipLocalCoordinates
Definition: fracture_propagation_3D_utilities.hpp:114
Kratos::FracturePropagation3DUtilities::PropagationLocalVariables::RotationMatrix
BoundedMatrix< double, 3, 3 > RotationMatrix
Definition: fracture_propagation_3D_utilities.hpp:112
Kratos::FracturePropagation3DUtilities::PropagationLocalVariables::TopFrontFracturePoints
std::vector< FracturePoint * > TopFrontFracturePoints
Definition: fracture_propagation_3D_utilities.hpp:115
Kratos::FracturePropagation3DUtilities::UtilityVariables
Basic Structs for the utility -----------------------------------------------------------------------...
Definition: fracture_propagation_3D_utilities.hpp:45
Kratos::FracturePropagation3DUtilities::UtilityVariables::Y_min
double Y_min
Definition: fracture_propagation_3D_utilities.hpp:46
Kratos::FracturePropagation3DUtilities::UtilityVariables::RowSize
double RowSize
Definition: fracture_propagation_3D_utilities.hpp:48
Kratos::FracturePropagation3DUtilities::UtilityVariables::X_min
double X_min
Definition: fracture_propagation_3D_utilities.hpp:46
Kratos::FracturePropagation3DUtilities::UtilityVariables::NColumns
int NColumns
Definition: fracture_propagation_3D_utilities.hpp:47
Kratos::FracturePropagation3DUtilities::UtilityVariables::Z_max
double Z_max
Definition: fracture_propagation_3D_utilities.hpp:46
Kratos::FracturePropagation3DUtilities::UtilityVariables::ColumnSize
double ColumnSize
Definition: fracture_propagation_3D_utilities.hpp:48
Kratos::FracturePropagation3DUtilities::UtilityVariables::SectionSize
double SectionSize
Definition: fracture_propagation_3D_utilities.hpp:48
Kratos::FracturePropagation3DUtilities::UtilityVariables::NRows
int NRows
Definition: fracture_propagation_3D_utilities.hpp:47
Kratos::FracturePropagation3DUtilities::UtilityVariables::Z_min
double Z_min
Definition: fracture_propagation_3D_utilities.hpp:46
Kratos::FracturePropagation3DUtilities::UtilityVariables::NSections
int NSections
Definition: fracture_propagation_3D_utilities.hpp:47
Kratos::FracturePropagation3DUtilities::UtilityVariables::Y_max
double Y_max
Definition: fracture_propagation_3D_utilities.hpp:46
Kratos::FracturePropagation3DUtilities::UtilityVariables::X_max
double X_max
Definition: fracture_propagation_3D_utilities.hpp:46
Generated on Sat Feb 3 2024 02:11:11 for KratosMultiphysics by doxygen 1.9.1