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.
generate_new_nodes_before_meshing_process.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosPfemFluidDynamicsApplication $
3 // Created by: $Author: AFranci $
4 // Last modified by: $Co-Author: $
5 // Date: $Date: October 2016 $
6 // Revision: $Revision: 0.0 $
7 //
8 //
9 
10 #if !defined(KRATOS_GENERATE_NEW_NODES_BEFORE_MESHING_PROCESS_H_INCLUDED)
11 #define KRATOS_GENERATE_NEW_NODES_BEFORE_MESHING_PROCESS_H_INCLUDED
12 
13 // External includes
14 
15 // System includes
16 
17 // Project includes
18 #include "containers/variables_list_data_value_container.h"
19 #include "spatial_containers/spatial_containers.h"
20 
21 #include "includes/model_part.h"
22 #include "custom_utilities/mesh_error_calculation_utilities.hpp"
23 #include "custom_utilities/mesher_utilities.hpp"
24 #include "custom_processes/mesher_process.hpp"
25 
27 // Data:
28 // Flags: (checked)
29 // (set)
30 // (modified)
31 // (reset)
32 //(set):=(set in this process)
33 
34 namespace Kratos
35 {
36 
39 
41 
45  class GenerateNewNodesBeforeMeshingProcess
46  : public MesherProcess
47  {
48  public:
51 
53  KRATOS_CLASS_POINTER_DEFINITION(GenerateNewNodesBeforeMeshingProcess);
54 
55  typedef ModelPart::NodeType NodeType;
56  typedef ModelPart::ConditionType ConditionType;
57  typedef ModelPart::PropertiesType PropertiesType;
58  typedef ConditionType::GeometryType GeometryType;
59  typedef std::size_t SizeType;
60 
64 
66  GenerateNewNodesBeforeMeshingProcess(ModelPart &rModelPart,
67  MesherUtilities::MeshingParameters &rRemeshingParameters,
68  int EchoLevel)
69  : mrModelPart(rModelPart),
70  mrRemesh(rRemeshingParameters)
71  {
72  mEchoLevel = EchoLevel;
73  }
74 
76  virtual ~GenerateNewNodesBeforeMeshingProcess() {}
77 
81 
83  void operator()()
84  {
85  Execute();
86  }
87 
91 
93  void Execute() override
94  {
95  KRATOS_TRY
96 
97  if (mEchoLevel > 1)
98  std::cout << " [ GENERATE NEW NODES for homomgeneous mesh: " << std::endl;
99 
100  if (mrModelPart.Name() != mrRemesh.SubModelPartName)
101  std::cout << " ModelPart Supplied do not corresponds to the Meshing Domain: (" << mrModelPart.Name() << " != " << mrRemesh.SubModelPartName << ")" << std::endl;
102 
103  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
104  double currentTime = rCurrentProcessInfo[TIME];
105  double timeInterval = rCurrentProcessInfo[DELTA_TIME];
106  const SizeType dimension = mrModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
107 
108  bool refiningBox = false;
109  for (SizeType index = 0; index < mrRemesh.UseRefiningBox.size(); index++)
110  {
111  if (mrRemesh.UseRefiningBox[index] == true && currentTime > mrRemesh.RefiningBoxInitialTime[index] && currentTime < mrRemesh.RefiningBoxFinalTime[index])
112  {
113  refiningBox = true;
114  }
115  }
116 
117  if (currentTime < 2 * timeInterval)
118  {
119  mrRemesh.Info->RemovedNodes = 0;
120  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes = 0;
121 
122  if (mEchoLevel > 1)
123  std::cout << " First meshes: I repare the mesh without adding new nodes" << std::endl;
124  mrRemesh.Info->InitialNumberOfNodes = mrRemesh.Info->NumberOfNodes;
125  }
126 
127  int ElementsToRefine = mrRemesh.Info->RemovedNodes;
128  SizeType eulerianInletNodes = mrRemesh.Info->NumberOfEulerianInletNodes;
129 
130  int initialNumberOfNodes = mrRemesh.Info->InitialNumberOfNodes;
131  int numberOfNodes = mrRemesh.Info->NumberOfNodes;
132  int extraNodes = initialNumberOfNodes - numberOfNodes;
133  int toleredExtraNodes = int(0.05 * mrRemesh.Info->InitialNumberOfNodes);
134 
135  if (mrRemesh.ExecutionOptions.Is(MesherUtilities::REFINE_WALL_CORNER))
136  {
137  if ((ElementsToRefine - extraNodes) > toleredExtraNodes && refiningBox == false)
138  {
139  ElementsToRefine = toleredExtraNodes + extraNodes;
140  if (ElementsToRefine < 0)
141  {
142  ElementsToRefine = 0;
143  }
144  }
145  }
146 
147  if (ElementsToRefine > 0 && mEchoLevel > 1)
148  std::cout << " I will look for " << ElementsToRefine << " new nodes" << std::endl;
149 
150  if (refiningBox == false)
151  {
152 
153  if (ElementsToRefine > 0 || eulerianInletNodes > 0)
154  {
155  std::vector<array_1d<double, 3>> new_positions;
156  std::vector<double> biggest_volumes;
157  std::vector<array_1d<SizeType, 4>> nodes_id_to_interpolate;
158 
159  int CountNodes = 0;
160 
161  new_positions.resize(ElementsToRefine);
162  biggest_volumes.resize(ElementsToRefine, false);
163  nodes_id_to_interpolate.resize(ElementsToRefine);
164 
165  for (int nn = 0; nn < ElementsToRefine; nn++)
166  {
167  biggest_volumes[nn] = -1.0;
168  }
169 
170  // std::vector<array_1d<double, 3>> CornerWallNewPositions;
171  // std::vector<array_1d<SizeType, 4>> CornerWallNodesIDToInterpolate;
172  // std::vector<Node::DofsContainerType> CornerWallNewDofs;
173  // int cornerWallNewNodes = 0;
174  // int maxOfNewWallNodes = toleredExtraNodes;
175  // if (mrRemesh.ExecutionOptions.Is(MesherUtilities::REFINE_WALL_CORNER))
176  // {
177  // CornerWallNewPositions.resize(maxOfNewWallNodes, false);
178  // CornerWallNodesIDToInterpolate.resize(maxOfNewWallNodes, false);
179  // CornerWallNewDofs.resize(maxOfNewWallNodes, false);
180  // }
181  SizeType addedNodesAtEulerianInlet = 0;
182  ModelPart::ElementsContainerType::iterator element_begin = mrModelPart.ElementsBegin();
183  for (ModelPart::ElementsContainerType::const_iterator ie = element_begin; ie != mrModelPart.ElementsEnd(); ie++)
184  {
185 
187  if (dimension == 2)
188  {
189  SelectEdgeToRefine2D(ie->GetGeometry(), new_positions, biggest_volumes, nodes_id_to_interpolate, CountNodes, ElementsToRefine, addedNodesAtEulerianInlet);
190  }
191  else if (dimension == 3)
192  {
193  SelectEdgeToRefine3D(ie->GetGeometry(), new_positions, biggest_volumes, nodes_id_to_interpolate, CountNodes, ElementsToRefine, addedNodesAtEulerianInlet);
194  }
195 
196  } // elements loop
197  ElementsToRefine += addedNodesAtEulerianInlet;
198 
199  mrRemesh.Info->RemovedNodes -= ElementsToRefine;
200  if (CountNodes < ElementsToRefine)
201  {
202  mrRemesh.Info->RemovedNodes += ElementsToRefine - CountNodes;
203  new_positions.resize(CountNodes);
204  biggest_volumes.resize(CountNodes, false);
205  nodes_id_to_interpolate.resize(CountNodes);
206  }
207  SizeType maxId = 0;
208  CreateAndAddNewNodes(new_positions, nodes_id_to_interpolate, ElementsToRefine, maxId);
209 
210  // if (mrRemesh.ExecutionOptions.Is(MesherUtilities::REFINE_WALL_CORNER))
211  // {
212  // if (cornerWallNewNodes < maxOfNewWallNodes)
213  // {
214  // CornerWallNewPositions.resize(cornerWallNewNodes, false);
215  // CornerWallNewDofs.resize(cornerWallNewNodes, false);
216  // CornerWallNodesIDToInterpolate.resize(cornerWallNewNodes, false);
217  // }
218  // CreateAndAddNewNodesInCornerWall(CornerWallNewPositions, CornerWallNodesIDToInterpolate, CornerWallNewDofs, cornerWallNewNodes, maxId);
219  // }
220  }
221  }
222  else
223  {
224 
225  std::vector<array_1d<double, 3>> new_positions;
226  std::vector<array_1d<SizeType, 4>> nodes_id_to_interpolate;
227 
228  int CountNodes = 0;
229 
230  new_positions.resize(0);
231  nodes_id_to_interpolate.resize(0);
232 
233  ModelPart::ElementsContainerType::iterator element_begin = mrModelPart.ElementsBegin();
234  int nodesInTransitionZone = 0;
235  for (ModelPart::ElementsContainerType::const_iterator ie = element_begin; ie != mrModelPart.ElementsEnd(); ie++)
236  {
237  const SizeType dimension = ie->GetGeometry().WorkingSpaceDimension();
239  if (dimension == 2)
240  {
241  SelectEdgeToRefine2DWithRefinement(ie->GetGeometry(), new_positions, nodes_id_to_interpolate, CountNodes, ElementsToRefine, nodesInTransitionZone);
242  }
243  else if (dimension == 3)
244  {
245  SelectEdgeToRefine3DWithRefinement(ie->GetGeometry(), new_positions, nodes_id_to_interpolate, CountNodes, ElementsToRefine, nodesInTransitionZone);
246  }
247 
248  } // elements loop
249 
250  mrRemesh.Info->RemovedNodes -= CountNodes;
251  if (CountNodes < ElementsToRefine)
252  {
253  mrRemesh.Info->RemovedNodes += ElementsToRefine - CountNodes;
254  new_positions.resize(CountNodes);
255  nodes_id_to_interpolate.resize(CountNodes);
256  }
257  SizeType maxId = 0;
258  CreateAndAddNewNodes(new_positions, nodes_id_to_interpolate, CountNodes, maxId);
259  }
260 
261  mrRemesh.InputInitializedFlag = false;
262 
263  if (mEchoLevel > 1)
264  std::cout << " GENERATE NEW NODES ]; " << std::endl;
265 
266  KRATOS_CATCH(" ")
267  }
268 
272 
276 
280 
282  std::string Info() const override
283  {
284  return "GenerateNewNodesBeforeMeshingProcess";
285  }
286 
288  void PrintInfo(std::ostream &rOStream) const override
289  {
290  rOStream << "GenerateNewNodesBeforeMeshingProcess";
291  }
292 
294  void PrintData(std::ostream &rOStream) const override
295  {
296  }
297 
301 
303 
304  private:
307 
311  ModelPart &mrModelPart;
312 
313  MesherUtilities::MeshingParameters &mrRemesh;
314 
315  MesherUtilities mMesherUtilities;
316 
317  int mEchoLevel;
318 
322 
326 
327  void CopyDofs(Node::DofsContainerType const &From, Node::DofsContainerType &To)
328  {
329  for (auto &p_dof : From)
330  {
331  To.push_back(Kratos::unique_ptr<Dof<double>>(new Dof<double>(*p_dof)));
332  }
333  }
334 
335  void InsertNodeInCornerElement2D(Element::GeometryType &Element,
336  std::vector<array_1d<double, 3>> &new_positions,
337  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
338  std::vector<Node::DofsContainerType> &NewDofs,
339  int &CountNodes)
340  {
341  KRATOS_TRY
342 
343  const SizeType nds = Element.size();
344 
345  SizeType rigidNodes = 0;
346  SizeType freesurfaceNodes = 0;
347 
348  for (SizeType pn = 0; pn < nds; ++pn)
349  {
350  if (Element[pn].Is(RIGID))
351  {
352  rigidNodes++;
353  }
354  if (Element[pn].Is(FREE_SURFACE))
355  {
356  freesurfaceNodes++;
357  }
358  }
359  double cosTolerance = 0.01;
360 
361  if (rigidNodes == 2 && freesurfaceNodes == 0)
362  {
363  array_1d<double, 2> NormalA(2, 0.0);
364  array_1d<double, 2> NormalB(2, 0.0);
365  double cosAngle = 1;
366 
367  if ((Element[0].Is(RIGID) && Element[1].Is(RIGID)) || (Element[0].Is(INLET) && Element[1].Is(INLET)))
368  {
369  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
370  NormalB = Element[1].FastGetSolutionStepValue(NORMAL);
371  cosAngle = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1];
372  if (cosAngle < cosTolerance && cosAngle > -cosTolerance)
373  {
374  array_1d<double, 3> new_position = (Element[0].Coordinates() + Element[1].Coordinates()) * 0.5;
375  nodes_id_to_interpolate[CountNodes][0] = Element[0].GetId();
376  nodes_id_to_interpolate[CountNodes][1] = Element[1].GetId();
377  if (Element[2].IsNot(TO_ERASE))
378  {
379  nodes_id_to_interpolate[CountNodes][2] = Element[2].GetId();
380  }
381  else
382  {
383  nodes_id_to_interpolate[CountNodes][2] = Element[0].GetId();
384  }
385  CopyDofs(Element[2].GetDofs(), NewDofs[CountNodes]);
386  new_positions[CountNodes] = new_position;
387  CountNodes++;
388  }
389  }
390  else if (Element[0].Is(RIGID) && Element[2].Is(RIGID))
391  {
392  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
393  NormalB = Element[2].FastGetSolutionStepValue(NORMAL);
394  cosAngle = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1];
395  if (cosAngle < cosTolerance && cosAngle > -cosTolerance)
396  {
397  array_1d<double, 3> new_position = (Element[0].Coordinates() + Element[1].Coordinates()) * 0.5;
398  nodes_id_to_interpolate[CountNodes][0] = Element[0].GetId();
399  nodes_id_to_interpolate[CountNodes][1] = Element[2].GetId();
400  if (Element[1].IsNot(TO_ERASE))
401  {
402  nodes_id_to_interpolate[CountNodes][2] = Element[1].GetId();
403  }
404  else
405  {
406  nodes_id_to_interpolate[CountNodes][2] = Element[0].GetId();
407  }
408  CopyDofs(Element[1].GetDofs(), NewDofs[CountNodes]);
409  new_positions[CountNodes] = new_position;
410  CountNodes++;
411  }
412  }
413  else if (Element[1].Is(RIGID) && Element[2].Is(RIGID))
414  {
415  NormalA = Element[1].FastGetSolutionStepValue(NORMAL);
416  NormalB = Element[2].FastGetSolutionStepValue(NORMAL);
417  cosAngle = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1];
418  if (cosAngle < cosTolerance && cosAngle > -cosTolerance)
419  {
420  array_1d<double, 3> new_position = (Element[2].Coordinates() + Element[1].Coordinates()) * 0.5;
421  nodes_id_to_interpolate[CountNodes][0] = Element[2].GetId();
422  nodes_id_to_interpolate[CountNodes][1] = Element[1].GetId();
423  if (Element[0].IsNot(TO_ERASE))
424  {
425  nodes_id_to_interpolate[CountNodes][2] = Element[0].GetId();
426  }
427  else
428  {
429  nodes_id_to_interpolate[CountNodes][2] = Element[2].GetId();
430  }
431  CopyDofs(Element[0].GetDofs(), NewDofs[CountNodes]);
432  new_positions[CountNodes] = new_position;
433  CountNodes++;
434  }
435  }
436  }
437 
438  KRATOS_CATCH("")
439  }
440 
441  void InsertNodeInCornerElement3D(Element::GeometryType &Element,
442  std::vector<array_1d<double, 3>> &new_positions,
443  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
444  std::vector<Node::DofsContainerType> &NewDofs,
445  int &CountNodes)
446  {
447  KRATOS_TRY
448 
449  const SizeType nds = Element.size();
450 
451  SizeType rigidNodes = 0;
452  SizeType freesurfaceNodes = 0;
453  SizeType toEraseNodes = 0;
454 
455  for (SizeType pn = 0; pn < nds; ++pn)
456  {
457  if (Element[pn].Is(RIGID))
458  {
459  rigidNodes++;
460  }
461  if (Element[pn].Is(FREE_SURFACE))
462  {
463  freesurfaceNodes++;
464  }
465  if (Element[pn].Is(TO_ERASE))
466  {
467  toEraseNodes++;
468  }
469  }
470 
471  if (rigidNodes == 3 && freesurfaceNodes == 0 && toEraseNodes == 0)
472  {
473  array_1d<double, 3> NormalA(3, 0.0);
474  array_1d<double, 3> NormalB(3, 0.0);
475  double normNormalA = 0;
476  double normNormalB = 0;
477  double cos = 1.0;
478  double minCos = 1.0;
479  array_1d<SizeType, 2> idsWallNodes(2, 0);
480  SizeType idFreeNode = 0;
481  double cosTolerance = 0.1;
482  if (Element[0].IsNot(RIGID))
483  {
484  NormalA = Element[1].FastGetSolutionStepValue(NORMAL);
485  NormalB = Element[2].FastGetSolutionStepValue(NORMAL);
486  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
487  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
488  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
489  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
490  {
491  minCos = cos;
492  idsWallNodes[0] = 1;
493  idsWallNodes[1] = 2;
494  idFreeNode = 0;
495  }
496 
497  NormalA = Element[1].FastGetSolutionStepValue(NORMAL);
498  NormalB = Element[3].FastGetSolutionStepValue(NORMAL);
499  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
500  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
501  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
502  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
503  {
504  minCos = cos;
505  idsWallNodes[0] = 1;
506  idsWallNodes[1] = 3;
507  idFreeNode = 0;
508  }
509 
510  NormalA = Element[2].FastGetSolutionStepValue(NORMAL);
511  NormalB = Element[3].FastGetSolutionStepValue(NORMAL);
512  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
513  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
514  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
515  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
516  {
517  minCos = cos;
518  idsWallNodes[0] = 2;
519  idsWallNodes[1] = 3;
520  idFreeNode = 0;
521  }
522  }
523  else if (Element[1].IsNot(RIGID))
524  {
525  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
526  NormalB = Element[2].FastGetSolutionStepValue(NORMAL);
527  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
528  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
529  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
530  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
531  {
532  minCos = cos;
533  idsWallNodes[0] = 0;
534  idsWallNodes[1] = 2;
535  idFreeNode = 1;
536  }
537 
538  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
539  NormalB = Element[3].FastGetSolutionStepValue(NORMAL);
540  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
541  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
542  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
543  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
544  {
545  minCos = cos;
546  idsWallNodes[0] = 0;
547  idsWallNodes[1] = 3;
548  idFreeNode = 1;
549  }
550 
551  NormalA = Element[2].FastGetSolutionStepValue(NORMAL);
552  NormalB = Element[3].FastGetSolutionStepValue(NORMAL);
553  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
554  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
555  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
556  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
557  {
558  minCos = cos;
559  idsWallNodes[0] = 2;
560  idsWallNodes[1] = 3;
561  idFreeNode = 1;
562  }
563  }
564  else if (Element[2].IsNot(RIGID))
565  {
566 
567  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
568  NormalB = Element[1].FastGetSolutionStepValue(NORMAL);
569  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
570  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
571  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
572  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
573  {
574  minCos = cos;
575  idsWallNodes[0] = 0;
576  idsWallNodes[1] = 1;
577  idFreeNode = 2;
578  }
579 
580  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
581  NormalB = Element[3].FastGetSolutionStepValue(NORMAL);
582  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
583  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
584  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
585  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
586  {
587  minCos = cos;
588  idsWallNodes[0] = 0;
589  idsWallNodes[1] = 3;
590  idFreeNode = 2;
591  }
592 
593  NormalA = Element[1].FastGetSolutionStepValue(NORMAL);
594  NormalB = Element[3].FastGetSolutionStepValue(NORMAL);
595  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
596  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
597  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
598  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
599  {
600  minCos = cos;
601  idsWallNodes[0] = 1;
602  idsWallNodes[1] = 3;
603  idFreeNode = 2;
604  }
605  }
606  else if (Element[3].IsNot(RIGID))
607  {
608 
609  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
610  NormalB = Element[1].FastGetSolutionStepValue(NORMAL);
611  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
612  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
613  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
614  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
615  {
616  minCos = cos;
617  idsWallNodes[0] = 0;
618  idsWallNodes[1] = 1;
619  idFreeNode = 3;
620  }
621 
622  NormalA = Element[0].FastGetSolutionStepValue(NORMAL);
623  NormalB = Element[2].FastGetSolutionStepValue(NORMAL);
624  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
625  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
626  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
627  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
628  {
629  minCos = cos;
630  idsWallNodes[0] = 0;
631  idsWallNodes[1] = 2;
632  idFreeNode = 3;
633  }
634 
635  NormalA = Element[1].FastGetSolutionStepValue(NORMAL);
636  NormalB = Element[2].FastGetSolutionStepValue(NORMAL);
637  normNormalA = NormalA[0] * NormalA[0] + NormalA[1] * NormalA[1] + NormalA[2] * NormalA[2];
638  normNormalB = NormalB[0] * NormalB[0] + NormalB[1] * NormalB[1] + NormalB[2] * NormalB[2];
639  cos = NormalA[0] * NormalB[0] + NormalA[1] * NormalB[1] + NormalA[2] * NormalB[2];
640  if (cos < minCos && (cos < cosTolerance && cos > -cosTolerance) && (normNormalA > 0.99 && normNormalA < 1.01) && (normNormalB > 0.99 && normNormalB < 1.01))
641  {
642  minCos = cos;
643  idsWallNodes[0] = 1;
644  idsWallNodes[1] = 2;
645  idFreeNode = 3;
646  }
647  }
648 
649  if (minCos < cosTolerance && minCos > -cosTolerance)
650  {
651 
652  bool alreadyAddedNode = false;
653  SizeType idA = Element[idsWallNodes[0]].GetId();
654  SizeType idB = Element[idsWallNodes[1]].GetId();
655  double minimumDistanceToInstert = 1.3 * mrRemesh.Refine->CriticalRadius;
656  array_1d<double, 3> CoorDifference = Element[idsWallNodes[0]].Coordinates() - Element[idsWallNodes[1]].Coordinates();
657  double SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
658  double separation = std::sqrt(SquaredLength);
659  SizeType idC = Element[idFreeNode].GetId();
660  if (separation > minimumDistanceToInstert)
661  {
662 
663  for (SizeType i = 0; i < unsigned(CountNodes); i++)
664  {
665  if (idA == nodes_id_to_interpolate[i][0] || idA == nodes_id_to_interpolate[i][1] || idB == nodes_id_to_interpolate[i][0] || idB == nodes_id_to_interpolate[i][1])
666  {
667  alreadyAddedNode = true;
668  break;
669  }
670  }
671  if (alreadyAddedNode == false)
672  {
673  array_1d<double, 3> new_position = (Element[idsWallNodes[0]].Coordinates() + Element[idsWallNodes[1]].Coordinates()) * 0.5;
674  nodes_id_to_interpolate[CountNodes][0] = idA;
675  nodes_id_to_interpolate[CountNodes][1] = idB;
676  if (Element[idFreeNode].IsNot(TO_ERASE))
677  {
678  nodes_id_to_interpolate[CountNodes][2] = idC;
679  }
680  else
681  {
682  nodes_id_to_interpolate[CountNodes][2] = idA;
683  }
684  CopyDofs(Element[idFreeNode].GetDofs(), NewDofs[CountNodes]);
685  new_positions[CountNodes] = new_position;
686  CountNodes++;
687  }
688  }
689  }
690  }
691 
692  KRATOS_CATCH("")
693  }
694 
695  void SelectEdgeToRefine2D(Element::GeometryType &Element,
696  std::vector<array_1d<double, 3>> &new_positions,
697  std::vector<double> &biggest_volumes,
698  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
699  int &CountNodes,
700  int &ElementsToRefine,
701  SizeType &addedNodesAtEulerianInlet)
702  {
703  KRATOS_TRY
704 
705  const SizeType nds = Element.size();
706  double meanMeshSize = mrRemesh.Refine->CriticalRadius;
707 
708  SizeType rigidNodes = 0;
709  SizeType boundaryNodes = 0;
710  SizeType freesurfaceNodes = 0;
711  SizeType lagrangianInletNodes = 0;
712  SizeType eulerianInletNodes = 0;
713  bool toEraseNodeFound = false;
714  double rigidNodeLocalMeshSize = 0;
715  double rigidNodeMeshCounter = 0;
716  bool suitableElementForSecondAdd = true;
717 
718  for (SizeType pn = 0; pn < nds; ++pn)
719  {
720  if (Element[pn].Is(RIGID))
721  {
722  rigidNodes++;
723  rigidNodeLocalMeshSize += Element[pn].FastGetSolutionStepValue(NODAL_H_WALL);
724  rigidNodeMeshCounter += 1.0;
725  }
726  if (Element[pn].Is(BOUNDARY))
727  {
728  boundaryNodes++;
729  }
730  if (Element[pn].Is(TO_ERASE))
731  {
732  toEraseNodeFound = true;
733  }
734  if (Element[pn].Is(FREE_SURFACE))
735  {
736  freesurfaceNodes++;
737  }
738  if (Element[pn].Is(PFEMFlags::LAGRANGIAN_INLET))
739  {
740  lagrangianInletNodes++;
741  }
742  if (Element[pn].Is(PFEMFlags::EULERIAN_INLET))
743  {
744  eulerianInletNodes++;
745  }
746  }
747 
748  if (rigidNodeMeshCounter > 0)
749  {
750  const double rigidWallMeshSize = rigidNodeLocalMeshSize / rigidNodeMeshCounter;
751  const double ratio = rigidWallMeshSize / meanMeshSize;
752  const double tolerance = 1.8;
753  if (ratio > tolerance)
754  {
755  meanMeshSize *= 0.5;
756  meanMeshSize += 0.5 * rigidWallMeshSize;
757  }
758  }
759 
760  const double limitEdgeLength = 1.4 * meanMeshSize;
761  const double safetyCoefficient2D = 1.5;
762  double penalization = 1.0; // to penalize adding node, penalization here should be smaller than 1
763  if (rigidNodes > 1)
764  {
765  penalization = 0.8;
766  if (lagrangianInletNodes > 0)
767  {
768  penalization = 0.9;
769  }
770  }
771  else if (rigidNodes > 0 && freesurfaceNodes > 0 && eulerianInletNodes == 0)
772  {
773  penalization = 0;
774  }
775  else if (freesurfaceNodes > 0)
776  {
777  penalization = 0.875;
778  }
779 
780  double ElementalVolume = Element.Area();
781 
782  array_1d<double, 3> Edges(3, 0.0);
783  array_1d<SizeType, 3> FirstEdgeNode(3, 0);
784  array_1d<SizeType, 3> SecondEdgeNode(3, 0);
785  double WallCharacteristicDistance = 0;
786  array_1d<double, 3> CoorDifference = Element[1].Coordinates() - Element[0].Coordinates();
787  double SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
788  Edges[0] = std::sqrt(SquaredLength);
789  FirstEdgeNode[0] = 0;
790  SecondEdgeNode[0] = 1;
791  if ((Element[0].Is(RIGID) && Element[1].Is(RIGID)) || (Element[0].Is(INLET) && Element[1].Is(INLET)))
792  {
793  WallCharacteristicDistance = Edges[0];
794  }
795  SizeType Counter = 0;
796  for (SizeType i = 2; i < nds; i++)
797  {
798  for (SizeType j = 0; j < i; j++)
799  {
800  noalias(CoorDifference) = Element[i].Coordinates() - Element[j].Coordinates();
801  SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
802  Counter += 1;
803  Edges[Counter] = std::sqrt(SquaredLength);
804  FirstEdgeNode[Counter] = j;
805  SecondEdgeNode[Counter] = i;
806  if (((Element[i].Is(RIGID) && Element[j].Is(RIGID)) || (Element[i].Is(INLET) && Element[j].Is(INLET))) && Edges[Counter] > WallCharacteristicDistance)
807  {
808  WallCharacteristicDistance = Edges[Counter];
809  }
810  }
811  }
812 
813  bool dangerousElement = false;
814 
815  for (SizeType i = 0; i < 3; i++)
816  {
817  if (rigidNodes > 1)
818  {
819  if ((Edges[i] < WallCharacteristicDistance * safetyCoefficient2D && (Element[FirstEdgeNode[i]].Is(RIGID) || Element[SecondEdgeNode[i]].Is(RIGID))) ||
820  ((Element[FirstEdgeNode[i]].Is(RIGID) && Element[SecondEdgeNode[i]].Is(RIGID)) || (Element[FirstEdgeNode[i]].Is(INLET) && Element[SecondEdgeNode[i]].Is(INLET))))
821  {
822  Edges[i] = 0;
823  }
824  if ((Element[FirstEdgeNode[i]].Is(FREE_SURFACE) || Element[FirstEdgeNode[i]].Is(RIGID)) &&
825  (Element[SecondEdgeNode[i]].Is(FREE_SURFACE) || Element[SecondEdgeNode[i]].Is(RIGID)) &&
826  (Element[FirstEdgeNode[i]].IsNot(PFEMFlags::EULERIAN_INLET) && Element[SecondEdgeNode[i]].IsNot(PFEMFlags::EULERIAN_INLET)))
827  {
828  Edges[i] = 0;
829  }
830  }
831  else if (rigidNodes == 0)
832  {
833  SizeType propertyIdFirstNode = Element[FirstEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
834  SizeType propertyIdSecondNode = Element[SecondEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
835  if (propertyIdFirstNode != propertyIdSecondNode)
836  {
837  penalization = 0.9; // 10% less than normal nodes
838  }
839  }
840  }
841  if ((Edges[0] == 0 && Edges[1] == 0 && Edges[2] == 0) || rigidNodes == 3)
842  {
843  dangerousElement = true;
844  }
845 
846  if (dangerousElement == false && toEraseNodeFound == false)
847  {
848  SizeType maxCount = 3;
849  double LargestEdge = 0;
850 
851  for (SizeType i = 0; i < 3; i++)
852  {
853  if (Edges[i] > LargestEdge)
854  {
855  maxCount = i;
856  LargestEdge = Edges[i];
857  }
858  }
859 
860  if (CountNodes < ElementsToRefine && LargestEdge > limitEdgeLength && eulerianInletNodes == 0)
861  {
862 
863  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
864 
865  bool suitableElement = true;
866  for (int j = 0; j < CountNodes; j++)
867  {
868  const double diffX = std::abs(new_positions[j][0] - new_position[0]) - meanMeshSize * 0.5;
869  const double diffY = std::abs(new_positions[j][1] - new_position[1]) - meanMeshSize * 0.5;
870  if (diffX < 0 && diffY < 0) // the node is in the same zone of a previously inserted node
871  {
872  suitableElement = false;
873  }
874  }
875 
876  if (suitableElement)
877  {
878  nodes_id_to_interpolate[CountNodes][0] = Element[FirstEdgeNode[maxCount]].GetId();
879  nodes_id_to_interpolate[CountNodes][1] = Element[SecondEdgeNode[maxCount]].GetId();
880  biggest_volumes[CountNodes] = ElementalVolume;
881  new_positions[CountNodes] = new_position;
882  CountNodes++;
883  }
884  }
885  else if (freesurfaceNodes < 3 && rigidNodes < 3 && eulerianInletNodes == 0)
886  {
887  ElementalVolume *= penalization;
888  for (int nn = 0; nn < ElementsToRefine; nn++)
889  {
890  if (ElementalVolume > biggest_volumes[nn])
891  {
892 
893  if (maxCount < 3 && LargestEdge > limitEdgeLength)
894  {
895  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
896  if (ElementsToRefine > 1 && CountNodes > 0)
897  {
898  for (int j = 0; j < ElementsToRefine; j++)
899  {
900  const double diffX = std::abs(new_positions[j][0] - new_position[0]) - meanMeshSize * 0.5;
901  const double diffY = std::abs(new_positions[j][1] - new_position[1]) - meanMeshSize * 0.5;
902  if (diffX < 0 && diffY < 0) // the node is in the same zone of a previously inserted node
903  {
904  suitableElementForSecondAdd = false;
905  }
906  }
907  }
908 
909  if (suitableElementForSecondAdd)
910  {
911  nodes_id_to_interpolate[nn][0] = Element[FirstEdgeNode[maxCount]].GetId();
912  nodes_id_to_interpolate[nn][1] = Element[SecondEdgeNode[maxCount]].GetId();
913  biggest_volumes[nn] = ElementalVolume;
914  new_positions[nn] = new_position;
915  }
916  }
917 
918  break;
919  }
920  }
921  }
922 
923  if (eulerianInletNodes > 0 && LargestEdge > (2.0 * meanMeshSize))
924  {
925 
926  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
927 
928  bool suitableElement = true;
929  for (int j = 0; j < CountNodes; j++)
930  {
931  const double diffX = std::abs(new_positions[j][0] - new_position[0]) - meanMeshSize * 0.5;
932  const double diffY = std::abs(new_positions[j][1] - new_position[1]) - meanMeshSize * 0.5;
933  // if (diffX < 0 && diffY < 0) // the node is in the same zone of a previously inserted node
934  if ((diffX < 0 && diffY < 0) || (Element[FirstEdgeNode[maxCount]].IsNot(INLET) && Element[SecondEdgeNode[maxCount]].IsNot(INLET))) // the node is in the same zone of a previously inserted node
935 
936  {
937  suitableElement = false;
938  }
939  }
940 
941  if (suitableElement)
942  {
943  if (CountNodes >= ElementsToRefine)
944  {
945  new_positions.resize(CountNodes + 1);
946  biggest_volumes.resize(CountNodes + 1, false);
947  nodes_id_to_interpolate.resize(CountNodes + 1);
948  addedNodesAtEulerianInlet++;
949  }
950  nodes_id_to_interpolate[CountNodes][0] = Element[FirstEdgeNode[maxCount]].GetId();
951  nodes_id_to_interpolate[CountNodes][1] = Element[SecondEdgeNode[maxCount]].GetId();
952  biggest_volumes[CountNodes] = ElementalVolume;
953  new_positions[CountNodes] = new_position;
954  CountNodes++;
955  }
956  }
957  }
958 
959  KRATOS_CATCH("")
960  }
961 
962  void SelectEdgeToRefine3D(Element::GeometryType &Element,
963  std::vector<array_1d<double, 3>> &new_positions,
964  std::vector<double> &biggest_volumes,
965  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
966  int &CountNodes,
967  int &ElementsToRefine,
968  SizeType &addedNodesAtEulerianInlet)
969  {
970  KRATOS_TRY
971 
972  const SizeType nds = Element.size();
973  double meanMeshSize = mrRemesh.Refine->CriticalRadius;
974 
975  SizeType rigidNodes = 0;
976  SizeType freesurfaceNodes = 0;
977  SizeType lagrangianInletNodes = 0;
978  SizeType eulerianInletNodes = 0;
979  bool toEraseNodeFound = false;
980  double rigidNodeLocalMeshSize = 0;
981  double rigidNodeMeshCounter = 0;
982  bool suitableElementForSecondAdd = true;
983 
984  for (SizeType pn = 0; pn < nds; ++pn)
985  {
986  if (Element[pn].Is(RIGID))
987  {
988  rigidNodes++;
989  rigidNodeLocalMeshSize += Element[pn].FastGetSolutionStepValue(NODAL_H_WALL);
990  rigidNodeMeshCounter += 1.0;
991  }
992  if (Element[pn].Is(TO_ERASE))
993  {
994  toEraseNodeFound = true;
995  }
996  if (Element[pn].Is(FREE_SURFACE))
997  {
998  freesurfaceNodes++;
999  }
1000  if (Element[pn].Is(PFEMFlags::LAGRANGIAN_INLET))
1001  {
1002  lagrangianInletNodes++;
1003  }
1004  if (Element[pn].Is(PFEMFlags::EULERIAN_INLET))
1005  {
1006  eulerianInletNodes++;
1007  }
1008  }
1009 
1010  if (rigidNodeMeshCounter > 0)
1011  {
1012  const double rigidWallMeshSize = rigidNodeLocalMeshSize / rigidNodeMeshCounter;
1013  const double ratio = rigidWallMeshSize / meanMeshSize;
1014  const double tolerance = 1.8;
1015  if (ratio > tolerance)
1016  {
1017  meanMeshSize *= 0.5;
1018  meanMeshSize += 0.5 * rigidWallMeshSize;
1019  }
1020  }
1021 
1022  const double limitEdgeLength = 1.25 * meanMeshSize;
1023  const double safetyCoefficient3D = 1.6;
1024  double penalization = 1.0; // penalization here should be smaller than 1
1025  if (rigidNodes > 2)
1026  {
1027  penalization = 0.7;
1028  if (lagrangianInletNodes > 0)
1029  {
1030  penalization = 0.9;
1031  }
1032  }
1033  else if (rigidNodes > 0 && freesurfaceNodes > 0 && eulerianInletNodes == 0)
1034  {
1035  penalization = 0;
1036  }
1037  else if (freesurfaceNodes > 0)
1038  {
1039  penalization = 0.95;
1040  }
1041 
1042  double ElementalVolume = Element.Volume();
1043 
1044  array_1d<double, 6> Edges(6, 0.0);
1045  array_1d<SizeType, 6> FirstEdgeNode(6, 0);
1046  array_1d<SizeType, 6> SecondEdgeNode(6, 0);
1047  double WallCharacteristicDistance = 0;
1048  array_1d<double, 3> CoorDifference = Element[1].Coordinates() - Element[0].Coordinates();
1049  double SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1] + CoorDifference[2] * CoorDifference[2];
1050  Edges[0] = std::sqrt(SquaredLength);
1051  FirstEdgeNode[0] = 0;
1052  SecondEdgeNode[0] = 1;
1053  if ((Element[0].Is(RIGID) && Element[1].Is(RIGID)) || (Element[0].Is(INLET) && Element[1].Is(INLET)))
1054  {
1055  WallCharacteristicDistance = Edges[0];
1056  }
1057  SizeType Counter = 0;
1058  for (SizeType i = 2; i < nds; i++)
1059  {
1060  for (SizeType j = 0; j < i; j++)
1061  {
1062  noalias(CoorDifference) = Element[i].Coordinates() - Element[j].Coordinates();
1063  SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1] + CoorDifference[2] * CoorDifference[2];
1064  Counter += 1;
1065  Edges[Counter] = std::sqrt(SquaredLength);
1066  FirstEdgeNode[Counter] = j;
1067  SecondEdgeNode[Counter] = i;
1068  if (((Element[i].Is(RIGID) && Element[j].Is(RIGID)) || (Element[i].Is(INLET) && Element[j].Is(INLET))) && Edges[Counter] > WallCharacteristicDistance)
1069  {
1070  WallCharacteristicDistance = Edges[Counter];
1071  }
1072  }
1073  }
1074  // Edges connectivity: Edges[0]=d01, Edges[1]=d20, Edges[2]=d21, Edges[3]=d30, Edges[4]=d31, Edges[5]=d32
1075  bool dangerousElement = false;
1076 
1077  for (SizeType i = 0; i < 6; i++)
1078  {
1079  if (rigidNodes > 1)
1080  {
1081  if ((Edges[i] < WallCharacteristicDistance * safetyCoefficient3D && (Element[FirstEdgeNode[i]].Is(RIGID) || Element[SecondEdgeNode[i]].Is(RIGID))) ||
1082  ((Element[FirstEdgeNode[i]].Is(RIGID) && Element[SecondEdgeNode[i]].Is(RIGID)) || (Element[FirstEdgeNode[i]].Is(INLET) && Element[SecondEdgeNode[i]].Is(INLET))))
1083  {
1084  Edges[i] = 0;
1085  }
1086  if ((Element[FirstEdgeNode[i]].Is(FREE_SURFACE) || Element[FirstEdgeNode[i]].Is(RIGID)) &&
1087  (Element[SecondEdgeNode[i]].Is(FREE_SURFACE) || Element[SecondEdgeNode[i]].Is(RIGID)) &&
1088  (Element[FirstEdgeNode[i]].IsNot(PFEMFlags::EULERIAN_INLET) && Element[SecondEdgeNode[i]].IsNot(PFEMFlags::EULERIAN_INLET)))
1089  {
1090  Edges[i] = 0;
1091  }
1092  }
1093  else if (rigidNodes == 0)
1094  {
1095  SizeType propertyIdFirstNode = Element[FirstEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1096  SizeType propertyIdSecondNode = Element[SecondEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1097  if (propertyIdFirstNode != propertyIdSecondNode)
1098  {
1099  penalization = 0.8; // 20% less than normal nodes
1100  }
1101  }
1102  }
1103  if (rigidNodes == 1)
1104  {
1105  if (Element[0].Is(RIGID))
1106  {
1107  Edges[0] = 0;
1108  Edges[1] = 0;
1109  Edges[3] = 0;
1110  }
1111  if (Element[1].Is(RIGID))
1112  {
1113  Edges[0] = 0;
1114  Edges[2] = 0;
1115  Edges[4] = 0;
1116  }
1117  if (Element[2].Is(RIGID))
1118  {
1119  Edges[1] = 0;
1120  Edges[2] = 0;
1121  Edges[5] = 0;
1122  }
1123  if (Element[3].Is(RIGID))
1124  {
1125  Edges[3] = 0;
1126  Edges[4] = 0;
1127  Edges[5] = 0;
1128  }
1129  }
1130 
1131  if ((Edges[0] == 0 && Edges[1] == 0 && Edges[2] == 0 && Edges[3] == 0 && Edges[4] == 0 && Edges[5] == 0) || rigidNodes > 2)
1132  {
1133  dangerousElement = true;
1134  }
1135 
1136  // just to fill the vector
1137  if (dangerousElement == false && toEraseNodeFound == false)
1138  {
1139  SizeType maxCount = 6;
1140  double LargestEdge = 0;
1141 
1142  for (SizeType i = 0; i < 6; i++)
1143  {
1144  if (Edges[i] > LargestEdge)
1145  {
1146  maxCount = i;
1147  LargestEdge = Edges[i];
1148  }
1149  }
1150  if ((CountNodes < ElementsToRefine && LargestEdge > limitEdgeLength && eulerianInletNodes == 0))
1151  {
1152 
1153  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
1154 
1155  bool suitableElement = true;
1156  for (int j = 0; j < CountNodes; j++)
1157  {
1158  const double diffX = std::abs(new_positions[j][0] - new_position[0]) - meanMeshSize * 0.5;
1159  const double diffY = std::abs(new_positions[j][1] - new_position[1]) - meanMeshSize * 0.5;
1160  const double diffZ = std::abs(new_positions[j][2] - new_position[2]) - meanMeshSize * 0.5;
1161  if (diffX < 0 && diffY < 0 && diffZ < 0) // the node is in the same zone of a previously inserted node
1162  {
1163  suitableElement = false;
1164  }
1165  }
1166 
1167  if (suitableElement)
1168  {
1169  nodes_id_to_interpolate[CountNodes][0] = Element[FirstEdgeNode[maxCount]].GetId();
1170  nodes_id_to_interpolate[CountNodes][1] = Element[SecondEdgeNode[maxCount]].GetId();
1171  biggest_volumes[CountNodes] = ElementalVolume;
1172  new_positions[CountNodes] = new_position;
1173  CountNodes++;
1174  }
1175  }
1176  else if (freesurfaceNodes < 4 && rigidNodes < 4 && eulerianInletNodes == 0)
1177  {
1178 
1179  ElementalVolume *= penalization;
1180  for (int nn = 0; nn < ElementsToRefine; nn++)
1181  {
1182  if (ElementalVolume > biggest_volumes[nn])
1183  {
1184  if (maxCount < 6 && LargestEdge > limitEdgeLength)
1185  {
1186  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
1187 
1188  if (ElementsToRefine > 1 && CountNodes > 0)
1189  {
1190  for (int j = 0; j < ElementsToRefine; j++)
1191  {
1192  const double diffX = std::abs(new_positions[j][0] - new_position[0]) - meanMeshSize * 0.5;
1193  const double diffY = std::abs(new_positions[j][1] - new_position[1]) - meanMeshSize * 0.5;
1194  const double diffZ = std::abs(new_positions[j][2] - new_position[2]) - meanMeshSize * 0.5;
1195  if (diffX < 0 && diffY < 0 && diffZ < 0) // the node is in the same zone of a previously inserted node
1196  {
1197  suitableElementForSecondAdd = false;
1198  }
1199  }
1200  }
1201 
1202  if (suitableElementForSecondAdd)
1203  {
1204  nodes_id_to_interpolate[nn][0] = Element[FirstEdgeNode[maxCount]].GetId();
1205  nodes_id_to_interpolate[nn][1] = Element[SecondEdgeNode[maxCount]].GetId();
1206  biggest_volumes[nn] = ElementalVolume;
1207  new_positions[nn] = new_position;
1208  }
1209  }
1210 
1211  break;
1212  }
1213  }
1214  }
1215 
1216  if (eulerianInletNodes > 0 && LargestEdge > (1.7 * meanMeshSize))
1217  {
1218 
1219  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
1220 
1221  bool suitableElement = true;
1222  for (int j = 0; j < CountNodes; j++)
1223  {
1224  const double diffX = std::abs(new_positions[j][0] - new_position[0]) - meanMeshSize * 0.5;
1225  const double diffY = std::abs(new_positions[j][1] - new_position[1]) - meanMeshSize * 0.5;
1226  const double diffZ = std::abs(new_positions[j][2] - new_position[2]) - meanMeshSize * 0.5;
1227  if ((diffX < 0 && diffY < 0 && diffZ < 0) || (Element[FirstEdgeNode[maxCount]].IsNot(INLET) && Element[SecondEdgeNode[maxCount]].IsNot(INLET))) // the node is in the same zone of a previously inserted node
1228  {
1229  suitableElement = false;
1230  }
1231  }
1232 
1233  if (suitableElement)
1234  {
1235  if (CountNodes >= ElementsToRefine)
1236  {
1237  new_positions.resize(CountNodes + 1);
1238  biggest_volumes.resize(CountNodes + 1, false);
1239  nodes_id_to_interpolate.resize(CountNodes + 1);
1240  addedNodesAtEulerianInlet++;
1241  }
1242  nodes_id_to_interpolate[CountNodes][0] = Element[FirstEdgeNode[maxCount]].GetId();
1243  nodes_id_to_interpolate[CountNodes][1] = Element[SecondEdgeNode[maxCount]].GetId();
1244  biggest_volumes[CountNodes] = ElementalVolume;
1245  new_positions[CountNodes] = new_position;
1246  CountNodes++;
1247  }
1248  }
1249  }
1250 
1251  KRATOS_CATCH("")
1252  }
1253 
1254  void SelectEdgeToRefine2DWithRefinement(Element::GeometryType &Element,
1255  std::vector<array_1d<double, 3>> &new_positions,
1256  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
1257  int &CountNodes,
1258  const int ElementsToRefine,
1259  int &nodesInTransitionZone)
1260  {
1261  KRATOS_TRY
1262 
1263  const SizeType nds = Element.size();
1264 
1265  SizeType rigidNodes = 0;
1266  SizeType freesurfaceNodes = 0;
1267  bool toEraseNodeFound = false;
1268  double rigidNodeLocalMeshSize = 0;
1269  double rigidNodeMeshCounter = 0;
1270  double meanMeshSize = mrRemesh.Refine->CriticalRadius;
1271  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
1272  const double currentTime = rCurrentProcessInfo[TIME];
1273  bool insideTransitionZone = false;
1274  for (SizeType pn = 0; pn < nds; ++pn)
1275  {
1276  mMesherUtilities.DefineMeshSizeInTransitionZones2D(mrRemesh, currentTime, Element[pn].Coordinates(), meanMeshSize, insideTransitionZone);
1277  }
1278 
1279  for (SizeType pn = 0; pn < nds; ++pn)
1280  {
1281  if (Element[pn].Is(RIGID))
1282  {
1283  rigidNodes++;
1284  rigidNodeLocalMeshSize += Element[pn].FastGetSolutionStepValue(NODAL_H_WALL);
1285  rigidNodeMeshCounter += 1.0;
1286  }
1287  if (Element[pn].Is(TO_ERASE))
1288  {
1289  toEraseNodeFound = true;
1290  }
1291  if (Element[pn].Is(FREE_SURFACE))
1292  {
1293  freesurfaceNodes++;
1294  }
1295  }
1296 
1297  if (rigidNodeMeshCounter > 0)
1298  {
1299  const double rigidWallMeshSize = rigidNodeLocalMeshSize / rigidNodeMeshCounter;
1300  const double ratio = rigidWallMeshSize / meanMeshSize;
1301  const double tolerance = 1.8;
1302  if (ratio > tolerance)
1303  {
1304  meanMeshSize *= 0.5;
1305  meanMeshSize += 0.5 * rigidWallMeshSize;
1306  }
1307  }
1308  double penalization = 1.0; // penalization here should be greater than 1
1309 
1310  if (freesurfaceNodes > 0)
1311  {
1312  penalization = 1.2; // to avoid to gain too much volume during remeshing step
1313  }
1314  const double safetyCoefficient2D = 1.5;
1315  array_1d<double, 3> Edges(3, 0.0);
1316  array_1d<SizeType, 3> FirstEdgeNode(3, 0);
1317  array_1d<SizeType, 3> SecondEdgeNode(3, 0);
1318  double WallCharacteristicDistance = 0;
1319  array_1d<double, 3> CoorDifference = Element[1].Coordinates() - Element[0].Coordinates();
1320  double SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
1321  Edges[0] = std::sqrt(SquaredLength);
1322  FirstEdgeNode[0] = 0;
1323  SecondEdgeNode[0] = 1;
1324  if ((Element[0].Is(RIGID) && Element[1].Is(RIGID)) || (Element[0].Is(INLET) && Element[1].Is(INLET)))
1325  {
1326  WallCharacteristicDistance = Edges[0];
1327  }
1328  SizeType Counter = 0;
1329  for (SizeType i = 2; i < nds; i++)
1330  {
1331  for (SizeType j = 0; j < i; j++)
1332  {
1333  noalias(CoorDifference) = Element[i].Coordinates() - Element[j].Coordinates();
1334  SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
1335  Counter += 1;
1336  Edges[Counter] = std::sqrt(SquaredLength);
1337  FirstEdgeNode[Counter] = j;
1338  SecondEdgeNode[Counter] = i;
1339  if (Element[i].Is(RIGID) && Element[j].Is(RIGID) && Edges[Counter] > WallCharacteristicDistance)
1340  {
1341  WallCharacteristicDistance = Edges[Counter];
1342  }
1343  }
1344  }
1345 
1346  bool dangerousElement = false;
1347 
1348  for (SizeType i = 0; i < 3; i++)
1349  {
1350  if (rigidNodes > 1)
1351  {
1352  if ((Edges[i] < WallCharacteristicDistance * safetyCoefficient2D && (Element[FirstEdgeNode[i]].Is(RIGID) || Element[SecondEdgeNode[i]].Is(RIGID))) ||
1353  (Element[FirstEdgeNode[i]].Is(RIGID) && Element[SecondEdgeNode[i]].Is(RIGID)))
1354  {
1355  Edges[i] = 0;
1356  }
1357  if ((Element[FirstEdgeNode[i]].Is(FREE_SURFACE) || Element[FirstEdgeNode[i]].Is(RIGID)) &&
1358  (Element[SecondEdgeNode[i]].Is(FREE_SURFACE) || Element[SecondEdgeNode[i]].Is(RIGID)))
1359  {
1360  Edges[i] = 0;
1361  }
1362  }
1363  else if (rigidNodes == 0)
1364  {
1365  SizeType propertyIdFirstNode = Element[FirstEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1366  SizeType propertyIdSecondNode = Element[SecondEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1367  if (propertyIdFirstNode != propertyIdSecondNode)
1368  {
1369  penalization = 1.1; // 10% more than normal nodes
1370  }
1371  }
1372  }
1373  if ((Edges[0] == 0 && Edges[1] == 0 && Edges[2] == 0) || rigidNodes == 3)
1374  {
1375  dangerousElement = true;
1376  }
1377  const double limitEdgeLength = 1.9 * meanMeshSize * penalization;
1378  const double extraLimitEdgeLength = 2.5 * meanMeshSize * penalization;
1379 
1380  if (dangerousElement == false && toEraseNodeFound == false)
1381  {
1382  SizeType maxCount = 3;
1383  double LargestEdge = 0;
1384 
1385  for (SizeType i = 0; i < 3; i++)
1386  {
1387  if (Edges[i] > LargestEdge)
1388  {
1389  maxCount = i;
1390  LargestEdge = Edges[i];
1391  }
1392  }
1393 
1394  if (((CountNodes < (ElementsToRefine + nodesInTransitionZone) || insideTransitionZone == true) && LargestEdge > limitEdgeLength) || LargestEdge > extraLimitEdgeLength)
1395  {
1396  bool newNode = true;
1397  for (SizeType i = 0; i < unsigned(CountNodes); i++)
1398  {
1399  if ((nodes_id_to_interpolate[i][0] == Element[FirstEdgeNode[maxCount]].GetId() && nodes_id_to_interpolate[i][1] == Element[SecondEdgeNode[maxCount]].GetId()) ||
1400  (nodes_id_to_interpolate[i][1] == Element[FirstEdgeNode[maxCount]].GetId() && nodes_id_to_interpolate[i][0] == Element[SecondEdgeNode[maxCount]].GetId()))
1401  {
1402  newNode = false;
1403  }
1404  }
1405  if (newNode == true)
1406  {
1407  new_positions.resize(CountNodes + 1);
1408  nodes_id_to_interpolate.resize(CountNodes + 1);
1409  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
1410  nodes_id_to_interpolate[CountNodes][0] = Element[FirstEdgeNode[maxCount]].GetId();
1411  nodes_id_to_interpolate[CountNodes][1] = Element[SecondEdgeNode[maxCount]].GetId();
1412  new_positions[CountNodes] = new_position;
1413  CountNodes++;
1414  if (insideTransitionZone)
1415  {
1416  nodesInTransitionZone++;
1417  }
1418  }
1419  }
1420  }
1421 
1422  KRATOS_CATCH("")
1423  }
1424 
1425  void SelectEdgeToRefine3DWithRefinement(Element::GeometryType &Element,
1426  std::vector<array_1d<double, 3>> &new_positions,
1427  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
1428  int &CountNodes,
1429  const int ElementsToRefine,
1430  int &nodesInTransitionZone)
1431  {
1432  KRATOS_TRY
1433 
1434  const SizeType nds = Element.size();
1435 
1436  SizeType rigidNodes = 0;
1437  SizeType freesurfaceNodes = 0;
1438  bool toEraseNodeFound = false;
1439 
1440  double meanMeshSize = mrRemesh.Refine->CriticalRadius;
1441  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
1442  double currentTime = rCurrentProcessInfo[TIME];
1443  bool insideTransitionZone = false;
1444  double rigidNodeLocalMeshSize = 0;
1445  double rigidNodeMeshCounter = 0;
1446  for (SizeType pn = 0; pn < nds; ++pn)
1447  {
1448  mMesherUtilities.DefineMeshSizeInTransitionZones3D(mrRemesh, currentTime, Element[pn].Coordinates(), meanMeshSize, insideTransitionZone);
1449  }
1450 
1451  for (SizeType pn = 0; pn < nds; ++pn)
1452  {
1453  if (Element[pn].Is(RIGID))
1454  {
1455  rigidNodes++;
1456  rigidNodeLocalMeshSize += Element[pn].FastGetSolutionStepValue(NODAL_H_WALL);
1457  rigidNodeMeshCounter += 1.0;
1458  }
1459  if (Element[pn].Is(TO_ERASE))
1460  {
1461  toEraseNodeFound = true;
1462  }
1463  if (Element[pn].Is(FREE_SURFACE))
1464  {
1465  freesurfaceNodes++;
1466  }
1467  }
1468 
1469  if (rigidNodeMeshCounter > 0)
1470  {
1471  const double rigidWallMeshSize = rigidNodeLocalMeshSize / rigidNodeMeshCounter;
1472  const double ratio = rigidWallMeshSize / meanMeshSize;
1473  const double tolerance = 1.8;
1474  if (ratio > tolerance)
1475  {
1476  meanMeshSize *= 0.5;
1477  meanMeshSize += 0.5 * rigidWallMeshSize;
1478  }
1479  }
1480 
1481  double penalization = 1.0; // penalization here should be greater than 1
1482 
1483  if (freesurfaceNodes > 0)
1484  {
1485  penalization = 1.2; // to avoid to gain too much volume during remeshing step
1486  }
1487 
1488  const double safetyCoefficient3D = 1.6;
1489  array_1d<double, 6> Edges(6, 0.0);
1490  array_1d<SizeType, 6> FirstEdgeNode(6, 0);
1491  array_1d<SizeType, 6> SecondEdgeNode(6, 0);
1492  double WallCharacteristicDistance = 0;
1493  array_1d<double, 3> CoorDifference = Element[1].Coordinates() - Element[0].Coordinates();
1494  double SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1] + CoorDifference[2] * CoorDifference[2];
1495  Edges[0] = std::sqrt(SquaredLength);
1496  FirstEdgeNode[0] = 0;
1497  SecondEdgeNode[0] = 1;
1498  if ((Element[0].Is(RIGID) && Element[1].Is(RIGID)) || (Element[0].Is(INLET) && Element[1].Is(INLET)))
1499  {
1500  WallCharacteristicDistance = Edges[0];
1501  }
1502  SizeType Counter = 0;
1503  for (SizeType i = 2; i < nds; i++)
1504  {
1505  for (SizeType j = 0; j < i; j++)
1506  {
1507  noalias(CoorDifference) = Element[i].Coordinates() - Element[j].Coordinates();
1508  SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1] + CoorDifference[2] * CoorDifference[2];
1509  Counter += 1;
1510  Edges[Counter] = std::sqrt(SquaredLength);
1511  FirstEdgeNode[Counter] = j;
1512  SecondEdgeNode[Counter] = i;
1513  if (Element[i].Is(RIGID) && Element[j].Is(RIGID) && Edges[Counter] > WallCharacteristicDistance)
1514  {
1515  WallCharacteristicDistance = Edges[Counter];
1516  }
1517  }
1518  }
1519  // Edges connectivity: Edges[0]=d01, Edges[1]=d20, Edges[2]=d21, Edges[3]=d30, Edges[4]=d31, Edges[5]=d32
1520  bool dangerousElement = false;
1521  for (SizeType i = 0; i < 6; i++)
1522  {
1523  if (rigidNodes > 1)
1524  {
1525  if ((Edges[i] < WallCharacteristicDistance * safetyCoefficient3D && (Element[FirstEdgeNode[i]].Is(RIGID) || Element[SecondEdgeNode[i]].Is(RIGID))) ||
1526  (Element[FirstEdgeNode[i]].Is(RIGID) && Element[SecondEdgeNode[i]].Is(RIGID)))
1527  {
1528  Edges[i] = 0;
1529  }
1530  if ((Element[FirstEdgeNode[i]].Is(FREE_SURFACE) || Element[FirstEdgeNode[i]].Is(RIGID)) &&
1531  (Element[SecondEdgeNode[i]].Is(FREE_SURFACE) || Element[SecondEdgeNode[i]].Is(RIGID)))
1532  {
1533  Edges[i] = 0;
1534  }
1535  }
1536  else if (rigidNodes == 0)
1537  {
1538  SizeType propertyIdFirstNode = Element[FirstEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1539  SizeType propertyIdSecondNode = Element[SecondEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1540  if (propertyIdFirstNode != propertyIdSecondNode)
1541  {
1542  penalization = 1.2; // 20% less than normal nodes
1543  }
1544  }
1545  }
1546  if (rigidNodes == 1)
1547  {
1548  if (Element[0].Is(RIGID))
1549  {
1550  Edges[0] = 0;
1551  Edges[1] = 0;
1552  Edges[3] = 0;
1553  }
1554  if (Element[1].Is(RIGID))
1555  {
1556  Edges[0] = 0;
1557  Edges[2] = 0;
1558  Edges[4] = 0;
1559  }
1560  if (Element[2].Is(RIGID))
1561  {
1562  Edges[1] = 0;
1563  Edges[2] = 0;
1564  Edges[5] = 0;
1565  }
1566  if (Element[3].Is(RIGID))
1567  {
1568  Edges[3] = 0;
1569  Edges[4] = 0;
1570  Edges[5] = 0;
1571  }
1572  }
1573 
1574  if ((Edges[0] == 0 && Edges[1] == 0 && Edges[2] == 0 && Edges[3] == 0 && Edges[4] == 0 && Edges[5] == 0) || rigidNodes > 2)
1575  {
1576  dangerousElement = true;
1577  }
1578  const double limitEdgeLength = 1.9 * meanMeshSize * penalization;
1579  const double extraLimitEdgeLength = 2.5 * meanMeshSize * penalization;
1580 
1581  // just to fill the vector
1582  if (dangerousElement == false && toEraseNodeFound == false)
1583  {
1584  SizeType maxCount = 6;
1585  double LargestEdge = 0;
1586  for (SizeType i = 0; i < 6; i++)
1587  {
1588  if (Edges[i] > LargestEdge)
1589  {
1590  maxCount = i;
1591  LargestEdge = Edges[i];
1592  }
1593  }
1594 
1595  if (((CountNodes < (ElementsToRefine + nodesInTransitionZone) || insideTransitionZone == true) && LargestEdge > limitEdgeLength) || LargestEdge > extraLimitEdgeLength)
1596  {
1597  bool newNode = true;
1598  for (SizeType i = 0; i < unsigned(CountNodes); i++)
1599  {
1600  if ((nodes_id_to_interpolate[i][0] == Element[FirstEdgeNode[maxCount]].GetId() && nodes_id_to_interpolate[i][1] == Element[SecondEdgeNode[maxCount]].GetId()) ||
1601  (nodes_id_to_interpolate[i][1] == Element[FirstEdgeNode[maxCount]].GetId() && nodes_id_to_interpolate[i][0] == Element[SecondEdgeNode[maxCount]].GetId()))
1602  {
1603  newNode = false;
1604  }
1605  }
1606  if (newNode == true)
1607  {
1608  new_positions.resize(CountNodes + 1);
1609  nodes_id_to_interpolate.resize(CountNodes + 1);
1610  array_1d<double, 3> new_position = (Element[FirstEdgeNode[maxCount]].Coordinates() + Element[SecondEdgeNode[maxCount]].Coordinates()) * 0.5;
1611  nodes_id_to_interpolate[CountNodes][0] = Element[FirstEdgeNode[maxCount]].GetId();
1612  nodes_id_to_interpolate[CountNodes][1] = Element[SecondEdgeNode[maxCount]].GetId();
1613  new_positions[CountNodes] = new_position;
1614  CountNodes++;
1615  if (insideTransitionZone)
1616  {
1617  nodesInTransitionZone++;
1618  }
1619  }
1620  }
1621  }
1622 
1623  KRATOS_CATCH("")
1624  }
1625 
1626  void CreateAndAddNewNodesInCornerWall(std::vector<array_1d<double, 3>> &new_positions,
1627  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
1628  std::vector<Node::DofsContainerType> &NewDofs,
1629  int ElementsToRefine,
1630  SizeType &maxId)
1631  {
1632  KRATOS_TRY
1633 
1634  const SizeType dimension = mrModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
1635 
1636  std::vector<Node::Pointer> list_of_new_nodes;
1637 
1638  // assign data to dofs
1639  VariablesList &VariablesList = mrModelPart.GetNodalSolutionStepVariablesList();
1640 
1641  for (SizeType nn = 0; nn < new_positions.size(); nn++)
1642  {
1643 
1644  SizeType id = maxId + 1 + nn;
1645 
1646  double x = new_positions[nn][0];
1647  double y = new_positions[nn][1];
1648  double z = 0;
1649  if (dimension == 3)
1650  z = new_positions[nn][2];
1651 
1652  Node::Pointer pnode = mrModelPart.CreateNewNode(id, x, y, z);
1653  pnode->Set(NEW_ENTITY); // not boundary
1654  list_of_new_nodes.push_back(pnode);
1655  if (mrRemesh.InputInitializedFlag)
1656  {
1657  mrRemesh.NodalPreIds.push_back(pnode->Id());
1658  pnode->SetId(id);
1659  }
1660 
1661  // //giving model part variables list to the node
1662  pnode->SetSolutionStepVariablesList(&VariablesList);
1663 
1664  // //set buffer size
1665  pnode->SetBufferSize(mrModelPart.GetBufferSize());
1666 
1667  Node::DofsContainerType &reference_dofs = NewDofs[nn];
1668 
1669  for (Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
1670  {
1671  Node::DofType &rDof = **iii;
1672  pnode->pAddDof(rDof);
1673  }
1674 
1675  auto slave_node_1 = mrModelPart.pGetNode(nodes_id_to_interpolate[nn][0]);
1676  auto slave_node_2 = mrModelPart.pGetNode(nodes_id_to_interpolate[nn][1]);
1677  auto slave_node_3 = mrModelPart.pGetNode(nodes_id_to_interpolate[nn][2]);
1678 
1679  InterpolateFromTwoNodes(pnode, slave_node_1, slave_node_2, VariablesList);
1680 
1681  TakeMaterialPropertiesFromNotRigidNode(pnode, slave_node_3);
1682  }
1683 
1684  // set the coordinates to the original value
1685  const array_1d<double, 3> ZeroNormal(3, 0.0);
1686  for (std::vector<Node::Pointer>::iterator it = list_of_new_nodes.begin(); it != list_of_new_nodes.end(); it++)
1687  {
1688  const array_1d<double, 3> &displacement = (*it)->FastGetSolutionStepValue(DISPLACEMENT);
1689  (*it)->X0() = (*it)->X() - displacement[0];
1690  (*it)->Y0() = (*it)->Y() - displacement[1];
1691  (*it)->Z0() = (*it)->Z() - displacement[2];
1692 
1693  (*it)->Set(FLUID);
1694  (*it)->Set(ACTIVE);
1695  (*it)->Reset(TO_ERASE);
1696  }
1697 
1698  KRATOS_CATCH("")
1699  }
1700 
1701  void CreateAndAddNewNodes(std::vector<array_1d<double, 3>> &new_positions,
1702  std::vector<array_1d<SizeType, 4>> &nodes_id_to_interpolate,
1703  int ElementsToRefine,
1704  SizeType &maxId)
1705  {
1706  KRATOS_TRY
1707 
1708  const SizeType dimension = mrModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
1709 
1710  std::vector<Node::Pointer> list_of_new_nodes;
1711  double NodeIdParent = MesherUtilities::GetMaxNodeId(mrModelPart.GetParentModelPart());
1712  double NodeId = MesherUtilities::GetMaxNodeId(mrModelPart);
1713 
1714  SizeType initial_node_size = NodeIdParent + 1 + ElementsToRefine; // total model part node size
1715 
1716  NodeType::Pointer pnode;
1717  NodeType::DofsContainerType &ReferenceDofs = mrModelPart.Nodes().front().GetDofs();
1718 
1719  if (NodeId > NodeIdParent)
1720  {
1721  initial_node_size = NodeId + 1 + ElementsToRefine;
1722  std::cout << "initial_node_size " << initial_node_size << std::endl;
1723  }
1724 
1725  // assign data to dofs
1726  VariablesList &VariablesList = mrModelPart.GetNodalSolutionStepVariablesList();
1727 
1728  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
1729  SizeType principalModelPartId = rCurrentProcessInfo[MAIN_MATERIAL_PROPERTY];
1730 
1731  for (SizeType nn = 0; nn < new_positions.size(); nn++)
1732  {
1733 
1734  SizeType id = initial_node_size + nn;
1735  maxId = id;
1736  double x = new_positions[nn][0];
1737  double y = new_positions[nn][1];
1738  double z = 0;
1739  if (dimension == 3)
1740  z = new_positions[nn][2];
1741 
1742  // Node::Pointer pnode = mrModelPart.CreateNewNode(id, x, y, z);
1743  pnode = Kratos::make_intrusive<Node>(id, x, y, z);
1744 
1745  pnode->Set(NEW_ENTITY); // not boundary
1746 
1747  if (mrRemesh.InputInitializedFlag)
1748  {
1749  mrRemesh.NodalPreIds.push_back(pnode->Id());
1750  pnode->SetId(id);
1751  }
1752 
1753  // //giving model part variables list to the node
1754  pnode->SetSolutionStepVariablesList(&VariablesList);
1755 
1756  // //set buffer size
1757  pnode->SetBufferSize(mrModelPart.GetBufferSize());
1758 
1759  // Node::DofsContainerType &reference_dofs = NewDofs[nn];
1760  // for (Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
1761  // {
1762  // Node::DofType &rDof = **iii;
1763  // pnode->pAddDof(rDof);
1764  // }
1765  // generating the dofs
1766  for (Node::DofsContainerType::iterator i_dof = ReferenceDofs.begin(); i_dof != ReferenceDofs.end(); ++i_dof)
1767  {
1768  NodeType::DofType &rDof = **i_dof;
1769  NodeType::DofType::Pointer pNewDof = pnode->pAddDof(rDof);
1770 
1771  (pNewDof)->FreeDof();
1772  }
1773 
1774  list_of_new_nodes.push_back(pnode);
1775 
1776  auto slave_node_1 = mrModelPart.pGetNode(nodes_id_to_interpolate[nn][0]);
1777  auto slave_node_2 = mrModelPart.pGetNode(nodes_id_to_interpolate[nn][1]);
1778  InterpolateFromTwoNodes(pnode, slave_node_1, slave_node_2, VariablesList);
1779  if (slave_node_1->Is(RIGID) || slave_node_1->Is(SOLID))
1780  {
1781  TakeMaterialPropertiesFromNotRigidNode(pnode, slave_node_2);
1782  }
1783  else
1784  {
1785 
1786  SizeType propertyIdNodeSlave1 = slave_node_1->FastGetSolutionStepValue(PROPERTY_ID);
1787  if ((mrRemesh.Info->BalancePrincipalSecondaryPartsNodes > 0 && propertyIdNodeSlave1 == principalModelPartId) ||
1788  (mrRemesh.Info->BalancePrincipalSecondaryPartsNodes < 0 && propertyIdNodeSlave1 != principalModelPartId) ||
1789  (slave_node_2->Is(RIGID) || slave_node_2->Is(SOLID)))
1790  {
1791  TakeMaterialPropertiesFromNotRigidNode(pnode, slave_node_1);
1792  }
1793  else
1794  {
1795  TakeMaterialPropertiesFromNotRigidNode(pnode, slave_node_2);
1796  }
1797  }
1798 
1799  SizeType propertyIdNode = pnode->FastGetSolutionStepValue(PROPERTY_ID);
1800  if (propertyIdNode != principalModelPartId)
1801  {
1802  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += 1;
1803  }
1804  }
1805 
1806  // set the coordinates to the original value
1807  const array_1d<double, 3> ZeroNormal(3, 0.0);
1808  for (std::vector<Node::Pointer>::iterator it = list_of_new_nodes.begin(); it != list_of_new_nodes.end(); it++)
1809  {
1810  const array_1d<double, 3> &displacement = (*it)->FastGetSolutionStepValue(DISPLACEMENT);
1811  (*it)->X0() = (*it)->X() - displacement[0];
1812  (*it)->Y0() = (*it)->Y() - displacement[1];
1813  (*it)->Z0() = (*it)->Z() - displacement[2];
1814 
1815  (*it)->Set(FLUID);
1816  (*it)->Set(ACTIVE);
1817  (*it)->Reset(TO_ERASE);
1818 
1819  mrModelPart.Nodes().push_back(*(it));
1820  }
1821 
1822  KRATOS_CATCH("")
1823  }
1824 
1825  void InterpolateFromTwoNodes(Node::Pointer master_node, Node::Pointer slave_node_1, Node::Pointer slave_node_2, VariablesList &rVariablesList)
1826  {
1827 
1828  KRATOS_TRY
1829 
1830  SizeType buffer_size = master_node->GetBufferSize();
1831 
1832  for (VariablesList::const_iterator i_variable = rVariablesList.begin(); i_variable != rVariablesList.end(); i_variable++)
1833  {
1834  std::string variable_name = i_variable->Name();
1835  if (KratosComponents<Variable<double>>::Has(variable_name))
1836  {
1837  const Variable<double> &variable = KratosComponents<Variable<double>>::Get(variable_name);
1838  for (SizeType step = 0; step < buffer_size; step++)
1839  {
1840  // getting the data of the solution step
1841  double &node_data = master_node->FastGetSolutionStepValue(variable, step);
1842 
1843  double node0_data = slave_node_1->FastGetSolutionStepValue(variable, step);
1844  double node1_data = slave_node_2->FastGetSolutionStepValue(variable, step);
1845 
1846  node_data = (0.5 * node0_data + 0.5 * node1_data);
1847  }
1848  }
1849  else if (KratosComponents<Variable<array_1d<double, 3>>>::Has(variable_name))
1850  {
1851  const Variable<array_1d<double, 3>> &variable = KratosComponents<Variable<array_1d<double, 3>>>::Get(variable_name);
1852  for (SizeType step = 0; step < buffer_size; step++)
1853  {
1854  // getting the data of the solution step
1855  array_1d<double, 3> &node_data = master_node->FastGetSolutionStepValue(variable, step);
1856 
1857  const array_1d<double, 3> &node0_data = slave_node_1->FastGetSolutionStepValue(variable, step);
1858  const array_1d<double, 3> &node1_data = slave_node_2->FastGetSolutionStepValue(variable, step);
1859 
1860  noalias(node_data) = (0.5 * node0_data + 0.5 * node1_data);
1861  // node_data = (0.5*node0_data + 0.5*node1_data);
1862  }
1863  }
1864  else if (KratosComponents<Variable<int>>::Has(variable_name))
1865  {
1866  // std::cout<<"int"<<std::endl;
1867  // NO INTERPOLATION
1868  }
1869  else if (KratosComponents<Variable<bool>>::Has(variable_name))
1870  {
1871  // std::cout<<"bool"<<std::endl;
1872  // NO INTERPOLATION
1873  }
1874  else if (KratosComponents<Variable<Matrix>>::Has(variable_name))
1875  {
1876  // std::cout<<"Matrix"<<std::endl;
1877  const Variable<Matrix> &variable = KratosComponents<Variable<Matrix>>::Get(variable_name);
1878  for (SizeType step = 0; step < buffer_size; step++)
1879  {
1880  // getting the data of the solution step
1881  Matrix &node_data = master_node->FastGetSolutionStepValue(variable, step);
1882 
1883  Matrix &node0_data = slave_node_1->FastGetSolutionStepValue(variable, step);
1884  Matrix &node1_data = slave_node_2->FastGetSolutionStepValue(variable, step);
1885 
1886  if (node_data.size1() > 0 && node_data.size2())
1887  {
1888  if (node_data.size1() == node0_data.size1() && node_data.size2() == node0_data.size2() &&
1889  node_data.size1() == node1_data.size1() && node_data.size2() == node1_data.size2())
1890  {
1891 
1892  noalias(node_data) = (0.5 * node0_data + 0.5 * node1_data);
1893  // node_data = (0.5*node0_data + 0.5*node1_data);
1894  }
1895  }
1896  }
1897  }
1898  else if (KratosComponents<Variable<Vector>>::Has(variable_name))
1899  {
1900  // std::cout<<"Vector"<<std::endl;
1901  const Variable<Vector> &variable = KratosComponents<Variable<Vector>>::Get(variable_name);
1902  for (SizeType step = 0; step < buffer_size; step++)
1903  {
1904  // getting the data of the solution step
1905  Vector &node_data = master_node->FastGetSolutionStepValue(variable, step);
1906 
1907  Vector &node0_data = slave_node_1->FastGetSolutionStepValue(variable, step);
1908  Vector &node1_data = slave_node_2->FastGetSolutionStepValue(variable, step);
1909 
1910  if (node_data.size() > 0)
1911  {
1912  if (node_data.size() == node0_data.size() &&
1913  node_data.size() == node1_data.size())
1914  {
1915 
1916  noalias(node_data) = (0.5 * node0_data + 0.5 * node1_data);
1917  // node_data = (0.5*node0_data + 0.5*node1_data);
1918  }
1919  }
1920  }
1921  }
1922  }
1923 
1924  KRATOS_CATCH("")
1925  }
1926 
1927  void TakeMaterialPropertiesFromNotRigidNode(Node::Pointer master_node, Node::Pointer SlaveNode)
1928  {
1929 
1930  KRATOS_TRY
1931  master_node->FastGetSolutionStepValue(PROPERTY_ID) = SlaveNode->FastGetSolutionStepValue(PROPERTY_ID);
1932  if (master_node->SolutionStepsDataHas(BULK_MODULUS) && SlaveNode->SolutionStepsDataHas(BULK_MODULUS))
1933  master_node->FastGetSolutionStepValue(BULK_MODULUS) = SlaveNode->FastGetSolutionStepValue(BULK_MODULUS);
1934  if (master_node->SolutionStepsDataHas(DENSITY) && SlaveNode->SolutionStepsDataHas(DENSITY))
1935  master_node->FastGetSolutionStepValue(DENSITY) = SlaveNode->FastGetSolutionStepValue(DENSITY);
1936  if (master_node->SolutionStepsDataHas(DYNAMIC_VISCOSITY) && SlaveNode->SolutionStepsDataHas(DYNAMIC_VISCOSITY))
1937  master_node->FastGetSolutionStepValue(DYNAMIC_VISCOSITY) = SlaveNode->FastGetSolutionStepValue(DYNAMIC_VISCOSITY);
1938 
1939  if (master_node->SolutionStepsDataHas(YIELD_SHEAR) && SlaveNode->SolutionStepsDataHas(YIELD_SHEAR))
1940  master_node->FastGetSolutionStepValue(YIELD_SHEAR) = SlaveNode->FastGetSolutionStepValue(YIELD_SHEAR);
1941 
1942  if (master_node->SolutionStepsDataHas(FLOW_INDEX) && SlaveNode->SolutionStepsDataHas(FLOW_INDEX))
1943  master_node->FastGetSolutionStepValue(FLOW_INDEX) = SlaveNode->FastGetSolutionStepValue(FLOW_INDEX);
1944  if (master_node->SolutionStepsDataHas(ADAPTIVE_EXPONENT) && SlaveNode->SolutionStepsDataHas(ADAPTIVE_EXPONENT))
1945  master_node->FastGetSolutionStepValue(ADAPTIVE_EXPONENT) = SlaveNode->FastGetSolutionStepValue(ADAPTIVE_EXPONENT);
1946  if (master_node->SolutionStepsDataHas(STATIC_FRICTION) && SlaveNode->SolutionStepsDataHas(STATIC_FRICTION))
1947  master_node->FastGetSolutionStepValue(STATIC_FRICTION) = SlaveNode->FastGetSolutionStepValue(STATIC_FRICTION);
1948  if (master_node->SolutionStepsDataHas(DYNAMIC_FRICTION) && SlaveNode->SolutionStepsDataHas(DYNAMIC_FRICTION))
1949  master_node->FastGetSolutionStepValue(DYNAMIC_FRICTION) = SlaveNode->FastGetSolutionStepValue(DYNAMIC_FRICTION);
1950  if (master_node->SolutionStepsDataHas(INERTIAL_NUMBER_ZERO) && SlaveNode->SolutionStepsDataHas(INERTIAL_NUMBER_ZERO))
1951  master_node->FastGetSolutionStepValue(INERTIAL_NUMBER_ZERO) = SlaveNode->FastGetSolutionStepValue(INERTIAL_NUMBER_ZERO);
1952  if (master_node->SolutionStepsDataHas(GRAIN_DIAMETER) && SlaveNode->SolutionStepsDataHas(GRAIN_DIAMETER))
1953  master_node->FastGetSolutionStepValue(GRAIN_DIAMETER) = SlaveNode->FastGetSolutionStepValue(GRAIN_DIAMETER);
1954  if (master_node->SolutionStepsDataHas(GRAIN_DENSITY) && SlaveNode->SolutionStepsDataHas(GRAIN_DENSITY))
1955  master_node->FastGetSolutionStepValue(GRAIN_DENSITY) = SlaveNode->FastGetSolutionStepValue(GRAIN_DENSITY);
1956  if (master_node->SolutionStepsDataHas(REGULARIZATION_COEFFICIENT) && SlaveNode->SolutionStepsDataHas(REGULARIZATION_COEFFICIENT))
1957  master_node->FastGetSolutionStepValue(REGULARIZATION_COEFFICIENT) = SlaveNode->FastGetSolutionStepValue(REGULARIZATION_COEFFICIENT);
1958 
1959  if (master_node->SolutionStepsDataHas(INTERNAL_FRICTION_ANGLE) && SlaveNode->SolutionStepsDataHas(INTERNAL_FRICTION_ANGLE))
1960  master_node->FastGetSolutionStepValue(INTERNAL_FRICTION_ANGLE) = SlaveNode->FastGetSolutionStepValue(INTERNAL_FRICTION_ANGLE);
1961  if (master_node->SolutionStepsDataHas(COHESION) && SlaveNode->SolutionStepsDataHas(COHESION))
1962  master_node->FastGetSolutionStepValue(COHESION) = SlaveNode->FastGetSolutionStepValue(COHESION);
1963 
1964  if (master_node->SolutionStepsDataHas(DEVIATORIC_COEFFICIENT) && SlaveNode->SolutionStepsDataHas(DEVIATORIC_COEFFICIENT))
1965  {
1966  master_node->FastGetSolutionStepValue(DEVIATORIC_COEFFICIENT) = SlaveNode->FastGetSolutionStepValue(DEVIATORIC_COEFFICIENT);
1967  master_node->FastGetSolutionStepValue(VOLUMETRIC_COEFFICIENT) = SlaveNode->FastGetSolutionStepValue(VOLUMETRIC_COEFFICIENT);
1968  }
1969 
1970  if (master_node->SolutionStepsDataHas(YOUNG_MODULUS) && SlaveNode->SolutionStepsDataHas(YOUNG_MODULUS))
1971  {
1972  master_node->FastGetSolutionStepValue(YOUNG_MODULUS) = 0;
1973  master_node->FastGetSolutionStepValue(POISSON_RATIO) = 0;
1974  }
1975  if (master_node->SolutionStepsDataHas(SOLID_DENSITY) && SlaveNode->SolutionStepsDataHas(SOLID_DENSITY))
1976  {
1977  master_node->FastGetSolutionStepValue(SOLID_DENSITY) = 0;
1978  master_node->Reset(SOLID);
1979  master_node->FastGetSolutionStepValue(INTERFACE_NODE) = false;
1980  }
1981 
1982  KRATOS_CATCH("")
1983  }
1984 
1988 
1992 
1996 
1998  GenerateNewNodesBeforeMeshingProcess &operator=(GenerateNewNodesBeforeMeshingProcess const &rOther);
1999 
2001 
2003  // Process(Process const& rOther);
2004 
2006 
2007  }; // Class Process
2008 
2010 
2013 
2017 
2019  inline std::istream &operator>>(std::istream &rIStream,
2020  GenerateNewNodesBeforeMeshingProcess &rThis);
2021 
2023  inline std::ostream &operator<<(std::ostream &rOStream,
2024  const GenerateNewNodesBeforeMeshingProcess &rThis)
2025  {
2026  rThis.PrintInfo(rOStream);
2027  rOStream << std::endl;
2028  rThis.PrintData(rOStream);
2029 
2030  return rOStream;
2031  }
2033 
2034 } // namespace Kratos.
2035 
2036 #endif // KRATOS_GENERATE_NEW_NODES_BEFORE_MESHING_PROCESS_H_INCLUDED defined
Kratos::Condition
Base class for all Conditions.
Definition: condition.h:59
Kratos::Dof
Dof represents a degree of freedom (DoF).
Definition: dof.h:86
Kratos::Dof::Pointer
Dof * Pointer
Pointer definition of Dof.
Definition: dof.h:93
Kratos::Element::GeometryType
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: element.h:83
Kratos::Flags::Is
bool Is(Flags const &rOther) const
Definition: flags.h:274
Kratos::Flags::IsNot
bool IsNot(Flags const &rOther) const
Definition: flags.h:291
Kratos::GenerateNewNodesBeforeMeshingProcess
Refine Mesh Elements Process 2D and 3D.
Definition: generate_new_nodes_before_meshing_process.hpp:47
Kratos::GenerateNewNodesBeforeMeshingProcess::PropertiesType
ModelPart::PropertiesType PropertiesType
Definition: generate_new_nodes_before_meshing_process.hpp:57
Kratos::GenerateNewNodesBeforeMeshingProcess::Info
std::string Info() const override
Turn back information as a string.
Definition: generate_new_nodes_before_meshing_process.hpp:282
Kratos::GenerateNewNodesBeforeMeshingProcess::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(GenerateNewNodesBeforeMeshingProcess)
Pointer definition of Process.
Kratos::GenerateNewNodesBeforeMeshingProcess::NodeType
ModelPart::NodeType NodeType
Definition: generate_new_nodes_before_meshing_process.hpp:55
Kratos::GenerateNewNodesBeforeMeshingProcess::~GenerateNewNodesBeforeMeshingProcess
virtual ~GenerateNewNodesBeforeMeshingProcess()
Destructor.
Definition: generate_new_nodes_before_meshing_process.hpp:76
Kratos::GenerateNewNodesBeforeMeshingProcess::GeometryType
ConditionType::GeometryType GeometryType
Definition: generate_new_nodes_before_meshing_process.hpp:58
Kratos::GenerateNewNodesBeforeMeshingProcess::ConditionType
ModelPart::ConditionType ConditionType
Definition: generate_new_nodes_before_meshing_process.hpp:56
Kratos::GenerateNewNodesBeforeMeshingProcess::Execute
void Execute() override
Execute method is used to execute the Process algorithms.
Definition: generate_new_nodes_before_meshing_process.hpp:93
Kratos::GenerateNewNodesBeforeMeshingProcess::GenerateNewNodesBeforeMeshingProcess
GenerateNewNodesBeforeMeshingProcess(ModelPart &rModelPart, MesherUtilities::MeshingParameters &rRemeshingParameters, int EchoLevel)
Default constructor.
Definition: generate_new_nodes_before_meshing_process.hpp:66
Kratos::GenerateNewNodesBeforeMeshingProcess::PrintData
void PrintData(std::ostream &rOStream) const override
Print object's data.s.
Definition: generate_new_nodes_before_meshing_process.hpp:294
Kratos::GenerateNewNodesBeforeMeshingProcess::PrintInfo
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: generate_new_nodes_before_meshing_process.hpp:288
Kratos::GenerateNewNodesBeforeMeshingProcess::SizeType
std::size_t SizeType
Definition: generate_new_nodes_before_meshing_process.hpp:59
Kratos::GenerateNewNodesBeforeMeshingProcess::operator()
void operator()()
This operator is provided to call the process as a function and simply calls the Execute method.
Definition: generate_new_nodes_before_meshing_process.hpp:83
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::MesherProcess
The base class for processes passed to the solution scheme.
Definition: mesher_process.hpp:37
Kratos::MesherUtilities
Short class definition.
Definition: mesher_utilities.hpp:49
Kratos::MesherUtilities::GetMaxNodeId
static unsigned int GetMaxNodeId(ModelPart &rModelPart)
Definition: mesher_utilities.hpp:1408
Kratos::MesherUtilities::DefineMeshSizeInTransitionZones2D
void DefineMeshSizeInTransitionZones2D(MeshingParameters &rMeshingVariables, double currentTime, array_1d< double, 3 > NodeCoordinates, double &meanMeshSize, bool &insideTransitionZone)
Definition: mesher_utilities.cpp:2153
Kratos::MesherUtilities::DefineMeshSizeInTransitionZones3D
void DefineMeshSizeInTransitionZones3D(MeshingParameters &rMeshingVariables, double currentTime, array_1d< double, 3 > NodeCoordinates, double &meanMeshSize, bool &insideTransitionZone)
Definition: mesher_utilities.cpp:2311
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::Name
std::string & Name()
Definition: model_part.h:1811
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::ElementsEnd
ElementIterator ElementsEnd(IndexType ThisIndex=0)
Definition: model_part.h:1179
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::Node::DofType
Dof< double > DofType
Dof type.
Definition: node.h:83
Kratos::Node::DofsContainerType
std::vector< std::unique_ptr< Dof< double > >> DofsContainerType
The DoF container type definition.
Definition: node.h:92
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::Properties
Properties encapsulates data shared by different Elements or Conditions. It can store any type of dat...
Definition: properties.h:69
Kratos::VariablesList::const_iterator
boost::indirect_iterator< VariablesContainerType::const_iterator > const_iterator
Definition: variables_list.h:77
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
mesh_error_calculation_utilities.hpp
mesher_process.hpp
mesher_utilities.hpp
model_part.h
EMPIRE_API_helpers::EchoLevel
static int EchoLevel
Definition: co_sim_EMPIRE_API.h:42
GenerateWind.z
z
Definition: GenerateWind.py:163
Kratos::ModelerFactory::Has
bool Has(const std::string &ModelerName)
Checks if the modeler is registered.
Definition: modeler_factory.cpp:24
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::Vector
Internals::Matrix< double, AMatrix::dynamic, 1 > Vector
Definition: amatrix_interface.h:472
Kratos::unique_ptr
std::unique_ptr< T > unique_ptr
Definition: smart_pointers.h:33
Kratos::Matrix
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
Kratos::operator>>
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::int
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
Kratos::operator<<
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
edgebased_PureConvection.ProcessInfo
ProcessInfo
Definition: edgebased_PureConvection.py:116
face_heat.step
int step
Definition: face_heat.py:88
generate_axisymmetric_navier_stokes_element.y
y
Other simbols definition.
Definition: generate_axisymmetric_navier_stokes_element.py:54
generate_droplet_dynamics.pn
pn
Definition: generate_droplet_dynamics.py:65
isotropic_damage_automatic_differentiation.dimension
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
quadrature.j
int j
Definition: quadrature.py:648
sensitivityMatrix.x
x
Definition: sensitivityMatrix.py:49
tensors::i
integer i
Definition: TensorModule.f:17
spatial_containers.h
Kratos::MesherUtilities::MeshingParameters
Definition: mesher_utilities.hpp:631
Kratos::MesherUtilities::MeshingParameters::Info
MeshingInfoParameters::Pointer Info
Definition: mesher_utilities.hpp:681
Kratos::MesherUtilities::MeshingParameters::Refine
RefiningParameters::Pointer Refine
Definition: mesher_utilities.hpp:684
Kratos::MesherUtilities::MeshingParameters::RefiningBoxFinalTime
std::vector< double > RefiningBoxFinalTime
Definition: mesher_utilities.hpp:706
Kratos::MesherUtilities::MeshingParameters::UseRefiningBox
std::vector< bool > UseRefiningBox
Definition: mesher_utilities.hpp:704
Kratos::MesherUtilities::MeshingParameters::InputInitializedFlag
bool InputInitializedFlag
Definition: mesher_utilities.hpp:661
Kratos::MesherUtilities::MeshingParameters::SubModelPartName
std::string SubModelPartName
Definition: mesher_utilities.hpp:642
Kratos::MesherUtilities::MeshingParameters::NodalPreIds
std::vector< int > NodalPreIds
Definition: mesher_utilities.hpp:668
Kratos::MesherUtilities::MeshingParameters::ExecutionOptions
Flags ExecutionOptions
Definition: mesher_utilities.hpp:648
Kratos::MesherUtilities::MeshingParameters::RefiningBoxInitialTime
std::vector< double > RefiningBoxInitialTime
Definition: mesher_utilities.hpp:705
variables_list_data_value_container.h
Generated on Sat Feb 3 2024 02:08:33 for KratosMultiphysics by doxygen 1.9.1