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_2D_utilities.hpp
Go to the documentation of this file.
1 
2 // | / |
3 // ' / __| _` | __| _ \ __|
4 // . \ | ( | | ( |\__ `
5 // _|\_\_| \__,_|\__|\___/ ____/
6 // Multi-Physics
7 //
8 // License: BSD License
9 // Kratos default license: kratos/license.txt
10 //
11 // Main authors: Ignasi de Pouplana
12 //
13 
14 
15 #if !defined(KRATOS_PROPAGATE_FRACTURES_2D_UTILITIES )
16 #define KRATOS_PROPAGATE_FRACTURES_2D_UTILITIES
17 
18 // System includes
19 #include <fstream>
20 #include <iostream>
21 #include <cmath>
22 
23 // Project includes
24 #include "geometries/geometry.h"
25 #include "includes/define.h"
26 #include "includes/model_part.h"
27 #include "includes/kratos_parameters.h"
28 #include "utilities/openmp_utils.h"
29 #include "utilities/parallel_utilities.h"
30 
31 // Application includes
32 #include "poromechanics_application_variables.h"
33 
34 namespace Kratos
35 {
36 
37 class FracturePropagation2DUtilities
38 {
39 
40 protected:
41 
43 
44  struct UtilityVariables
45  {
46  double X_max, X_min, Y_max, Y_min;
47  int NRows, NColumns;
48  double RowSize, ColumnSize;
49  };
50 
52 
53  struct FracturePoint
54  {
55  array_1d<double,2> 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> TipCoordinates;
68  };
69 
71 
72  struct Bifurcation
73  {
74  int MotherFractureId;
75  array_1d<double,3> TopInitCoordinates;
76  array_1d<double,3> BotInitCoordinates;
77  array_1d<double,3> TopTipCoordinates;
78  array_1d<double,3> BotTipCoordinates;
79  };
80 
82 
83  struct PropagationGlobalVariables
84  {
85  std::vector< std::vector< std::vector<FracturePoint> > > FracturePointsCellMatrix;
86  ProcessInfo::Pointer pProcessInfo;
87  bool PropagateFractures;
88  std::vector<Propagation> PropagationVector;
89  std::vector<Bifurcation> BifurcationVector;
90  };
91 
93 
94  struct PropagationLocalVariables
95  {
96  BoundedMatrix<double,2,2> RotationMatrix;
97  array_1d<double,2> TipCoordinates;
98  array_1d<double,2> TipLocalCoordinates;
99  std::vector<FracturePoint*> TopFrontFracturePoints;
100  std::vector<FracturePoint*> BotFrontFracturePoints;
101  };
102 
104 
105  struct GaussPointOld
106  {
107  array_1d<double,3> Coordinates;
108  double StateVariable, Weight;
109  };
110 
111 public:
112 
113  KRATOS_CLASS_POINTER_DEFINITION( FracturePropagation2DUtilities );
114 
116  FracturePropagation2DUtilities() {}
117 
119 
121  virtual ~FracturePropagation2DUtilities() {}
122 
124 
125  bool CheckFracturePropagation (Parameters& rParameters, ModelPart& rModelPart, const bool& move_mesh_flag)
126  {
127  // Define necessary variables
128  UtilityVariables AuxVariables;
129  PropagationGlobalVariables PropagationData;
130 
131  this->InitializeCheckFracture(PropagationData, AuxVariables, rParameters, rModelPart, move_mesh_flag);
132 
133  //Loop for all pre-existing fractures
134  for(unsigned int i = 0; i < rParameters["fractures_list"].size(); i++)
135  {
136  this->CheckFracture(i, PropagationData, AuxVariables, rParameters);
137  }
138 
139  this->FinalizeCheckFracture(PropagationData, rParameters, rModelPart, move_mesh_flag);
140 
141  return PropagationData.PropagateFractures;
142  }
143 
145 
146  void MappingModelParts (Parameters& rParameters, ModelPart& rModelPartOld, ModelPart& rModelPartNew, const bool& move_mesh_flag)
147  {
148  // Define necessary variables
149  UtilityVariables AuxVariables;
150 
151  this->InitializeMapping(AuxVariables,rModelPartNew, move_mesh_flag);
152 
153  this->NodalVariablesMapping(AuxVariables,rParameters,rModelPartOld,rModelPartNew);
154 
155  this->GaussPointStateVariableMapping(AuxVariables,rParameters,rModelPartOld,rModelPartNew);
156 
157  if(move_mesh_flag==true)
158  this->UpdateMeshPosition(rModelPartNew);
159  }
160 
161 protected:
162 
164 
166 
168 
169  void InitializeCheckFracture(
170  PropagationGlobalVariables& rPropagationData,
171  UtilityVariables& rAuxVariables,
172  Parameters& rParameters,
173  ModelPart& rModelPart,
174  const bool& move_mesh_flag)
175  {
176  // Set PropagateFractures bool to false
177  rPropagationData.PropagateFractures = false;
178 
179  // Move mesh to the original position to work in the reference state
180  if(move_mesh_flag==true)
181  this->ResetMeshPosition(rModelPart);
182 
183  // Locate fracture points inside CellMatrix
184  this->SetFracturePoints(rPropagationData, rAuxVariables, rParameters, rModelPart);
185  }
186 
188 
189  void CheckFracture(
190  const unsigned int& itFracture,
191  PropagationGlobalVariables& rPropagationData,
192  const UtilityVariables& AuxVariables,
193  Parameters& rParameters)
194  {
195  PropagationLocalVariables AuxPropagationVariables;
196 
197  // Tip Coordinates
198  for(unsigned int i = 0; i < 2; i++)
199  AuxPropagationVariables.TipCoordinates[i] = rParameters["fractures_list"][itFracture]["tip_point"]["coordinates"][i].GetDouble();
200 
201  // Tip Rotation Matrix
202  this->CalculateTipRotationMatrix(itFracture,AuxPropagationVariables.RotationMatrix,rParameters);
203 
204  // Tip search area
205  const double PropagationLength = rParameters["fracture_data"]["propagation_length"].GetDouble();
206 
207  double X_left = AuxPropagationVariables.TipCoordinates[0] - PropagationLength;
208  double X_right = AuxPropagationVariables.TipCoordinates[0] + PropagationLength;
209  double Y_top = AuxPropagationVariables.TipCoordinates[1] + PropagationLength;
210  double Y_bot = AuxPropagationVariables.TipCoordinates[1] - PropagationLength;
211 
212  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
213  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
214  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
215  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
216 
217  if(Column_left < 0) Column_left = 0;
218  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
219  if(Row_top < 0) Row_top = 0;
220  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
221 
222  // Search FracturePoints neighbours around the tip
223  std::vector<FracturePoint*> TipNeighbours;
224  noalias(AuxPropagationVariables.TipLocalCoordinates) = prod(AuxPropagationVariables.RotationMatrix,AuxPropagationVariables.TipCoordinates);
225  array_1d<double,2> OtherLocalCoordinates;
226  for(int i = Row_top; i <= Row_bot; i++)
227  {
228  for(int j = Column_left; j<= Column_right; j++)
229  {
230  for(unsigned int k = 0; k < rPropagationData.FracturePointsCellMatrix[i][j].size(); k++)
231  {
232  FracturePoint& rOtherPoint = rPropagationData.FracturePointsCellMatrix[i][j][k];
233 
234  rOtherPoint.TipDistance = sqrt((rOtherPoint.Coordinates[0]-AuxPropagationVariables.TipCoordinates[0])*(rOtherPoint.Coordinates[0]-AuxPropagationVariables.TipCoordinates[0]) +
235  (rOtherPoint.Coordinates[1]-AuxPropagationVariables.TipCoordinates[1])*(rOtherPoint.Coordinates[1]-AuxPropagationVariables.TipCoordinates[1]));
236 
237  if(rOtherPoint.TipDistance <= PropagationLength)
238  {
239  TipNeighbours.push_back(&rOtherPoint);
240 
241  noalias(OtherLocalCoordinates) = prod(AuxPropagationVariables.RotationMatrix,rOtherPoint.Coordinates);
242 
243  // FrontFracturePoints
244  if(OtherLocalCoordinates[0] >= AuxPropagationVariables.TipLocalCoordinates[0])
245  {
246  // TopFrontFracturePoints
247  if(OtherLocalCoordinates[1] >= AuxPropagationVariables.TipLocalCoordinates[1])
248  {
249  AuxPropagationVariables.TopFrontFracturePoints.push_back(&rOtherPoint);
250  }
251  // BotFrontFracturePoints
252  else
253  {
254  AuxPropagationVariables.BotFrontFracturePoints.push_back(&rOtherPoint);
255  }
256  }
257  }
258  }
259  }
260  }
261 
262  // Calculate Non-local damage around the tip
263  double NonlocalDamage = 0.0;
264  double WeightingFunctionDenominator = 0.0;
265  int NNeighbours = TipNeighbours.size();
266 
267  #pragma omp parallel for reduction(+:NonlocalDamage,WeightingFunctionDenominator)
268  for(int i = 0; i < NNeighbours; i++)
269  {
270  const FracturePoint& MyPoint = *(TipNeighbours[i]);
271  NonlocalDamage += (MyPoint.Weight)*
272  exp(-4.0*(MyPoint.TipDistance)*(MyPoint.TipDistance)/(PropagationLength*PropagationLength))*
273  (MyPoint.Damage);
274  WeightingFunctionDenominator += (MyPoint.Weight)*exp(-4.0*(MyPoint.TipDistance)*(MyPoint.TipDistance)/(PropagationLength*PropagationLength));
275  }
276 
277  if(WeightingFunctionDenominator > 1.0e-20)
278  NonlocalDamage = NonlocalDamage/WeightingFunctionDenominator;
279  else
280  NonlocalDamage = 0.0;
281 
282  // Check fracture propagation
283  if(NonlocalDamage >= rParameters["fracture_data"]["propagation_damage"].GetDouble())
284  this->PropagateFracture(itFracture,rPropagationData,AuxPropagationVariables,rParameters);
285  }
286 
288 
289  void FinalizeCheckFracture(
290  const PropagationGlobalVariables& PropagationData,
291  Parameters& rParameters,
292  ModelPart& rModelPart,
293  const bool& move_mesh_flag)
294  {
295  if(PropagationData.PropagateFractures==false)
296  {
297  // Move mesh to the current position
298  if(move_mesh_flag==true)
299  this->UpdateMeshPosition(rModelPart);
300  }
301  else
302  {
303  // Write PropagationData.tcl file
304  this->WritePropagationData(PropagationData,rParameters);
305  }
306  }
307 
309 
310  void WritePropagationData(
311  const PropagationGlobalVariables& PropagationData,
312  Parameters& rParameters)
313  {
314 
315  std::fstream PropDataFile;
316  PropDataFile.open ("PropagationData.tcl", std::fstream::out | std::fstream::trunc);
317  PropDataFile.precision(12);
318 
319  PropDataFile << "set PropagationData [list]" << std::endl;
320 
321  // GiDPath
322  PropDataFile << "lappend PropagationData \"" << rParameters["fracture_data"]["gid_path"].GetString() << "\"" << std::endl;
323 
324  int Id;
325 
326  // BodySurfacesDict
327  for(unsigned int i = 0; i < rParameters["body_surfaces_list"].size(); i++)
328  {
329  Id = rParameters["body_surfaces_list"][i]["id"].GetInt();
330 
331  PropDataFile << "set Groups [list]" << std::endl;
332  for(unsigned int j = 0; j < rParameters["body_surfaces_list"][i]["groups"].size(); j++)
333  {
334  PropDataFile << "lappend Groups \"" << rParameters["body_surfaces_list"][i]["groups"][j].GetString() << "\"" << std::endl;
335  }
336  PropDataFile << "dict set BodySurfacesDict " << Id << " Groups $Groups" << std::endl;
337 
338  PropDataFile << "set Lines [list]" << std::endl;
339  for(unsigned int j = 0; j < rParameters["body_surfaces_list"][i]["lines"].size(); j++)
340  {
341  PropDataFile << "lappend Lines " << rParameters["body_surfaces_list"][i]["lines"][j].GetInt() << std::endl;
342  }
343  PropDataFile << "dict set BodySurfacesDict " << Id << " Lines $Lines" << std::endl;
344  PropDataFile << "dict set BodySurfacesDict " << Id << " ElemType " << rParameters["body_surfaces_list"][i]["elem_type"].GetString() << std::endl;
345  PropDataFile << "dict set BodySurfacesDict " << Id << " MeshSize " << rParameters["body_surfaces_list"][i]["mesh_size"].GetDouble() << std::endl;
346  }
347  PropDataFile << "lappend PropagationData $BodySurfacesDict" << std::endl;
348 
349  // FracturesDict
350  for(unsigned int i = 0; i < rParameters["fractures_list"].size(); i++)
351  {
352  Id = rParameters["fractures_list"][i]["id"].GetInt();
353 
354  // TipPoint
355  PropDataFile << "dict set FracturesDict " << Id << " TipPoint Id "
356  << rParameters["fractures_list"][i]["tip_point"]["id"].GetInt() << std::endl;
357  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["tip_point"]["coordinates"][0].GetDouble()
358  << " " << rParameters["fractures_list"][i]["tip_point"]["coordinates"][1].GetDouble() << " "
359  << rParameters["fractures_list"][i]["tip_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
360  PropDataFile << "dict set FracturesDict " << Id << " TipPoint Coordinates $Coordinates" << std::endl;
361  // TopPoint
362  PropDataFile << "dict set FracturesDict " << Id << " TopPoint Id "
363  << rParameters["fractures_list"][i]["top_point"]["id"].GetInt() << std::endl;
364  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["top_point"]["coordinates"][0].GetDouble()
365  << " " << rParameters["fractures_list"][i]["top_point"]["coordinates"][1].GetDouble() << " "
366  << rParameters["fractures_list"][i]["top_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
367  PropDataFile << "dict set FracturesDict " << Id << " TopPoint Coordinates $Coordinates" << std::endl;
368  // BotPoint
369  PropDataFile << "dict set FracturesDict " << Id << " BotPoint Id "
370  << rParameters["fractures_list"][i]["bot_point"]["id"].GetInt() << std::endl;
371  PropDataFile << "set Coordinates \"" << rParameters["fractures_list"][i]["bot_point"]["coordinates"][0].GetDouble()
372  << " " << rParameters["fractures_list"][i]["bot_point"]["coordinates"][1].GetDouble() << " "
373  << rParameters["fractures_list"][i]["bot_point"]["coordinates"][2].GetDouble() << "\"" << std::endl;
374  PropDataFile << "dict set FracturesDict " << Id << " BotPoint Coordinates $Coordinates" << std::endl;
375  // TopLine
376  PropDataFile << "dict set FracturesDict " << Id << " TopLine Id "
377  << rParameters["fractures_list"][i]["top_line"]["id"].GetInt() << std::endl;
378  // BotLine
379  PropDataFile << "dict set FracturesDict " << Id << " BotLine Id "
380  << rParameters["fractures_list"][i]["bot_line"]["id"].GetInt() << std::endl;
381  // InterfaceSurface
382  PropDataFile << "dict set FracturesDict " << Id << " InterfaceSurface Id "
383  << rParameters["fractures_list"][i]["interface_surface"]["id"].GetInt() << std::endl;
384  PropDataFile << "dict set FracturesDict " << Id << " InterfaceSurface Layer \""
385  << rParameters["fractures_list"][i]["interface_surface"]["layer"].GetString() << "\"" << std::endl;
386  PropDataFile << "set Groups [list]" << std::endl;
387  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["interface_surface"]["groups"].size(); j++)
388  {
389  PropDataFile << "lappend Groups \"" << rParameters["fractures_list"][i]["interface_surface"]["groups"][j].GetString() << "\"" << std::endl;
390  }
391  PropDataFile << "dict set FracturesDict " << Id << " InterfaceSurface Groups $Groups" << std::endl;
392  // BodySurfaces
393  PropDataFile << "set BodySurfaces [list]" << std::endl;
394  for(unsigned int j = 0; j < rParameters["fractures_list"][i]["body_surfaces"].size(); j++)
395  {
396  PropDataFile << "lappend BodySurfaces " << rParameters["fractures_list"][i]["body_surfaces"][j].GetInt() << std::endl;
397  }
398  PropDataFile << "dict set FracturesDict " << Id << " BodySurfaces $BodySurfaces" << std::endl;
399  }
400  PropDataFile << "lappend PropagationData $FracturesDict" << std::endl;
401 
402  // PropagationDict
403  PropDataFile << "set PropagationDict [dict create]" << std::endl;
404  for(unsigned int i = 0; i < PropagationData.PropagationVector.size(); i++)
405  {
406  PropDataFile << "dict set PropagationDict " << i << " MotherFractureId "
407  << PropagationData.PropagationVector[i].MotherFractureId << std::endl;
408 
409  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].TopInitCoordinates[0]
410  << " " << PropagationData.PropagationVector[i].TopInitCoordinates[1] << " "
411  << PropagationData.PropagationVector[i].TopInitCoordinates[2] << "\"" << std::endl;
412  PropDataFile << "dict set PropagationDict " << i << " TopInitCoordinates $Coordinates" << std::endl;
413 
414  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].BotInitCoordinates[0]
415  << " " << PropagationData.PropagationVector[i].BotInitCoordinates[1] << " "
416  << PropagationData.PropagationVector[i].BotInitCoordinates[2] << "\"" << std::endl;
417  PropDataFile << "dict set PropagationDict " << i << " BotInitCoordinates $Coordinates" << std::endl;
418 
419  PropDataFile << "set Coordinates \"" << PropagationData.PropagationVector[i].TipCoordinates[0]
420  << " " << PropagationData.PropagationVector[i].TipCoordinates[1] << " "
421  << PropagationData.PropagationVector[i].TipCoordinates[2] << "\"" << std::endl;
422  PropDataFile << "dict set PropagationDict " << i << " TipCoordinates $Coordinates" << std::endl;
423  }
424  PropDataFile << "lappend PropagationData $PropagationDict" << std::endl;
425 
426  // BifurcationDict
427  PropDataFile << "set BifurcationDict [dict create]" << std::endl;
428  for(unsigned int i = 0; i < PropagationData.BifurcationVector.size(); i++)
429  {
430  PropDataFile << "dict set BifurcationDict " << i << " MotherFractureId "
431  << PropagationData.BifurcationVector[i].MotherFractureId << std::endl;
432 
433  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopInitCoordinates[0]
434  << " " << PropagationData.BifurcationVector[i].TopInitCoordinates[1] << " "
435  << PropagationData.BifurcationVector[i].TopInitCoordinates[2] << "\"" << std::endl;
436  PropDataFile << "dict set BifurcationDict " << i << " TopInitCoordinates $Coordinates" << std::endl;
437 
438  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotInitCoordinates[0]
439  << " " << PropagationData.BifurcationVector[i].BotInitCoordinates[1] << " "
440  << PropagationData.BifurcationVector[i].BotInitCoordinates[2] << "\"" << std::endl;
441  PropDataFile << "dict set BifurcationDict " << i << " BotInitCoordinates $Coordinates" << std::endl;
442 
443  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].TopTipCoordinates[0]
444  << " " << PropagationData.BifurcationVector[i].TopTipCoordinates[1] << " "
445  << PropagationData.BifurcationVector[i].TopTipCoordinates[2] << "\"" << std::endl;
446  PropDataFile << "dict set BifurcationDict " << i << " TopTipCoordinates $Coordinates" << std::endl;
447 
448  PropDataFile << "set Coordinates \"" << PropagationData.BifurcationVector[i].BotTipCoordinates[0]
449  << " " << PropagationData.BifurcationVector[i].BotTipCoordinates[1] << " "
450  << PropagationData.BifurcationVector[i].BotTipCoordinates[2] << "\"" << std::endl;
451  PropDataFile << "dict set BifurcationDict " << i << " BotTipCoordinates $Coordinates" << std::endl;
452  }
453  PropDataFile << "lappend PropagationData $BifurcationDict" << std::endl;
454 
455  PropDataFile << "return $PropagationData" << std::endl;
456 
457  PropDataFile.close();
458  }
459 
461 
462  void InitializeMapping(
463  UtilityVariables& rAuxVariables,
464  ModelPart& rModelPartNew,
465  const bool& move_mesh_flag)
466  {
467  // Move mesh to the original position to work in the reference state
468  if(move_mesh_flag==true)
469  this->ResetMeshPosition(rModelPartNew);
470 
471  this->ComputeCellMatrixDimensions(rAuxVariables,rModelPartNew);
472  }
473 
475 
476  void NodalVariablesMapping(
477  const UtilityVariables& AuxVariables,
478  Parameters& rParameters,
479  ModelPart& rModelPartOld,
480  ModelPart& rModelPartNew)
481  {
482  // Define ElementOld Cell matrix
483  std::vector< std::vector< std::vector<Element::Pointer> > > ElementOldCellMatrix;
484  ElementOldCellMatrix.resize(AuxVariables.NRows);
485  for(int i = 0; i < AuxVariables.NRows; i++) ElementOldCellMatrix[i].resize(AuxVariables.NColumns);
486 
487  // Locate Old Elments in cells
488  double X_me;
489  double Y_me;
490  int PointsNumber;
491 
492  unsigned int NumBodySubModelParts = rParameters["fracture_data"]["body_domain_sub_model_part_list"].size();
493 
494  // Loop through all BodySubModelParts
495  for(unsigned int m = 0; m < NumBodySubModelParts; m++)
496  {
497  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][m].GetString());
498 
499  int NElems = static_cast<int>(SubModelPart.Elements().size());
500  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
501 
502  #pragma omp parallel for private(X_me,Y_me,PointsNumber)
503  for(int i = 0; i < NElems; i++)
504  {
505  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
506 
507  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
508  double X_right = X_left;
509  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
510  double Y_bot = Y_top;
511  PointsNumber = itElemOld->GetGeometry().PointsNumber();
512 
513  for(int j = 1; j < PointsNumber; j++)
514  {
515  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
516  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
517 
518  if(X_me > X_right) X_right = X_me;
519  else if(X_me < X_left) X_left = X_me;
520  if(Y_me > Y_top) Y_top = Y_me;
521  else if(Y_me < Y_bot) Y_bot = Y_me;
522  }
523 
524  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
525  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
526  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
527  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
528 
529  if(Column_left < 0) Column_left = 0;
530  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
531  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
532  else if(Column_right < 0) Column_right = 0;
533 
534  if(Row_top < 0) Row_top = 0;
535  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
536  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
537  else if(Row_bot < 0) Row_bot = 0;
538 
539  for(int k = Row_top; k <= Row_bot; k++)
540  {
541  for(int l = Column_left; l <= Column_right; l++)
542  {
543  #pragma omp critical
544  {
545  ElementOldCellMatrix[k][l].push_back((*(itElemOld.base())));
546  }
547  }
548  }
549  }
550  }
551 
552  unsigned int NumInterfaceSubModelPartsOld = rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"].size();
553 
554  // Loop through all InterfaceSubModelParts
555  for(unsigned int m = 0; m < NumInterfaceSubModelPartsOld; m++)
556  {
557  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"][m].GetString());
558 
559  int NElems = static_cast<int>(SubModelPart.Elements().size());
560  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
561 
562  #pragma omp parallel for private(X_me,Y_me,PointsNumber)
563  for(int i = 0; i < NElems; i++)
564  {
565  ModelPart::ElementsContainerType::iterator itElemOld = el_begin + i;
566 
567  double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
568  double X_right = X_left;
569  double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
570  double Y_bot = Y_top;
571  PointsNumber = itElemOld->GetGeometry().PointsNumber();
572 
573  for(int j = 1; j < PointsNumber; j++)
574  {
575  X_me = itElemOld->GetGeometry().GetPoint(j).X0();
576  Y_me = itElemOld->GetGeometry().GetPoint(j).Y0();
577 
578  if(X_me > X_right) X_right = X_me;
579  else if(X_me < X_left) X_left = X_me;
580  if(Y_me > Y_top) Y_top = Y_me;
581  else if(Y_me < Y_bot) Y_bot = Y_me;
582  }
583 
584  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
585  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
586  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
587  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
588 
589  if(Column_left < 0) Column_left = 0;
590  else if(Column_left >= AuxVariables.NColumns) Column_left = AuxVariables.NColumns-1;
591  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
592  else if(Column_right < 0) Column_right = 0;
593 
594  if(Row_top < 0) Row_top = 0;
595  else if(Row_top >= AuxVariables.NRows) Row_top = AuxVariables.NRows-1;
596  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
597  else if(Row_bot < 0) Row_bot = 0;
598 
599  for(int k = Row_top; k <= Row_bot; k++)
600  {
601  for(int l = Column_left; l <= Column_right; l++)
602  {
603  #pragma omp critical
604  {
605  ElementOldCellMatrix[k][l].push_back((*(itElemOld.base())));
606  }
607  }
608  }
609  }
610  }
611 
612  // Locate new nodes inside old elements and interpolate nodal variables
613  const int NNodes = static_cast<int>(rModelPartNew.Nodes().size());
614  ModelPart::NodesContainerType::iterator node_begin = rModelPartNew.NodesBegin();
615 
616  array_1d<double,3> GlobalCoordinates;
617  array_1d<double,3> LocalCoordinates;
618 
619  #pragma omp parallel for private(X_me,Y_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
620  for(int i = 0; i < NNodes; i++)
621  {
622  ModelPart::NodesContainerType::iterator itNodeNew = node_begin + i;
623 
624  X_me = itNodeNew->X0();
625  Y_me = itNodeNew->Y0();
626 
627  int Column = int((X_me-AuxVariables.X_min)/AuxVariables.ColumnSize);
628  int Row = int((AuxVariables.Y_max-Y_me)/AuxVariables.RowSize);
629 
630  if(Column >= AuxVariables.NColumns) Column = AuxVariables.NColumns-1;
631  else if(Column < 0) Column = 0;
632  if(Row >= AuxVariables.NRows) Row = AuxVariables.NRows-1;
633  else if(Row < 0) Row = 0;
634 
635  noalias(GlobalCoordinates) = itNodeNew->Coordinates(); //Coordinates of new nodes are still in the original position
636  bool IsInside = false;
637  Element::Pointer pElementOld;
638 
639  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
640  {
641  pElementOld = ElementOldCellMatrix[Row][Column][m];
642  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates); //Checks whether the global coordinates fall inside the original old element
643  if(IsInside) break;
644  }
645  if(IsInside==false)
646  {
647  for(unsigned int m = 0; m < (ElementOldCellMatrix[Row][Column]).size(); m++)
648  {
649  pElementOld = ElementOldCellMatrix[Row][Column][m];
650  IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,1.0e-5);
651  if(IsInside) break;
652  }
653  }
654  if(IsInside == false)
655  std::cout << "ERROR!!, NONE OF THE OLD ELEMENTS CONTAINS NODE: " << itNodeNew->Id() << std::endl;
656 
657  PointsNumber = pElementOld->GetGeometry().PointsNumber();
658  Vector ShapeFunctionsValuesVector(PointsNumber);
659  Vector NodalVariableVector(PointsNumber);
660 
661  pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
662 
663  // Interpolation of nodal variables
664  if( itNodeNew->IsFixed(DISPLACEMENT_X)==false )
665  {
666  for(int j = 0; j < PointsNumber; j++)
667  {
668  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_X);
669  }
670  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
671  }
672  if( itNodeNew->IsFixed(VELOCITY_X)==false )
673  {
674  for(int j = 0; j < PointsNumber; j++)
675  {
676  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(VELOCITY_X);
677  }
678  itNodeNew->FastGetSolutionStepValue(VELOCITY_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
679  }
680  if( itNodeNew->IsFixed(ACCELERATION_X)==false )
681  {
682  for(int j = 0; j < PointsNumber; j++)
683  {
684  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(ACCELERATION_X);
685  }
686  itNodeNew->FastGetSolutionStepValue(ACCELERATION_X) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
687  }
688  if( itNodeNew->IsFixed(DISPLACEMENT_Y)==false )
689  {
690  for(int j = 0; j < PointsNumber; j++)
691  {
692  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DISPLACEMENT_Y);
693  }
694  itNodeNew->FastGetSolutionStepValue(DISPLACEMENT_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
695  }
696  if( itNodeNew->IsFixed(VELOCITY_Y)==false )
697  {
698  for(int j = 0; j < PointsNumber; j++)
699  {
700  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(VELOCITY_Y);
701  }
702  itNodeNew->FastGetSolutionStepValue(VELOCITY_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
703  }
704  if( itNodeNew->IsFixed(ACCELERATION_Y)==false )
705  {
706  for(int j = 0; j < PointsNumber; j++)
707  {
708  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(ACCELERATION_Y);
709  }
710  itNodeNew->FastGetSolutionStepValue(ACCELERATION_Y) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
711  }
712  if( itNodeNew->IsFixed(WATER_PRESSURE)==false )
713  {
714  for(int j = 0; j < PointsNumber; j++)
715  {
716  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(WATER_PRESSURE);
717  }
718  itNodeNew->FastGetSolutionStepValue(WATER_PRESSURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
719  for(int j = 0; j < PointsNumber; j++)
720  {
721  NodalVariableVector[j] = pElementOld->GetGeometry().GetPoint(j).FastGetSolutionStepValue(DT_WATER_PRESSURE);
722  }
723  itNodeNew->FastGetSolutionStepValue(DT_WATER_PRESSURE) = inner_prod(ShapeFunctionsValuesVector,NodalVariableVector);
724  }
725  }
726  }
727 
729 
730  void GaussPointStateVariableMapping(
731  const UtilityVariables& AuxVariables,
732  Parameters& rParameters,
733  ModelPart& rModelPartOld,
734  ModelPart& rModelPartNew)
735  {
736  // Define GaussPointOld Cell matrix
737  std::vector< std::vector< std::vector<GaussPointOld> > > BodyGaussPointOldCellMatrix;
738  BodyGaussPointOldCellMatrix.resize(AuxVariables.NRows);
739  for(int i = 0; i < AuxVariables.NRows; i++) BodyGaussPointOldCellMatrix[i].resize(AuxVariables.NColumns);
740 
741  std::vector< std::vector< std::vector<GaussPointOld> > > InterfaceGaussPointOldCellMatrix;
742  InterfaceGaussPointOldCellMatrix.resize(AuxVariables.NRows);
743  for(int i = 0; i < AuxVariables.NRows; i++) InterfaceGaussPointOldCellMatrix[i].resize(AuxVariables.NColumns);
744 
745  // Locate Old Gauss Points in cells
746  GaussPointOld MyGaussPointOld;
747  GeometryData::IntegrationMethod MyIntegrationMethod;
748  const ProcessInfo& CurrentProcessInfoOld = rModelPartOld.GetProcessInfo();
749  array_1d<double,3> AuxLocalCoordinates;
750 
751  unsigned int NumBodySubModelParts = rParameters["fracture_data"]["body_domain_sub_model_part_list"].size();
752 
753  // Loop through all OLD BodySubModelParts
754  for(unsigned int i = 0; i < NumBodySubModelParts; i++)
755  {
756  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][i].GetString());
757 
758  int NElems = static_cast<int>(SubModelPart.Elements().size());
759  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
760 
761  #pragma omp parallel for private(MyGaussPointOld,MyIntegrationMethod,AuxLocalCoordinates)
762  for(int j = 0; j < NElems; j++)
763  {
764  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
765 
766  Element::GeometryType& rGeom = itElem->GetGeometry();
767  MyIntegrationMethod = itElem->GetIntegrationMethod();
768  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
769  unsigned int NumGPoints = IntegrationPoints.size();
770  Vector detJContainer(NumGPoints);
771  rGeom.DeterminantOfJacobian(detJContainer,MyIntegrationMethod);
772  std::vector<double> StateVariableVector(NumGPoints);
773  itElem->CalculateOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoOld);
774  int Row;
775  int Column;
776 
777  // Loop through GaussPoints
778  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
779  {
780  // GaussPointOld Coordinates
781  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
782  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
783  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
784  rGeom.GlobalCoordinates(MyGaussPointOld.Coordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
785 
786  // GaussPointOld Weight
787  MyGaussPointOld.Weight = detJContainer[GPoint]*IntegrationPoints[GPoint].Weight();
788 
789  // GaussPointOld StateVariable and Damage
790  MyGaussPointOld.StateVariable = StateVariableVector[GPoint];
791 
792  // GaussPointOld Row and Column
793  Row = int((AuxVariables.Y_max-MyGaussPointOld.Coordinates[1])/AuxVariables.RowSize);
794  Column = int((MyGaussPointOld.Coordinates[0]-AuxVariables.X_min)/AuxVariables.ColumnSize);
795  #pragma omp critical
796  {
797  BodyGaussPointOldCellMatrix[Row][Column].push_back(MyGaussPointOld);
798  }
799  }
800  }
801  }
802 
803  unsigned int NumInterfaceSubModelPartsOld = rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"].size();
804 
805  // Loop through all OLD InterfaceSubModelParts
806  for(unsigned int i = 0; i < NumInterfaceSubModelPartsOld; i++)
807  {
808  ModelPart& SubModelPart = rModelPartOld.GetSubModelPart(rParameters["fracture_data"]["interface_domain_sub_model_part_old_list"][i].GetString());
809 
810  int NElems = static_cast<int>(SubModelPart.Elements().size());
811  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
812 
813  #pragma omp parallel for private(MyGaussPointOld,MyIntegrationMethod,AuxLocalCoordinates)
814  for(int j = 0; j < NElems; j++)
815  {
816  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
817 
818  Element::GeometryType& rGeom = itElem->GetGeometry();
819  MyIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1;
820  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
821  unsigned int NumGPoints = IntegrationPoints.size();
822  Vector detJContainer(NumGPoints);
823  rGeom.DeterminantOfJacobian(detJContainer,MyIntegrationMethod);
824  std::vector<double> StateVariableVector(NumGPoints);
825  itElem->CalculateOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoOld);
826  int Row;
827  int Column;
828 
829  // Loop through GaussPoints
830  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
831  {
832  // GaussPointOld Coordinates
833  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
834  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
835  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
836  rGeom.GlobalCoordinates(MyGaussPointOld.Coordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
837 
838  // GaussPointOld Weight
839  MyGaussPointOld.Weight = detJContainer[GPoint]*IntegrationPoints[GPoint].Weight();
840 
841  // GaussPointOld StateVariable
842  MyGaussPointOld.StateVariable = StateVariableVector[GPoint];
843 
844  // GaussPointOld Row and Column
845  Row = int((AuxVariables.Y_max-MyGaussPointOld.Coordinates[1])/AuxVariables.RowSize);
846  Column = int((MyGaussPointOld.Coordinates[0]-AuxVariables.X_min)/AuxVariables.ColumnSize);
847  #pragma omp critical
848  {
849  InterfaceGaussPointOldCellMatrix[Row][Column].push_back(MyGaussPointOld);
850  }
851  }
852  }
853  }
854 
855  // Transfer state variables from old Gauss points to new Gauss Points (nonlocal average)
856  const ProcessInfo& CurrentProcessInfoNew = rModelPartNew.GetProcessInfo();
857  const double PropagationLength = rParameters["fracture_data"]["propagation_length"].GetDouble();
858  array_1d<double,3> AuxGlobalCoordinates;
859 
860  // Loop through all NEW BodySubModelParts
861  for(unsigned int i = 0; i < NumBodySubModelParts; i++)
862  {
863  ModelPart& SubModelPart = rModelPartNew.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][i].GetString());
864 
865  int NElems = static_cast<int>(SubModelPart.Elements().size());
866  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
867 
868  double DamageThreshold = el_begin->GetProperties()[DAMAGE_THRESHOLD];
869 
870  #pragma omp parallel for private(MyIntegrationMethod,AuxLocalCoordinates,AuxGlobalCoordinates)
871  for(int j = 0; j < NElems; j++)
872  {
873  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
874 
875  Element::GeometryType& rGeom = itElem->GetGeometry();
876  MyIntegrationMethod = itElem->GetIntegrationMethod();
877  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
878  unsigned int NumGPoints = IntegrationPoints.size();
879  std::vector<double> StateVariableVector(NumGPoints);
880 
881  // Loop through NEW GaussPoints
882  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
883  {
884  // GaussPointNew Coordinates
885  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
886  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
887  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
888  rGeom.GlobalCoordinates(AuxGlobalCoordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
889  double X_me = AuxGlobalCoordinates[0];
890  double Y_me = AuxGlobalCoordinates[1];
891 
892  // GaussPointNew search area
893  double X_left = X_me - PropagationLength;
894  double X_right = X_me + PropagationLength;
895  double Y_top = Y_me + PropagationLength;
896  double Y_bot = Y_me - PropagationLength;
897 
898  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
899  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
900  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
901  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
902 
903  if(Column_left < 0) Column_left = 0;
904  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
905  if(Row_top < 0) Row_top = 0;
906  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
907 
908  // Search GaussPointOld neighbours around the GaussPointNew and compute nonlocal state variable
909  double Numerator = 0.0;
910  double WeightingFunctionDenominator = 0.0;
911  double Distance;
912  for(int k = Row_top; k <= Row_bot; k++)
913  {
914  for(int l = Column_left; l<= Column_right; l++)
915  {
916  for(unsigned int m = 0; m < BodyGaussPointOldCellMatrix[k][l].size(); m++)
917  {
918  GaussPointOld& rOtherGaussPointOld = BodyGaussPointOldCellMatrix[k][l][m];
919 
920  Distance = sqrt((rOtherGaussPointOld.Coordinates[0]-X_me)*(rOtherGaussPointOld.Coordinates[0]-X_me) +
921  (rOtherGaussPointOld.Coordinates[1]-Y_me)*(rOtherGaussPointOld.Coordinates[1]-Y_me));
922 
923  if(Distance <= PropagationLength)
924  {
925  Numerator += rOtherGaussPointOld.Weight
926  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength))
927  *rOtherGaussPointOld.StateVariable;
928  WeightingFunctionDenominator += rOtherGaussPointOld.Weight
929  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength));
930  }
931  }
932  }
933  }
934 
935  // Save computed stateVariable
936  if(WeightingFunctionDenominator > 0.0)
937  StateVariableVector[GPoint] = Numerator/WeightingFunctionDenominator;
938  else
939  StateVariableVector[GPoint] = DamageThreshold;
940  }
941  // Set stateVariable of new GaussPoints
942  itElem->SetValuesOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoNew);
943  }
944  }
945 
946  unsigned int NumInterfaceSubModelParts = rParameters["fracture_data"]["interface_domain_sub_model_part_list"].size();
947  const double PropagationDamage = rParameters["fracture_data"]["propagation_damage"].GetDouble();
948 
949  // Loop through all NEW InterfaceSubModelParts
950  for(unsigned int i = 0; i < NumInterfaceSubModelParts; i++)
951  {
952  ModelPart& SubModelPart = rModelPartNew.GetSubModelPart(rParameters["fracture_data"]["interface_domain_sub_model_part_list"][i].GetString());
953 
954  int NElems = static_cast<int>(SubModelPart.Elements().size());
955  ModelPart::ElementsContainerType::iterator el_begin = SubModelPart.ElementsBegin();
956 
957  //double DamageThreshold = el_begin->GetProperties()[DAMAGE_THRESHOLD];
958 
959  #pragma omp parallel for private(MyIntegrationMethod,AuxLocalCoordinates,AuxGlobalCoordinates)
960  for(int j = 0; j < NElems; j++)
961  {
962  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
963 
964  Element::GeometryType& rGeom = itElem->GetGeometry();
965  MyIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1;
966  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
967  unsigned int NumGPoints = IntegrationPoints.size();
968  std::vector<double> StateVariableVector(NumGPoints);
969 
970  // Loop through NEW GaussPoints
971  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
972  {
973  // GaussPointNew Coordinates
974  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
975  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
976  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
977  rGeom.GlobalCoordinates(AuxGlobalCoordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
978  double X_me = AuxGlobalCoordinates[0];
979  double Y_me = AuxGlobalCoordinates[1];
980 
981  // GaussPointNew search area
982  double X_left = X_me - PropagationLength;
983  double X_right = X_me + PropagationLength;
984  double Y_top = Y_me + PropagationLength;
985  double Y_bot = Y_me - PropagationLength;
986 
987  int Column_left = int((X_left-AuxVariables.X_min)/AuxVariables.ColumnSize);
988  int Column_right = int((X_right-AuxVariables.X_min)/AuxVariables.ColumnSize);
989  int Row_top = int((AuxVariables.Y_max-Y_top)/AuxVariables.RowSize);
990  int Row_bot = int((AuxVariables.Y_max-Y_bot)/AuxVariables.RowSize);
991 
992  if(Column_left < 0) Column_left = 0;
993  if(Column_right >= AuxVariables.NColumns) Column_right = AuxVariables.NColumns-1;
994  if(Row_top < 0) Row_top = 0;
995  if(Row_bot >= AuxVariables.NRows) Row_bot = AuxVariables.NRows-1;
996 
997  // Search GaussPointOld neighbours around the GaussPointNew and compute nonlocal state variable
998  double Numerator = 0.0;
999  double WeightingFunctionDenominator = 0.0;
1000  double Distance;
1001  for(int k = Row_top; k <= Row_bot; k++)
1002  {
1003  for(int l = Column_left; l<= Column_right; l++)
1004  {
1005  for(unsigned int m = 0; m < InterfaceGaussPointOldCellMatrix[k][l].size(); m++)
1006  {
1007  GaussPointOld& rOtherGaussPointOld = InterfaceGaussPointOldCellMatrix[k][l][m];
1008 
1009  Distance = sqrt((rOtherGaussPointOld.Coordinates[0]-X_me)*(rOtherGaussPointOld.Coordinates[0]-X_me) +
1010  (rOtherGaussPointOld.Coordinates[1]-Y_me)*(rOtherGaussPointOld.Coordinates[1]-Y_me));
1011 
1012  if(Distance <= PropagationLength)
1013  {
1014  Numerator += rOtherGaussPointOld.Weight
1015  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength))
1016  *rOtherGaussPointOld.StateVariable;
1017  WeightingFunctionDenominator += rOtherGaussPointOld.Weight
1018  *exp(-4.0*Distance*Distance/(PropagationLength*PropagationLength));
1019  }
1020  }
1021  }
1022  }
1023 
1024  // Save computed stateVariable
1025  if(WeightingFunctionDenominator > 0.0)
1026  {
1027  StateVariableVector[GPoint] = Numerator/WeightingFunctionDenominator;
1028  if(StateVariableVector[GPoint] < PropagationDamage)
1029  {
1030  StateVariableVector[GPoint] = PropagationDamage;
1031  }
1032  }
1033  else
1034  StateVariableVector[GPoint] = PropagationDamage;
1035  }
1036  // Set stateVariable of new GaussPoints
1037  itElem->SetValuesOnIntegrationPoints(STATE_VARIABLE,StateVariableVector,CurrentProcessInfoNew);
1038  }
1039  }
1040  }
1041 
1043 
1044  void ResetMeshPosition(ModelPart& rModelPart)
1045  {
1046  // Move mesh to the original position
1047  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
1048  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
1049 
1050  #pragma omp parallel for
1051  for(int i = 0; i < NNodes; i++)
1052  {
1053  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
1054 
1055  itNode->X() = itNode->X0();
1056  itNode->Y() = itNode->Y0();
1057  itNode->Z() = itNode->Z0();
1058  }
1059  }
1060 
1062 
1063  void UpdateMeshPosition(ModelPart& rModelPart)
1064  {
1065  // Move mesh to the current position
1066  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
1067  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
1068 
1069  #pragma omp parallel for
1070  for(int i = 0; i < NNodes; i++)
1071  {
1072  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
1073 
1074  itNode->X() = itNode->X0() + itNode->FastGetSolutionStepValue(DISPLACEMENT_X);
1075  itNode->Y() = itNode->Y0() + itNode->FastGetSolutionStepValue(DISPLACEMENT_Y);
1076  itNode->Z() = itNode->Z0() + itNode->FastGetSolutionStepValue(DISPLACEMENT_Z);
1077  }
1078  }
1079 
1080 private:
1081 
1083 
1084  void ComputeCellMatrixDimensions(
1085  UtilityVariables& rAuxVariables,
1086  ModelPart& rModelPart)
1087  {
1088  // Compute X and Y limits of the current geometry
1089  unsigned int NumThreads = ParallelUtilities::GetNumThreads();
1090  std::vector<double> X_max_partition(NumThreads);
1091  std::vector<double> X_min_partition(NumThreads);
1092  std::vector<double> Y_max_partition(NumThreads);
1093  std::vector<double> Y_min_partition(NumThreads);
1094 
1095  const int NNodes = static_cast<int>(rModelPart.Nodes().size());
1096  ModelPart::NodesContainerType::iterator node_begin = rModelPart.NodesBegin();
1097 
1098  #pragma omp parallel
1099  {
1100  int k = OpenMPUtils::ThisThread();
1101 
1102  X_max_partition[k] = node_begin->X();
1103  X_min_partition[k] = X_max_partition[k];
1104  Y_max_partition[k] = node_begin->Y();
1105  Y_min_partition[k] = Y_max_partition[k];
1106 
1107  double X_me, Y_me;
1108 
1109  #pragma omp for
1110  for(int i = 0; i < NNodes; i++)
1111  {
1112  ModelPart::NodesContainerType::iterator itNode = node_begin + i;
1113 
1114  X_me = itNode->X();
1115  Y_me = itNode->Y();
1116 
1117  if( X_me > X_max_partition[k] ) X_max_partition[k] = X_me;
1118  else if( X_me < X_min_partition[k] ) X_min_partition[k] = X_me;
1119 
1120  if( Y_me > Y_max_partition[k] ) Y_max_partition[k] = Y_me;
1121  else if( Y_me < Y_min_partition[k] ) Y_min_partition[k] = Y_me;
1122  }
1123  }
1124 
1125  rAuxVariables.X_max = X_max_partition[0];
1126  rAuxVariables.X_min = X_min_partition[0];
1127  rAuxVariables.Y_max = Y_max_partition[0];
1128  rAuxVariables.Y_min = Y_min_partition[0];
1129 
1130  for(unsigned int i=1; i < NumThreads; i++)
1131  {
1132  if(X_max_partition[i] > rAuxVariables.X_max) rAuxVariables.X_max = X_max_partition[i];
1133  if(X_min_partition[i] < rAuxVariables.X_min) rAuxVariables.X_min = X_min_partition[i];
1134  if(Y_max_partition[i] > rAuxVariables.Y_max) rAuxVariables.Y_max = Y_max_partition[i];
1135  if(Y_min_partition[i] < rAuxVariables.Y_min) rAuxVariables.Y_min = Y_min_partition[i];
1136  }
1137 
1138  // Calculate Average Element Length
1139  double AverageElementLength = 0.0;
1140 
1141  int NElems = static_cast<int>(rModelPart.Elements().size());
1142  ModelPart::ElementsContainerType::iterator el_begin = rModelPart.ElementsBegin();
1143 
1144  #pragma omp parallel for reduction(+:AverageElementLength)
1145  for(int i = 0; i < NElems; i++)
1146  {
1147  ModelPart::ElementsContainerType::iterator itElem = el_begin + i;
1148 
1149  AverageElementLength += itElem->GetGeometry().Length();
1150  }
1151 
1152  AverageElementLength = AverageElementLength/NElems;
1153 
1154  // Compute FracturePoints CellMatrix dimensions
1155  rAuxVariables.NRows = int((rAuxVariables.Y_max-rAuxVariables.Y_min)/(AverageElementLength));
1156  rAuxVariables.NColumns = int((rAuxVariables.X_max-rAuxVariables.X_min)/(AverageElementLength));
1157  if(rAuxVariables.NRows <= 0) rAuxVariables.NRows = 1;
1158  if(rAuxVariables.NColumns <= 0) rAuxVariables.NColumns = 1;
1159  rAuxVariables.RowSize = (rAuxVariables.Y_max-rAuxVariables.Y_min)/rAuxVariables.NRows;
1160  rAuxVariables.ColumnSize = (rAuxVariables.X_max-rAuxVariables.X_min)/rAuxVariables.NColumns;
1161  }
1162 
1164 
1165  void SetFracturePoints(
1166  PropagationGlobalVariables& rPropagationData,
1167  UtilityVariables& rAuxVariables,
1168  Parameters& rParameters,
1169  ModelPart& rModelPart)
1170  {
1171  // Compute FracturePointsCellMatrix dimensions
1172  this->ComputeCellMatrixDimensions(rAuxVariables,rModelPart);
1173 
1174  rPropagationData.FracturePointsCellMatrix.resize(rAuxVariables.NRows);
1175  for(int i = 0; i < rAuxVariables.NRows; i++) rPropagationData.FracturePointsCellMatrix[i].resize(rAuxVariables.NColumns);
1176 
1177  // Locate FracturePoints inside CellMatrix
1178  FracturePoint MyFracturePoint;
1179  GeometryData::IntegrationMethod MyIntegrationMethod;
1180  const ProcessInfo& CurrentProcessInfo = rModelPart.GetProcessInfo();
1181  rPropagationData.pProcessInfo = rModelPart.pGetProcessInfo();
1182  array_1d<double,3> AuxLocalCoordinates;
1183  array_1d<double,3> AuxGlobalCoordinates;
1184 
1185  unsigned int NumBodySubModelParts = rParameters["fracture_data"]["body_domain_sub_model_part_list"].size();
1186 
1187  // Loop through all BodySubModelParts
1188  for(unsigned int i = 0; i < NumBodySubModelParts; i++)
1189  {
1190  ModelPart& BodySubModelPart = rModelPart.GetSubModelPart(rParameters["fracture_data"]["body_domain_sub_model_part_list"][i].GetString());
1191 
1192  int NElems = static_cast<int>(BodySubModelPart.Elements().size());
1193  ModelPart::ElementsContainerType::iterator el_begin = BodySubModelPart.ElementsBegin();
1194 
1195  // Loop through all body elements
1196  #pragma omp parallel for private(MyFracturePoint,MyIntegrationMethod,AuxLocalCoordinates,AuxGlobalCoordinates)
1197  for(int j = 0; j < NElems; j++)
1198  {
1199  ModelPart::ElementsContainerType::iterator itElem = el_begin + j;
1200 
1201  Element::GeometryType& rGeom = itElem->GetGeometry();
1202  MyIntegrationMethod = itElem->GetIntegrationMethod();
1203  const Element::GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(MyIntegrationMethod);
1204  unsigned int NumGPoints = IntegrationPoints.size();
1205  Vector detJContainer(NumGPoints);
1206  rGeom.DeterminantOfJacobian(detJContainer,MyIntegrationMethod);
1207  std::vector<double> DamageVector(NumGPoints);
1208  itElem->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1209  int Row;
1210  int Column;
1211 
1212  // Loop through GaussPoints
1213  for ( unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++ )
1214  {
1215  // FracturePoint Coordinates
1216  AuxLocalCoordinates[0] = IntegrationPoints[GPoint][0];
1217  AuxLocalCoordinates[1] = IntegrationPoints[GPoint][1];
1218  AuxLocalCoordinates[2] = IntegrationPoints[GPoint][2];
1219  rGeom.GlobalCoordinates(AuxGlobalCoordinates,AuxLocalCoordinates); //Note: these are the CURRENT global coordinates
1220  MyFracturePoint.Coordinates[0] = AuxGlobalCoordinates[0];
1221  MyFracturePoint.Coordinates[1] = AuxGlobalCoordinates[1];
1222 
1223  // FracturePoint Weight
1224  MyFracturePoint.Weight = detJContainer[GPoint]*IntegrationPoints[GPoint].Weight();
1225 
1226  // FracturePoint Damage
1227  MyFracturePoint.Damage = DamageVector[GPoint];
1228 
1229  // FracturePoint Row and Column
1230  Row = int((rAuxVariables.Y_max-MyFracturePoint.Coordinates[1])/rAuxVariables.RowSize);
1231  Column = int((MyFracturePoint.Coordinates[0]-rAuxVariables.X_min)/rAuxVariables.ColumnSize);
1232 
1233  // Element containing the FracturePoint
1234  MyFracturePoint.pElement = (*(itElem.base()));
1235 
1236  #pragma omp critical
1237  {
1238  rPropagationData.FracturePointsCellMatrix[Row][Column].push_back(MyFracturePoint);
1239  }
1240  }
1241  }
1242  }
1243  }
1244 
1246 
1247  void PropagateFracture(
1248  const unsigned int& itFracture,
1249  PropagationGlobalVariables& rPropagationData,
1250  PropagationLocalVariables& rAuxPropagationVariables,
1251  Parameters& rParameters)
1252  {
1253  // Compute Propagation Coordinates
1254  double TipX = 0.0;
1255  double TipY = 0.0;
1256  double TipDen = 0.0;
1257 
1258  #pragma omp parallel sections reduction(+:TipX,TipY,TipDen)
1259  {
1260  #pragma omp section
1261  {
1262  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1263  {
1264  this->AverageTipCoordinates(TipX,TipY,TipDen,*(rAuxPropagationVariables.TopFrontFracturePoints[i]));
1265  }
1266  }
1267  #pragma omp section
1268  {
1269  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1270  {
1271  this->AverageTipCoordinates(TipX,TipY,TipDen,*(rAuxPropagationVariables.BotFrontFracturePoints[i]));
1272  }
1273  }
1274  }
1275 
1276  array_1d<double,2> AuxArray1;
1277  array_1d<double,2> AuxArray2;
1278  const double PropagationWidth = rParameters["fracture_data"]["propagation_width"].GetDouble();
1279  int MotherFractureId = rParameters["fractures_list"][itFracture]["id"].GetInt();
1280  Propagation MyPropagation;
1281 
1282  MyPropagation.MotherFractureId = MotherFractureId;
1283 
1284  MyPropagation.TipCoordinates[0] = TipX/TipDen;
1285  MyPropagation.TipCoordinates[1] = TipY/TipDen;
1286  MyPropagation.TipCoordinates[2] = 0.0;
1287 
1288  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1289  AuxArray1[1] += 0.5*PropagationWidth;
1290  noalias(AuxArray2) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1291  MyPropagation.TopInitCoordinates[0] = AuxArray2[0];
1292  MyPropagation.TopInitCoordinates[1] = AuxArray2[1];
1293  MyPropagation.TopInitCoordinates[2] = 0.0;
1294 
1295  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1296  AuxArray1[1] -= 0.5*PropagationWidth;
1297  noalias(AuxArray2) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1298  MyPropagation.BotInitCoordinates[0] = AuxArray2[0];
1299  MyPropagation.BotInitCoordinates[1] = AuxArray2[1];
1300  MyPropagation.BotInitCoordinates[2] = 0.0;
1301 
1302  // Check straight propagation
1303  const double PropagationLength = rParameters["fracture_data"]["propagation_length"].GetDouble();
1304  const double CorrectionTol = rParameters["fracture_data"]["correction_tolerance"].GetDouble();
1305  AuxArray2[0] = MyPropagation.TipCoordinates[0];
1306  AuxArray2[1] = MyPropagation.TipCoordinates[1];
1307  noalias(AuxArray1) = prod(rAuxPropagationVariables.RotationMatrix,AuxArray2);
1308  AuxArray1[1] = rAuxPropagationVariables.TipLocalCoordinates[1];
1309  noalias(AuxArray2) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1310  double Distance = sqrt((MyPropagation.TipCoordinates[0]-AuxArray2[0])*(MyPropagation.TipCoordinates[0]-AuxArray2[0])+
1311  (MyPropagation.TipCoordinates[1]-AuxArray2[1])*(MyPropagation.TipCoordinates[1]-AuxArray2[1]));
1312  if (Distance <= PropagationLength*CorrectionTol)
1313  {
1314  MyPropagation.TipCoordinates[0] = AuxArray2[0];
1315  MyPropagation.TipCoordinates[1] = AuxArray2[1];
1316  }
1317 
1318  // Check whether new tip falls inside a valid element
1319  array_1d<double,3> LocalCoordinates;
1320  Element::Pointer pElement;
1321  bool IsInside = false;
1322 
1323  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1324  {
1325  pElement = rAuxPropagationVariables.TopFrontFracturePoints[i]->pElement;
1326  IsInside = pElement->GetGeometry().IsInside(MyPropagation.TipCoordinates,LocalCoordinates);
1327  if(IsInside) break;
1328  }
1329  if(IsInside == false)
1330  {
1331  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1332  {
1333  pElement = rAuxPropagationVariables.BotFrontFracturePoints[i]->pElement;
1334  IsInside = pElement->GetGeometry().IsInside(MyPropagation.TipCoordinates,LocalCoordinates);
1335  if(IsInside) break;
1336  }
1337  }
1338 
1339  const double PropagationDamage = rParameters["fracture_data"]["propagation_damage"].GetDouble();
1340  const ProcessInfo& CurrentProcessInfo = *(rPropagationData.pProcessInfo);
1341 
1342  if (IsInside == true)
1343  {
1344  std::vector<double> DamageVector;
1345  pElement->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1346  unsigned int NumGPoints = DamageVector.size();
1347  double InvNumGPoints = 1.0/static_cast<double>(NumGPoints);
1348  double ElementDamage = 0.0;
1349  for (unsigned int i = 0; i < NumGPoints; i++)
1350  {
1351  ElementDamage += DamageVector[i];
1352  }
1353  ElementDamage *= InvNumGPoints;
1354  if (ElementDamage >= 0.5*PropagationDamage)
1355  {
1356  rPropagationData.PropagationVector.push_back(MyPropagation);
1357  rPropagationData.PropagateFractures = true;
1358  return;
1359  }
1360  }
1361 
1362  // Being here means that the new tip does not fall inside a valid element. We need to check top and bot fracture points
1363  bool PropagateTop = false;
1364  bool PropagateBot = false;
1365  double TopEndX = 0.0;
1366  double TopEndY = 0.0;
1367  double TopEndDen = 0.0;
1368  double BotEndX = 0.0;
1369  double BotEndY = 0.0;
1370  double BotEndDen = 0.0;
1371  array_1d<double,3> GlobalCoordinates;
1372 
1373  #pragma omp parallel sections private(GlobalCoordinates,LocalCoordinates,pElement,IsInside)
1374  {
1375  #pragma omp section
1376  {
1377  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1378  {
1379  this->AverageTipCoordinates(TopEndX,TopEndY,TopEndDen,*(rAuxPropagationVariables.TopFrontFracturePoints[i]));
1380  }
1381  GlobalCoordinates[0] = TopEndX/TopEndDen;
1382  GlobalCoordinates[1] = TopEndY/TopEndDen;
1383  GlobalCoordinates[2] = 0.0;
1384 
1385  // Check whether new tip falls inside a valid element
1386  IsInside = false;
1387  for(unsigned int i = 0; i < rAuxPropagationVariables.TopFrontFracturePoints.size(); i++)
1388  {
1389  pElement = rAuxPropagationVariables.TopFrontFracturePoints[i]->pElement;
1390  IsInside = pElement->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates);
1391  if(IsInside) break;
1392  }
1393 
1394  if (IsInside == true)
1395  {
1396  std::vector<double> DamageVector;
1397  pElement->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1398  unsigned int NumGPoints = DamageVector.size();
1399  double InvNumGPoints = 1.0/static_cast<double>(NumGPoints);
1400  double ElementDamage = 0.0;
1401  for (unsigned int i = 0; i < NumGPoints; i++)
1402  {
1403  ElementDamage += DamageVector[i];
1404  }
1405  ElementDamage *= InvNumGPoints;
1406  if (ElementDamage >= PropagationDamage)
1407  PropagateTop = true;
1408  }
1409  }
1410  #pragma omp section
1411  {
1412  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1413  {
1414  this->AverageTipCoordinates(BotEndX,BotEndY,BotEndDen,*(rAuxPropagationVariables.BotFrontFracturePoints[i]));
1415  }
1416  GlobalCoordinates[0] = BotEndX/BotEndDen;
1417  GlobalCoordinates[1] = BotEndY/BotEndDen;
1418  GlobalCoordinates[2] = 0.0;
1419 
1420  // Check whether new tip falls inside a valid element
1421  IsInside = false;
1422  for(unsigned int i = 0; i < rAuxPropagationVariables.BotFrontFracturePoints.size(); i++)
1423  {
1424  pElement = rAuxPropagationVariables.BotFrontFracturePoints[i]->pElement;
1425  IsInside = pElement->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates);
1426  if(IsInside) break;
1427  }
1428 
1429  if (IsInside == true)
1430  {
1431  std::vector<double> DamageVector;
1432  pElement->CalculateOnIntegrationPoints(DAMAGE_VARIABLE,DamageVector,CurrentProcessInfo);
1433  unsigned int NumGPoints = DamageVector.size();
1434  double InvNumGPoints = 1.0/static_cast<double>(NumGPoints);
1435  double ElementDamage = 0.0;
1436  for (unsigned int i = 0; i < NumGPoints; i++)
1437  {
1438  ElementDamage += DamageVector[i];
1439  }
1440  ElementDamage *= InvNumGPoints;
1441  if (ElementDamage >= PropagationDamage)
1442  PropagateBot = true;
1443  }
1444  }
1445  }
1446 
1447  if (PropagateTop == true && PropagateBot == true) // Bifurcation
1448  {
1449  Bifurcation MyBifurcation;
1450  MyBifurcation.MotherFractureId = MotherFractureId;
1451 
1452  MyBifurcation.TopTipCoordinates[0] = TopEndX/TopEndDen;
1453  MyBifurcation.TopTipCoordinates[1] = TopEndY/TopEndDen;
1454  MyBifurcation.TopTipCoordinates[2] = 0.0;
1455 
1456  MyBifurcation.BotTipCoordinates[0] = BotEndX/BotEndDen;
1457  MyBifurcation.BotTipCoordinates[1] = BotEndY/BotEndDen;
1458  MyBifurcation.BotTipCoordinates[2] = 0.0;
1459 
1460  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1461  AuxArray1[0] -= PropagationWidth;
1462  AuxArray1[1] += 0.5*PropagationWidth;
1463  noalias(AuxArray2) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1464  MyBifurcation.TopInitCoordinates[0] = AuxArray2[0];
1465  MyBifurcation.TopInitCoordinates[1] = AuxArray2[1];
1466  MyBifurcation.TopInitCoordinates[2] = 0.0;
1467 
1468  noalias(AuxArray1) = rAuxPropagationVariables.TipLocalCoordinates;
1469  AuxArray1[0] -= PropagationWidth;
1470  AuxArray1[1] -= 0.5*PropagationWidth;
1471  noalias(AuxArray2) = prod(trans(rAuxPropagationVariables.RotationMatrix),AuxArray1);
1472  MyBifurcation.BotInitCoordinates[0] = AuxArray2[0];
1473  MyBifurcation.BotInitCoordinates[1] = AuxArray2[1];
1474  MyBifurcation.BotInitCoordinates[2] = 0.0;
1475 
1476  rPropagationData.BifurcationVector.push_back(MyBifurcation);
1477  rPropagationData.PropagateFractures = true;
1478  return;
1479  }
1480  else if (PropagateTop == true) // Top Propagation
1481  {
1482  MyPropagation.TipCoordinates[0] = TopEndX/TopEndDen;
1483  MyPropagation.TipCoordinates[1] = TopEndY/TopEndDen;
1484  MyPropagation.TipCoordinates[2] = 0.0;
1485 
1486  rPropagationData.PropagationVector.push_back(MyPropagation);
1487  rPropagationData.PropagateFractures = true;
1488  return;
1489  }
1490  else if (PropagateBot == true) // Bot Propagation
1491  {
1492  MyPropagation.TipCoordinates[0] = BotEndX/BotEndDen;
1493  MyPropagation.TipCoordinates[1] = BotEndY/BotEndDen;
1494  MyPropagation.TipCoordinates[2] = 0.0;
1495 
1496  rPropagationData.PropagationVector.push_back(MyPropagation);
1497  rPropagationData.PropagateFractures = true;
1498  return;
1499  }
1500  }
1501 
1503 
1504  void CalculateTipRotationMatrix(
1505  const unsigned int& itFracture,
1506  BoundedMatrix<double,2,2>& rRotationMatrix,
1507  Parameters& rParameters)
1508  {
1509  array_1d<double,3> TipPointCoordinates;
1510  array_1d<double,3> TopPointCoordinates;
1511  array_1d<double,3> BotPointCoordinates;
1512  for(unsigned int i = 0; i < 3; i++)
1513  {
1514  TipPointCoordinates[i] = rParameters["fractures_list"][itFracture]["tip_point"]["coordinates"][i].GetDouble();
1515  TopPointCoordinates[i] = rParameters["fractures_list"][itFracture]["top_point"]["coordinates"][i].GetDouble();
1516  BotPointCoordinates[i] = rParameters["fractures_list"][itFracture]["bot_point"]["coordinates"][i].GetDouble();
1517  }
1518 
1519  array_1d<double, 3> pmid0;
1520  noalias(pmid0) = 0.5 * (TopPointCoordinates + BotPointCoordinates);
1521 
1522  array_1d<double, 3> Vx;
1523  noalias(Vx) = TipPointCoordinates - pmid0;
1524  double inv_norm_x = 1.0/norm_2(Vx);
1525  Vx[0] *= inv_norm_x;
1526  Vx[1] *= inv_norm_x;
1527 
1528  //Rotation Matrix
1529  // We need to determine the unitary vector in local y direction pointing towards the TOP face of the joint
1530  rRotationMatrix(0,0) = Vx[0];
1531  rRotationMatrix(0,1) = Vx[1];
1532 
1533  // NOTE. Assuming that the nodes in quadrilateral_interface_2d_4 are
1534  // ordered counter-clockwise (GiD does so now), the rotation matrix is build like follows:
1535  rRotationMatrix(1,0) = -Vx[1];
1536  rRotationMatrix(1,1) = Vx[0];
1537  // NOTE. Assuming that the nodes in quadrilateral_interface_2d_4 are
1538  // ordered clockwise, the rotation matrix is build like follows:
1539  // rRotationMatrix(1,0) = Vx[1];
1540  // rRotationMatrix(1,1) = -Vx[0];
1541 
1542  // NOTE. In zero-thickness quadrilateral_interface_2d_4 elements we are not able to know
1543  // whether the nodes are ordered clockwise or counter-clockwise.
1544  }
1545 
1547 
1548  void AverageTipCoordinates(
1549  double& rTipX,
1550  double& rTipY,
1551  double& rTipDenominator,
1552  const FracturePoint& MyFracturePoint)
1553  {
1554  rTipX += MyFracturePoint.Weight * MyFracturePoint.Damage * (MyFracturePoint.Coordinates[0]);
1555  rTipY += MyFracturePoint.Weight * MyFracturePoint.Damage * (MyFracturePoint.Coordinates[1]);
1556  rTipDenominator += MyFracturePoint.Weight * MyFracturePoint.Damage;
1557  }
1558 
1559 }; // Class FracturePropagation2DUtilities
1560 
1561 } // namespace Kratos.
1562 
1563 #endif /* KRATOS_PROPAGATE_FRACTURES_2D_UTILITIES defined */
Kratos::Element::GeometryType
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: element.h:83
Kratos::FracturePropagation2DUtilities
Definition: fracture_propagation_2D_utilities.hpp:38
Kratos::FracturePropagation2DUtilities::CheckFracturePropagation
bool CheckFracturePropagation(Parameters &rParameters, ModelPart &rModelPart, const bool &move_mesh_flag)
Definition: fracture_propagation_2D_utilities.hpp:125
Kratos::FracturePropagation2DUtilities::InitializeMapping
void InitializeMapping(UtilityVariables &rAuxVariables, ModelPart &rModelPartNew, const bool &move_mesh_flag)
Mapping ---------------------------------------------------------------------------------------------...
Definition: fracture_propagation_2D_utilities.hpp:462
Kratos::FracturePropagation2DUtilities::~FracturePropagation2DUtilities
virtual ~FracturePropagation2DUtilities()
Destructor.
Definition: fracture_propagation_2D_utilities.hpp:121
Kratos::FracturePropagation2DUtilities::CheckFracture
void CheckFracture(const unsigned int &itFracture, PropagationGlobalVariables &rPropagationData, const UtilityVariables &AuxVariables, Parameters &rParameters)
Definition: fracture_propagation_2D_utilities.hpp:189
Kratos::FracturePropagation2DUtilities::InitializeCheckFracture
void InitializeCheckFracture(PropagationGlobalVariables &rPropagationData, UtilityVariables &rAuxVariables, Parameters &rParameters, ModelPart &rModelPart, const bool &move_mesh_flag)
Member Variables.
Definition: fracture_propagation_2D_utilities.hpp:169
Kratos::FracturePropagation2DUtilities::ResetMeshPosition
void ResetMeshPosition(ModelPart &rModelPart)
Common ----------------------------------------------------------------------------------------------...
Definition: fracture_propagation_2D_utilities.hpp:1044
Kratos::FracturePropagation2DUtilities::GaussPointStateVariableMapping
void GaussPointStateVariableMapping(const UtilityVariables &AuxVariables, Parameters &rParameters, ModelPart &rModelPartOld, ModelPart &rModelPartNew)
Definition: fracture_propagation_2D_utilities.hpp:730
Kratos::FracturePropagation2DUtilities::FracturePropagation2DUtilities
FracturePropagation2DUtilities()
Constructor.
Definition: fracture_propagation_2D_utilities.hpp:116
Kratos::FracturePropagation2DUtilities::WritePropagationData
void WritePropagationData(const PropagationGlobalVariables &PropagationData, Parameters &rParameters)
Definition: fracture_propagation_2D_utilities.hpp:310
Kratos::FracturePropagation2DUtilities::MappingModelParts
void MappingModelParts(Parameters &rParameters, ModelPart &rModelPartOld, ModelPart &rModelPartNew, const bool &move_mesh_flag)
Definition: fracture_propagation_2D_utilities.hpp:146
Kratos::FracturePropagation2DUtilities::FinalizeCheckFracture
void FinalizeCheckFracture(const PropagationGlobalVariables &PropagationData, Parameters &rParameters, ModelPart &rModelPart, const bool &move_mesh_flag)
Definition: fracture_propagation_2D_utilities.hpp:289
Kratos::FracturePropagation2DUtilities::NodalVariablesMapping
void NodalVariablesMapping(const UtilityVariables &AuxVariables, Parameters &rParameters, ModelPart &rModelPartOld, ModelPart &rModelPartNew)
Definition: fracture_propagation_2D_utilities.hpp:476
Kratos::FracturePropagation2DUtilities::UpdateMeshPosition
void UpdateMeshPosition(ModelPart &rModelPart)
Definition: fracture_propagation_2D_utilities.hpp:1063
Kratos::FracturePropagation2DUtilities::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(FracturePropagation2DUtilities)
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::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, 2 >
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::FracturePropagation2DUtilities::Bifurcation
Definition: fracture_propagation_2D_utilities.hpp:73
Kratos::FracturePropagation2DUtilities::Bifurcation::MotherFractureId
int MotherFractureId
Definition: fracture_propagation_2D_utilities.hpp:74
Kratos::FracturePropagation2DUtilities::Bifurcation::TopTipCoordinates
array_1d< double, 3 > TopTipCoordinates
Definition: fracture_propagation_2D_utilities.hpp:77
Kratos::FracturePropagation2DUtilities::Bifurcation::BotInitCoordinates
array_1d< double, 3 > BotInitCoordinates
Definition: fracture_propagation_2D_utilities.hpp:76
Kratos::FracturePropagation2DUtilities::Bifurcation::TopInitCoordinates
array_1d< double, 3 > TopInitCoordinates
Definition: fracture_propagation_2D_utilities.hpp:75
Kratos::FracturePropagation2DUtilities::Bifurcation::BotTipCoordinates
array_1d< double, 3 > BotTipCoordinates
Definition: fracture_propagation_2D_utilities.hpp:78
Kratos::FracturePropagation2DUtilities::FracturePoint
Structs for fracture propagation check --------------------------------------------------------------...
Definition: fracture_propagation_2D_utilities.hpp:54
Kratos::FracturePropagation2DUtilities::FracturePoint::TipDistance
double TipDistance
Definition: fracture_propagation_2D_utilities.hpp:56
Kratos::FracturePropagation2DUtilities::FracturePoint::Coordinates
array_1d< double, 2 > Coordinates
Definition: fracture_propagation_2D_utilities.hpp:55
Kratos::FracturePropagation2DUtilities::FracturePoint::Weight
double Weight
Definition: fracture_propagation_2D_utilities.hpp:56
Kratos::FracturePropagation2DUtilities::FracturePoint::Damage
double Damage
Definition: fracture_propagation_2D_utilities.hpp:56
Kratos::FracturePropagation2DUtilities::FracturePoint::pElement
Element::Pointer pElement
Definition: fracture_propagation_2D_utilities.hpp:57
Kratos::FracturePropagation2DUtilities::GaussPointOld
Structs for mapping model parts ---------------------------------------------------------------------...
Definition: fracture_propagation_2D_utilities.hpp:106
Kratos::FracturePropagation2DUtilities::GaussPointOld::Coordinates
array_1d< double, 3 > Coordinates
Definition: fracture_propagation_2D_utilities.hpp:107
Kratos::FracturePropagation2DUtilities::GaussPointOld::StateVariable
double StateVariable
Definition: fracture_propagation_2D_utilities.hpp:108
Kratos::FracturePropagation2DUtilities::GaussPointOld::Weight
double Weight
Definition: fracture_propagation_2D_utilities.hpp:108
Kratos::FracturePropagation2DUtilities::PropagationGlobalVariables
Definition: fracture_propagation_2D_utilities.hpp:84
Kratos::FracturePropagation2DUtilities::PropagationGlobalVariables::BifurcationVector
std::vector< Bifurcation > BifurcationVector
Definition: fracture_propagation_2D_utilities.hpp:89
Kratos::FracturePropagation2DUtilities::PropagationGlobalVariables::PropagationVector
std::vector< Propagation > PropagationVector
Definition: fracture_propagation_2D_utilities.hpp:88
Kratos::FracturePropagation2DUtilities::PropagationGlobalVariables::PropagateFractures
bool PropagateFractures
Definition: fracture_propagation_2D_utilities.hpp:87
Kratos::FracturePropagation2DUtilities::PropagationGlobalVariables::FracturePointsCellMatrix
std::vector< std::vector< std::vector< FracturePoint > > > FracturePointsCellMatrix
Definition: fracture_propagation_2D_utilities.hpp:85
Kratos::FracturePropagation2DUtilities::PropagationGlobalVariables::pProcessInfo
ProcessInfo::Pointer pProcessInfo
Definition: fracture_propagation_2D_utilities.hpp:86
Kratos::FracturePropagation2DUtilities::Propagation
Definition: fracture_propagation_2D_utilities.hpp:63
Kratos::FracturePropagation2DUtilities::Propagation::TopInitCoordinates
array_1d< double, 3 > TopInitCoordinates
Definition: fracture_propagation_2D_utilities.hpp:65
Kratos::FracturePropagation2DUtilities::Propagation::BotInitCoordinates
array_1d< double, 3 > BotInitCoordinates
Definition: fracture_propagation_2D_utilities.hpp:66
Kratos::FracturePropagation2DUtilities::Propagation::MotherFractureId
int MotherFractureId
Definition: fracture_propagation_2D_utilities.hpp:64
Kratos::FracturePropagation2DUtilities::Propagation::TipCoordinates
array_1d< double, 3 > TipCoordinates
Definition: fracture_propagation_2D_utilities.hpp:67
Kratos::FracturePropagation2DUtilities::PropagationLocalVariables
Definition: fracture_propagation_2D_utilities.hpp:95
Kratos::FracturePropagation2DUtilities::PropagationLocalVariables::TipCoordinates
array_1d< double, 2 > TipCoordinates
Definition: fracture_propagation_2D_utilities.hpp:97
Kratos::FracturePropagation2DUtilities::PropagationLocalVariables::RotationMatrix
BoundedMatrix< double, 2, 2 > RotationMatrix
Definition: fracture_propagation_2D_utilities.hpp:96
Kratos::FracturePropagation2DUtilities::PropagationLocalVariables::BotFrontFracturePoints
std::vector< FracturePoint * > BotFrontFracturePoints
Definition: fracture_propagation_2D_utilities.hpp:100
Kratos::FracturePropagation2DUtilities::PropagationLocalVariables::TopFrontFracturePoints
std::vector< FracturePoint * > TopFrontFracturePoints
Definition: fracture_propagation_2D_utilities.hpp:99
Kratos::FracturePropagation2DUtilities::PropagationLocalVariables::TipLocalCoordinates
array_1d< double, 2 > TipLocalCoordinates
Definition: fracture_propagation_2D_utilities.hpp:98
Kratos::FracturePropagation2DUtilities::UtilityVariables
Basic Structs for the utility -----------------------------------------------------------------------...
Definition: fracture_propagation_2D_utilities.hpp:45
Kratos::FracturePropagation2DUtilities::UtilityVariables::RowSize
double RowSize
Definition: fracture_propagation_2D_utilities.hpp:48
Kratos::FracturePropagation2DUtilities::UtilityVariables::Y_min
double Y_min
Definition: fracture_propagation_2D_utilities.hpp:46
Kratos::FracturePropagation2DUtilities::UtilityVariables::X_min
double X_min
Definition: fracture_propagation_2D_utilities.hpp:46
Kratos::FracturePropagation2DUtilities::UtilityVariables::NRows
int NRows
Definition: fracture_propagation_2D_utilities.hpp:47
Kratos::FracturePropagation2DUtilities::UtilityVariables::ColumnSize
double ColumnSize
Definition: fracture_propagation_2D_utilities.hpp:48
Kratos::FracturePropagation2DUtilities::UtilityVariables::X_max
double X_max
Definition: fracture_propagation_2D_utilities.hpp:46
Kratos::FracturePropagation2DUtilities::UtilityVariables::NColumns
int NColumns
Definition: fracture_propagation_2D_utilities.hpp:47
Kratos::FracturePropagation2DUtilities::UtilityVariables::Y_max
double Y_max
Definition: fracture_propagation_2D_utilities.hpp:46
Generated on Sat Feb 3 2024 02:11:11 for KratosMultiphysics by doxygen 1.9.1