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.
remove_mesh_nodes_for_fluids_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_REMOVE_MESH_NODES_FOR_FLUIDS_PROCESS_H_INCLUDED)
11 #define KRATOS_REMOVE_MESH_NODES_FOR_FLUIDS_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 
25 #include "pfem_fluid_dynamics_application_variables.h"
26 #include "custom_processes/mesher_process.hpp"
27 
29 // Data: NORMAL, MASTER_NODES, NEIGHBOUR_NODES, NEIGBOUR_ELEMENTS
30 // Flags: (checked) TO_ERASE, BOUNDARY, STRUCTURE, TO_SPLIT, CONTACT, NEW_ENTITY, BLOCKED
31 // (set) TO_ERASE(conditions,nodes)(set), NEW_ENTITY(conditions,nodes)(set), BLOCKED(nodes)->locally, VISITED(nodes)(set)
32 // (modified)
33 // (reset) BLOCKED->locally
34 //(set):=(set in this process)
35 
36 namespace Kratos
37 {
38 
41 
43 
53  class RemoveMeshNodesForFluidsProcess
54  : public MesherProcess
55  {
56  public:
59 
61  KRATOS_CLASS_POINTER_DEFINITION(RemoveMeshNodesForFluidsProcess);
62 
63  typedef ModelPart::ConditionType ConditionType;
64  typedef ModelPart::PropertiesType PropertiesType;
65  typedef ConditionType::GeometryType GeometryType;
66  typedef Bucket<3, Node, std::vector<Node::Pointer>, Node::Pointer, std::vector<Node::Pointer>::iterator, std::vector<double>::iterator> BucketType;
67  typedef Tree<KDTreePartition<BucketType>> KdtreeType; // Kdtree
68  typedef ModelPart::MeshType::GeometryType::PointsArrayType PointsArrayType;
69  typedef GlobalPointersVector<Node> NodeWeakPtrVectorType;
70  typedef GlobalPointersVector<Element> ElementWeakPtrVectorType;
71  typedef std::size_t SizeType;
75 
77  RemoveMeshNodesForFluidsProcess(ModelPart &rModelPart,
78  MesherUtilities::MeshingParameters &rRemeshingParameters,
79  int EchoLevel)
80  : mrModelPart(rModelPart),
81  mrRemesh(rRemeshingParameters)
82  {
83  KRATOS_INFO("RemoveMeshNodesForFluidsProcess") << " activated " << std::endl;
84 
85  mEchoLevel = EchoLevel;
86  }
87 
89  virtual ~RemoveMeshNodesForFluidsProcess() {}
90 
94 
96  void operator()()
97  {
98  Execute();
99  }
100 
104 
106  void Execute() override
107  {
108 
109  KRATOS_TRY
110 
111  if (mEchoLevel > 1)
112  {
113  std::cout << " [ REMOVE CLOSE NODES: " << std::endl;
114  }
115 
116  SizeType inside_nodes_removed = 0;
117  SizeType boundary_nodes_removed = 0;
118  SizeType nodes_removed_inlet_zone = 0;
119 
120  // if the remove_node switch is activated, we check if the nodes got too close
121  if (mrRemesh.Refine->RemovingOptions.Is(MesherUtilities::REMOVE_NODES))
122  {
123  bool some_node_is_removed = false;
124 
125  if (mEchoLevel > 1)
126  std::cout << " REMOVE_NODES is TRUE " << std::endl;
127 
128  if (mrRemesh.Refine->RemovingOptions.Is(MesherUtilities::REMOVE_NODES_ON_DISTANCE))
129  {
130  if (mEchoLevel > 1)
131  std::cout << " REMOVE_NODES_ON_DISTANCE is TRUE " << std::endl;
132 
133  some_node_is_removed = RemoveNodesOnDistance(inside_nodes_removed, boundary_nodes_removed, nodes_removed_inlet_zone);
134  }
135 
136  if (some_node_is_removed || mrRemesh.UseBoundingBox == true)
137  this->CleanRemovedNodes(mrModelPart);
138  }
139 
140  // number of removed nodes:
141  mrRemesh.Info->RemovedNodes += inside_nodes_removed + boundary_nodes_removed - nodes_removed_inlet_zone;
142  int distance_remove = inside_nodes_removed + boundary_nodes_removed;
143 
144  if (mEchoLevel > 1)
145  {
146  std::cout << " [ NODES ( removed : " << mrRemesh.Info->RemovedNodes << " ) ]" << std::endl;
147  std::cout << " [ Distance(removed: " << distance_remove << "; inside: " << inside_nodes_removed << "; boundary: " << boundary_nodes_removed << ") ]" << std::endl;
148  std::cout << " REMOVE CLOSE NODES ]; " << std::endl;
149  }
150 
151  KRATOS_CATCH(" ")
152  }
153 
157 
161 
165 
167  std::string Info() const override
168  {
169  return "RemoveMeshNodesForFluidsProcess";
170  }
171 
173  void PrintInfo(std::ostream &rOStream) const override
174  {
175  rOStream << "RemoveMeshNodesForFluidsProcess";
176  }
177 
179  void PrintData(std::ostream &rOStream) const override
180  {
181  }
182 
186 
188 
189  private:
192 
196  ModelPart &mrModelPart;
197 
198  MesherUtilities::MeshingParameters &mrRemesh;
199 
200  MesherUtilities mMesherUtilities;
201 
202  int mEchoLevel;
203 
207 
208  //**************************************************************************
209  //**************************************************************************
210 
211  void CleanRemovedNodes(ModelPart &rModelPart)
212  {
213  KRATOS_TRY
214 
215  // MESH 0 total domain mesh
216  ModelPart::NodesContainerType temporal_nodes;
217  temporal_nodes.reserve(rModelPart.Nodes().size());
218 
219  temporal_nodes.swap(rModelPart.Nodes());
220 
221  for (ModelPart::NodesContainerType::iterator i_node = temporal_nodes.begin(); i_node != temporal_nodes.end(); i_node++)
222  {
223  if (i_node->IsNot(TO_ERASE))
224  {
225 
227  bool boundingBox = mrRemesh.UseBoundingBox;
228  if (boundingBox == true && i_node->IsNot(RIGID))
229  {
230  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
231  double currentTime = rCurrentProcessInfo[TIME];
232  double initialTime = mrRemesh.BoundingBoxInitialTime;
233  double finalTime = mrRemesh.BoundingBoxFinalTime;
234  if (currentTime > initialTime && currentTime < finalTime)
235  {
236  array_1d<double, 3> BoundingBoxLowerPoint = mrRemesh.BoundingBoxLowerPoint;
237  array_1d<double, 3> BoundingBoxUpperPoint = mrRemesh.BoundingBoxUpperPoint;
238  if (i_node->X() < BoundingBoxLowerPoint[0] || i_node->Y() < BoundingBoxLowerPoint[1] || i_node->Z() < BoundingBoxLowerPoint[2] ||
239  i_node->X() > BoundingBoxUpperPoint[0] || i_node->Y() > BoundingBoxUpperPoint[1] || i_node->Z() > BoundingBoxUpperPoint[2])
240  {
241  i_node->Set(TO_ERASE);
242  }
243  else
244  {
245  (rModelPart.Nodes()).push_back(*(i_node.base()));
246  }
247  }
248  else
249  {
250  (rModelPart.Nodes()).push_back(*(i_node.base()));
251  }
252  }
253  else
254  {
255  (rModelPart.Nodes()).push_back(*(i_node.base()));
256  }
257 
259  }
260  }
261 
262  rModelPart.Nodes().Sort();
263 
264  KRATOS_CATCH("")
265  }
266 
267  //**************************************************************************
268  //**************************************************************************
269 
270  bool RemoveNodesOnDistance(SizeType &inside_nodes_removed,
271  SizeType &boundary_nodes_removed,
272  SizeType &nodes_removed_inlet_zone)
273  {
274  KRATOS_TRY
275 
276  const SizeType dimension = mrModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
277 
278  bool some_node_is_removed = false;
279 
280  // bucket size definition:
281  SizeType bucket_size = 20;
282 
283  // create the list of the nodes to be check during the search
284  std::vector<Node::Pointer> list_of_nodes;
285  list_of_nodes.reserve(mrModelPart.NumberOfNodes());
286  for (ModelPart::NodesContainerType::iterator i_node = mrModelPart.NodesBegin(); i_node != mrModelPart.NodesEnd(); i_node++)
287  {
288  (list_of_nodes).push_back(*(i_node.base()));
289  }
290 
291  KdtreeType nodes_tree(list_of_nodes.begin(), list_of_nodes.end(), bucket_size);
293  // all of the nodes in this list will be preserved
294  SizeType num_neighbours = 100;
295  std::vector<Node::Pointer> neighbours(num_neighbours);
296  std::vector<double> neighbour_distances(num_neighbours);
297 
298  // radius means the distance, if the distance between two nodes is closer to radius -> mark for removing
299  double radius = 0;
300  Node work_point(0, 0.0, 0.0, 0.0);
301  SizeType n_points_in_radius;
302 
303  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
304  const double currentTime = rCurrentProcessInfo[TIME];
305  double meshSize = mrRemesh.Refine->CriticalRadius;
306 
307  bool refiningBox = false;
308  for (SizeType index = 0; index < mrRemesh.UseRefiningBox.size(); index++)
309  {
310  if (mrRemesh.UseRefiningBox[index] == true && currentTime > mrRemesh.RefiningBoxInitialTime[index] && currentTime < mrRemesh.RefiningBoxFinalTime[index])
311  {
312  refiningBox = true;
313  }
314  }
315 
316  const SizeType principalModelPartId = rCurrentProcessInfo[MAIN_MATERIAL_PROPERTY];
317  for (ModelPart::ElementsContainerType::const_iterator ie = mrModelPart.ElementsBegin();
318  ie != mrModelPart.ElementsEnd(); ie++)
319  {
320  SizeType rigidNodes = 0;
321  // coordinates
322  for (SizeType i = 0; i < ie->GetGeometry().size(); i++)
323  {
324  if ((ie->GetGeometry()[i].Is(RIGID) && ie->GetGeometry()[i].IsNot(PFEMFlags::LAGRANGIAN_INLET)) || ie->GetGeometry()[i].Is(SOLID) || ie->GetGeometry()[i].Is(PFEMFlags::EULERIAN_INLET))
325  {
326  rigidNodes++;
327  }
328  }
329  if (dimension == 2)
330  {
331  if (rigidNodes > 0)
332  EraseCriticalNodes2D(ie->GetGeometry(), inside_nodes_removed, nodes_removed_inlet_zone);
333  }
334  else if (dimension == 3)
335  {
336  if (rigidNodes > 1)
337  EraseCriticalNodes3D(ie->GetGeometry(), inside_nodes_removed, nodes_removed_inlet_zone, rigidNodes);
338  }
339  }
340 
341  for (ModelPart::NodesContainerType::const_iterator in = mrModelPart.NodesBegin(); in != mrModelPart.NodesEnd(); in++)
342  {
343 
344  const SizeType propertyIdNode = in->FastGetSolutionStepValue(PROPERTY_ID);
345  meshSize = mrRemesh.Refine->CriticalRadius;
346  double rigidNodeLocalMeshSize = 0;
347  double rigidNodeMeshCounter = 0;
348  NodeWeakPtrVectorType &rN = in->GetValue(NEIGHBOUR_NODES);
349 
350  for (SizeType i = 0; i < rN.size(); i++)
351  {
352  if (rN[i].Is(RIGID))
353  {
354  rigidNodeLocalMeshSize += rN[i].FastGetSolutionStepValue(NODAL_H_WALL);
355  rigidNodeMeshCounter += 1.0;
356  }
357  }
358 
359  if (refiningBox == true)
360  {
361  array_1d<double, 3> NodeCoordinates = in->Coordinates();
362  if (dimension == 2)
363  {
364  bool insideTransitionZone = false;
365  mMesherUtilities.DefineMeshSizeInTransitionZones2D(mrRemesh, currentTime, NodeCoordinates, meshSize, insideTransitionZone);
366  }
367  else if (dimension == 3)
368  {
369  bool insideTransitionZone = false;
370  mMesherUtilities.DefineMeshSizeInTransitionZones3D(mrRemesh, currentTime, NodeCoordinates, meshSize, insideTransitionZone);
371  }
372  }
373 
374  if (rigidNodeMeshCounter > 0)
375  {
376  const double rigidWallMeshSize = rigidNodeLocalMeshSize / rigidNodeMeshCounter;
377  const double ratio = rigidWallMeshSize / meshSize;
378  const double tolerance = 2.0;
379  if (ratio > tolerance)
380  {
381  meshSize *= 0.5;
382  meshSize += 0.5 * rigidWallMeshSize;
383  }
384  }
385 
386  const double size_for_distance_boundary = 0.6 * meshSize;
387  const double size_for_wall_tip_contact_side = 0.15 * mrRemesh.Refine->CriticalSide;
388 
389  if (in->Is(TO_ERASE))
390  {
391  some_node_is_removed = true;
392  }
393  bool on_contact_tip = false;
394 
395  if (in->Is(TO_SPLIT) || in->Is(CONTACT))
396  on_contact_tip = true;
397 
398  if (in->IsNot(NEW_ENTITY) && in->IsNot(INLET) && in->IsNot(ISOLATED) && in->IsNot(RIGID) && in->IsNot(SOLID))
399  {
400  SizeType neighErasedNodes = 0;
401  SizeType freeSurfaceNeighNodes = 0;
402  bool inletElement = false;
403  bool interfaceElement = false;
404 
405  radius = 0.6 * meshSize;
406  work_point[0] = in->X();
407  work_point[1] = in->Y();
408  work_point[2] = in->Z();
409 
410  if (in->Is(FREE_SURFACE))
411  {
412  // it must be more difficult to erase a free_surface node, otherwise, lot of volume is lost
413  // this value has a strong effect on volume variation due to remeshing
414  radius = 0.475 * meshSize; // compared with element radius
415  NodeWeakPtrVectorType &neighb_nodes = in->GetValue(NEIGHBOUR_NODES);
416  SizeType countRigid = 0;
417 
418  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
419  {
420  if ((nn)->Is(RIGID) || (nn)->Is(SOLID))
421  {
422  countRigid++;
423  }
424  if ((nn)->Is(TO_ERASE))
425  {
426  neighErasedNodes++;
427  }
428  }
429  if (countRigid == neighb_nodes.size())
430  {
431  radius = 0.15 * meshSize;
432  }
433  }
434  else
435  {
436  NodeWeakPtrVectorType &neighb_nodes = in->GetValue(NEIGHBOUR_NODES);
437  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
438  {
439  const SizeType propertyIdSecondNode = (nn)->FastGetSolutionStepValue(PROPERTY_ID);
440  if ((nn)->Is(FREE_SURFACE))
441  {
442  freeSurfaceNeighNodes++;
443  }
444  if ((nn)->Is(TO_ERASE))
445  {
446  neighErasedNodes++;
447  }
448  if ((nn)->Is(PFEMFlags::EULERIAN_INLET))
449  {
450  inletElement = true;
451  }
452  if (propertyIdNode != propertyIdSecondNode && (nn)->IsNot(RIGID))
453  {
454  interfaceElement = true;
455  }
456  }
457  }
458 
459  if (freeSurfaceNeighNodes > 1)
460  {
461  radius = 0.5 * meshSize;
462  }
463  else if (interfaceElement == true)
464  {
465  if (dimension == 2)
466  radius = 0.54 * meshSize; // 10% less than normal nodes
467  if (dimension == 3)
468  radius = 0.48 * meshSize; // 20% less than normal nodes
469  }
470 
471  n_points_in_radius = nodes_tree.SearchInRadius(work_point, radius, neighbours.begin(), neighbour_distances.begin(), num_neighbours);
472 
473  if (n_points_in_radius > 1 && neighErasedNodes == 0 && in->IsNot(INLET))
474  {
475 
476  if (in->IsNot(RIGID) && in->IsNot(SOLID) && in->IsNot(ISOLATED))
477  {
478  if (mrRemesh.Refine->RemovingOptions.Is(MesherUtilities::REMOVE_NODES_ON_DISTANCE))
479  {
480 
481  if (in->IsNot(FREE_SURFACE) && in->IsNot(RIGID) && freeSurfaceNeighNodes == dimension)
482  {
483  NodeWeakPtrVectorType &neighb_nodes = in->GetValue(NEIGHBOUR_NODES);
484  array_1d<double, 3> sumOfCoordinates = in->Coordinates();
485  array_1d<double, 3> sumOfCurrentVelocities = in->FastGetSolutionStepValue(VELOCITY, 0);
486  array_1d<double, 3> sumOfPreviousVelocities = in->FastGetSolutionStepValue(VELOCITY, 1);
487  double sumOfPressures = in->FastGetSolutionStepValue(PRESSURE, 0);
488  double counter = 1.0;
489 
490  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
491  {
492  counter += 1.0;
493  noalias(sumOfCoordinates) += (nn)->Coordinates();
494  noalias(sumOfCurrentVelocities) += (nn)->FastGetSolutionStepValue(VELOCITY, 0);
495  noalias(sumOfPreviousVelocities) += (nn)->FastGetSolutionStepValue(VELOCITY, 1);
496  sumOfPressures += (nn)->FastGetSolutionStepValue(PRESSURE, 0);
497  }
498  in->X() = sumOfCoordinates[0] / counter;
499  in->Y() = sumOfCoordinates[1] / counter;
500  in->X0() = sumOfCoordinates[0] / counter;
501  in->Y0() = sumOfCoordinates[1] / counter;
502  in->FastGetSolutionStepValue(DISPLACEMENT_X, 0) = 0;
503  in->FastGetSolutionStepValue(DISPLACEMENT_Y, 0) = 0;
504  in->FastGetSolutionStepValue(DISPLACEMENT_X, 1) = 0;
505  in->FastGetSolutionStepValue(DISPLACEMENT_Y, 1) = 0;
506  in->FastGetSolutionStepValue(VELOCITY_X, 0) = sumOfCurrentVelocities[0] / counter;
507  in->FastGetSolutionStepValue(VELOCITY_Y, 0) = sumOfCurrentVelocities[1] / counter;
508  in->FastGetSolutionStepValue(VELOCITY_X, 1) = sumOfPreviousVelocities[0] / counter;
509  in->FastGetSolutionStepValue(VELOCITY_Y, 1) = sumOfPreviousVelocities[1] / counter;
510  in->FastGetSolutionStepValue(PRESSURE, 0) = sumOfPressures / counter;
511  if (dimension == 3)
512  {
513  in->Z() = sumOfCoordinates[2] / counter;
514  in->Z0() = sumOfCoordinates[2] / counter;
515  in->FastGetSolutionStepValue(DISPLACEMENT_Z, 0) = 0;
516  in->FastGetSolutionStepValue(DISPLACEMENT_Z, 1) = 0;
517  in->FastGetSolutionStepValue(VELOCITY_Z, 0) = sumOfCurrentVelocities[2] / counter;
518  in->FastGetSolutionStepValue(VELOCITY_Z, 1) = sumOfPreviousVelocities[2] / counter;
519  }
520  }
521  else
522  {
523  in->Set(TO_ERASE);
524  some_node_is_removed = true;
525  inside_nodes_removed++;
526  if (inletElement)
527  {
528  nodes_removed_inlet_zone++;
529  }
530  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
531  {
532  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
533  }
534  }
535  }
536  }
537  else
538  {
539  SizeType k = 0;
540  SizeType counter = 0;
541  for (std::vector<Node::Pointer>::iterator nn = neighbours.begin(); nn != neighbours.begin() + n_points_in_radius; nn++)
542  {
543  bool nn_on_contact_tip = false;
544 
545  if ((*nn)->Is(TO_SPLIT) || (*nn)->Is(CONTACT))
546  nn_on_contact_tip = true;
547 
548  if ((*nn)->Is(BOUNDARY) && !nn_on_contact_tip && neighbour_distances[k] < size_for_distance_boundary && neighbour_distances[k] > 0.0)
549  {
550  if ((*nn)->IsNot(TO_ERASE))
551  {
552  counter += 1;
553  }
554  }
555 
556  if ((*nn)->Is(BOUNDARY) && nn_on_contact_tip && neighbour_distances[k] < size_for_wall_tip_contact_side)
557  {
558  if ((*nn)->IsNot(TO_ERASE))
559  {
560  counter += 1;
561  }
562  }
563  k++;
564  }
565 
566  if (counter > 1 && in->IsNot(RIGID) && in->IsNot(SOLID) && in->IsNot(NEW_ENTITY) && !on_contact_tip)
567  { // Can be inserted in the boundary refine
568  in->Set(TO_ERASE);
569  if (mEchoLevel > 1)
570  std::cout << " Removed Boundary Node [" << in->Id() << "] on Distance " << std::endl;
571  some_node_is_removed = true;
572  boundary_nodes_removed++;
573 
574  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
575  {
576  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
577  }
578  }
579  }
580  }
581  }
582  }
583 
584  if (boundary_nodes_removed > 0 || inside_nodes_removed>0)
585  {
586  some_node_is_removed = true;
587  }
588  // Build boundary after removing boundary nodes due distance criterion
589  if (mEchoLevel > 1)
590  {
591  std::cout << "boundary_nodes_removed " << boundary_nodes_removed << std::endl;
592  std::cout << "inside_nodes_removed " << inside_nodes_removed << std::endl;
593  }
594  return some_node_is_removed;
595 
596  KRATOS_CATCH(" ")
597  }
598 
599  void EraseCriticalNodes2D(Element::GeometryType &eElement, SizeType &inside_nodes_removed, SizeType &nodes_removed_inlet_zone)
600  {
601 
602  KRATOS_TRY
603 
604  double safetyCoefficient2D = 0.5;
605  double elementVolume = eElement.Area();
606  SizeType numNodes = eElement.size();
607 
608  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
609  const SizeType principalModelPartId = rCurrentProcessInfo[MAIN_MATERIAL_PROPERTY];
610 
611  array_1d<double, 3> Edges(3, 0.0);
612  array_1d<SizeType, 3> FirstEdgeNode(3, 0);
613  array_1d<SizeType, 3> SecondEdgeNode(3, 0);
614  double wallLength = 0;
615  array_1d<double, 3> CoorDifference = eElement[1].Coordinates() - eElement[0].Coordinates();
616  double SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
617  Edges[0] = std::sqrt(SquaredLength);
618  FirstEdgeNode[0] = 0;
619  SecondEdgeNode[0] = 1;
620  if (eElement[0].Is(RIGID) && eElement[1].Is(RIGID))
621  {
622  wallLength = Edges[0];
623  }
624  bool inletElement = false;
625  if (eElement[0].Is(PFEMFlags::EULERIAN_INLET) || eElement[1].Is(PFEMFlags::EULERIAN_INLET) || eElement[2].Is(PFEMFlags::EULERIAN_INLET))
626  {
627  inletElement = true;
628  }
629  SizeType counter = 0;
630  for (SizeType i = 2; i < eElement.size(); i++)
631  {
632  for (SizeType j = 0; j < i; j++)
633  {
634  noalias(CoorDifference) = eElement[i].Coordinates() - eElement[j].Coordinates();
635  SquaredLength = CoorDifference[0] * CoorDifference[0] + CoorDifference[1] * CoorDifference[1];
636  counter += 1;
637  Edges[counter] = std::sqrt(SquaredLength);
638  FirstEdgeNode[counter] = j;
639  SecondEdgeNode[counter] = i;
640  if (eElement[i].Is(RIGID) && eElement[j].Is(RIGID) && Edges[counter] > wallLength)
641  {
642  wallLength = Edges[counter];
643  }
644  }
645  }
646 
648  for (SizeType i = 0; i < eElement.size(); i++)
649  {
650  if (eElement[i].IsNot(RIGID) && eElement[i].IsNot(TO_ERASE) && eElement[i].IsNot(SOLID) && eElement[i].IsNot(ISOLATED))
651  {
652  double height = elementVolume * 2.0 / wallLength;
653 
655  if (eElement[i].Is(FREE_SURFACE))
656  {
657  NodeWeakPtrVectorType &neighb_nodes = eElement[i].GetValue(NEIGHBOUR_NODES);
658  SizeType countRigid = 0;
659  SizeType countFreeSurface = 0;
660  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
661  {
662  if ((nn)->Is(RIGID) || (nn)->Is(SOLID))
663  {
664  countRigid++;
665  }
666  if ((nn)->Is(FREE_SURFACE) && (nn)->IsNot(RIGID) && (nn)->IsNot(SOLID))
667  {
668  countFreeSurface++;
669  }
670  }
671  if ((countRigid + countFreeSurface) == neighb_nodes.size() && countRigid > 0)
672  {
673  safetyCoefficient2D = 0.25;
674  }
675  }
676 
678  if (height < (0.5 * safetyCoefficient2D * wallLength))
679  {
680  eElement[i].Set(TO_ERASE);
681  inside_nodes_removed++;
682 
683  if (inletElement)
684  {
685  nodes_removed_inlet_zone++;
686  }
687 
688  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
689  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
690  {
691  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
692  }
693  }
694  }
695  }
696 
697  bool longDamBreak = false; // to attivate in case of long dam breaks to avoid separeted elements in the water front
698  if (longDamBreak == true)
699  {
700  for (SizeType i = 0; i < numNodes; i++)
701  {
702  if (eElement[i].Is(FREE_SURFACE) && eElement[i].IsNot(RIGID))
703  {
704 
705  GlobalPointersVector<Element> &neighb_elems = eElement[i].GetValue(NEIGHBOUR_ELEMENTS);
706  GlobalPointersVector<Node> &neighb_nodes = eElement[i].GetValue(NEIGHBOUR_NODES);
707 
708  if (neighb_elems.size() < 2)
709  {
710  eElement[i].Set(TO_ERASE);
711  std::cout << "erased an isolated element node" << std::endl;
712  inside_nodes_removed++;
713  if (inletElement)
714  {
715  nodes_removed_inlet_zone++;
716  }
717 
718  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
719  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
720  {
721  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
722  }
723  }
724  else
725  {
726  if (neighb_nodes.size() < 4)
727  {
728  for (SizeType j = 0; j < neighb_nodes.size(); j++)
729  {
730  if (neighb_nodes[j].IsNot(FREE_SURFACE) && neighb_nodes[j].IsNot(RIGID))
731  {
732  break;
733  }
734  if (j == (neighb_nodes.size() - 1))
735  {
736  eElement[i].Set(TO_ERASE);
737  std::cout << "_________________________ erased an isolated element node" << std::endl;
738  inside_nodes_removed++;
739 
740  if (inletElement)
741  {
742  nodes_removed_inlet_zone++;
743  }
744 
745  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
746  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
747  {
748  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
749  }
750  }
751  }
752  }
753  }
754  }
755  }
756  }
757 
758  KRATOS_CATCH("")
759  }
760 
761  void EraseCriticalNodes3D(Element::GeometryType &eElement, SizeType &inside_nodes_removed, SizeType &nodes_removed_inlet_zone, SizeType rigidNodes)
762  {
763 
764  KRATOS_TRY
765  double safetyCoefficient3D = 0.6;
766 
767  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
768  SizeType principalModelPartId = rCurrentProcessInfo[MAIN_MATERIAL_PROPERTY];
769 
770  SizeType freeSurfaceNodes = 0;
771  SizeType numNodes = eElement.size();
772  double elementVolume = eElement.Volume();
773  double criticalVolume = 0.1 * mrRemesh.Refine->MeanVolume;
774 
775  std::vector<array_1d<double, 3>> rigidNodesCoordinates;
776  std::vector<array_1d<double, 3>> rigidNodesNormals;
777  array_1d<double, 3> notRigidNodeCoordinates(3, 0.0);
778  SizeType notRigidNodeId = 0;
779  rigidNodesCoordinates.resize(3);
780  rigidNodesNormals.resize(3);
781  double baricenterX = 0.25 * (eElement[0].X() + eElement[1].X() + eElement[2].X() + eElement[3].X());
782  double baricenterY = 0.25 * (eElement[0].Y() + eElement[1].Y() + eElement[2].Y() + eElement[3].Y());
783  double baricenterZ = 0.25 * (eElement[0].Z() + eElement[1].Z() + eElement[2].Z() + eElement[3].Z());
784 
785  const SizeType numberOfRefiningBoxes = mrRemesh.UseRefiningBox.size();
786  double currentTime = rCurrentProcessInfo[TIME];
787  for (SizeType index = 0; index < numberOfRefiningBoxes; index++)
788  {
789  double initialTime = mrRemesh.RefiningBoxInitialTime[index];
790  double finalTime = mrRemesh.RefiningBoxFinalTime[index];
791  bool refiningBox = mrRemesh.UseRefiningBox[index];
792 
793  if (refiningBox == true && currentTime > initialTime && currentTime < finalTime)
794  {
795  double meanMeshSize = mrRemesh.RefiningBoxMeshSize[index] * 0.5 + mrRemesh.Refine->CriticalRadius * 0.5;
796  array_1d<double, 3> minExternalPoint = mrRemesh.RefiningBoxShiftedMinimumPoint[index];
797  array_1d<double, 3> maxExternalPoint = mrRemesh.RefiningBoxShiftedMaximumPoint[index];
798  array_1d<double, 3> RefiningBoxMinimumPoint = mrRemesh.RefiningBoxMinimumPoint[index];
799  array_1d<double, 3> RefiningBoxMaximumPoint = mrRemesh.RefiningBoxMaximumPoint[index];
800  if (baricenterX > RefiningBoxMinimumPoint[0] && baricenterX < RefiningBoxMaximumPoint[0] &&
801  baricenterY > RefiningBoxMinimumPoint[1] && baricenterY < RefiningBoxMaximumPoint[1] &&
802  baricenterZ > RefiningBoxMinimumPoint[2] && baricenterZ < RefiningBoxMaximumPoint[2])
803  {
804  criticalVolume = 0.01 * (std::pow(mrRemesh.RefiningBoxMeshSize[index], 3) / (6.0 * std::sqrt(2))); // mean Volume of a regular tetrahedral per node with 0.01 of penalization
805  }
806  else if ((baricenterX < RefiningBoxMinimumPoint[0] && baricenterX > minExternalPoint[0]) || (baricenterX > RefiningBoxMaximumPoint[0] && baricenterX < maxExternalPoint[0]) ||
807  (baricenterY < RefiningBoxMinimumPoint[1] && baricenterY > minExternalPoint[1]) || (baricenterY > RefiningBoxMaximumPoint[1] && baricenterY < maxExternalPoint[1]) ||
808  (baricenterZ < RefiningBoxMinimumPoint[2] && baricenterZ > minExternalPoint[2]) || (baricenterZ > RefiningBoxMaximumPoint[2] && baricenterZ < maxExternalPoint[2])) // transition zone
809  {
810  criticalVolume = 0.005 * (std::pow(meanMeshSize, 3) / (6.0 * std::sqrt(2)));
811  safetyCoefficient3D *= 0.5;
812  }
813  else
814  {
815  criticalVolume = 0.01 * (std::pow(meanMeshSize, 3) / (6.0 * std::sqrt(2)));
816  }
817  }
818  }
819 
820  bool inletElement = false;
821  if (eElement[0].Is(PFEMFlags::EULERIAN_INLET) || eElement[1].Is(PFEMFlags::EULERIAN_INLET) || eElement[2].Is(PFEMFlags::EULERIAN_INLET) || eElement[3].Is(PFEMFlags::EULERIAN_INLET))
822  {
823  inletElement = true;
824  }
825 
826  SizeType rigidNode = 0;
827  array_1d<double, 3> WallBaricenter = ZeroVector(3);
828 
829  for (SizeType i = 0; i < numNodes; i++)
830  {
831  if (eElement[i].Is(FREE_SURFACE))
832  {
833  freeSurfaceNodes++;
834  }
835  if (rigidNodes == 3)
836  {
837  if (eElement[i].Is(RIGID))
838  {
839  WallBaricenter += eElement[i].Coordinates() / 3.0;
840  rigidNodesCoordinates[rigidNode] = eElement[i].Coordinates();
841  rigidNodesNormals[rigidNode] = eElement[i].FastGetSolutionStepValue(NORMAL);
842  rigidNode++;
843  }
844  else
845  {
846  notRigidNodeCoordinates = eElement[i].Coordinates();
847  notRigidNodeId = i;
848  }
849  }
850  if (elementVolume < criticalVolume)
851  {
852 
853  if (eElement[i].IsNot(RIGID) && eElement[i].IsNot(SOLID) && eElement[i].IsNot(TO_ERASE))
854  {
855  eElement[i].Set(TO_ERASE);
856  if (mEchoLevel > 1)
857  std::cout << "erase this layer node because it may be potentially dangerous and pass through the solid contour" << std::endl;
858  inside_nodes_removed++;
859 
860  if (inletElement)
861  {
862  nodes_removed_inlet_zone++;
863  }
864 
865  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
866  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
867  {
868  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
869  }
870  }
871  }
872  }
873  double wallNodeDistance = -1;
874  if (rigidNode == 3 && eElement[notRigidNodeId].IsNot(TO_ERASE))
875  {
876 
877  double a1 = 0; // slope x,y,z for the plane composed by rigid nodes only
878  double b1 = 0;
879  double c1 = 0;
880  a1 = (rigidNodesCoordinates[1][1] - rigidNodesCoordinates[0][1]) * (rigidNodesCoordinates[2][2] - rigidNodesCoordinates[0][2]) - (rigidNodesCoordinates[2][1] - rigidNodesCoordinates[0][1]) * (rigidNodesCoordinates[1][2] - rigidNodesCoordinates[0][2]);
881  b1 = (rigidNodesCoordinates[1][2] - rigidNodesCoordinates[0][2]) * (rigidNodesCoordinates[2][0] - rigidNodesCoordinates[0][0]) - (rigidNodesCoordinates[2][2] - rigidNodesCoordinates[0][2]) * (rigidNodesCoordinates[1][0] - rigidNodesCoordinates[0][0]);
882  c1 = (rigidNodesCoordinates[1][0] - rigidNodesCoordinates[0][0]) * (rigidNodesCoordinates[2][1] - rigidNodesCoordinates[0][1]) - (rigidNodesCoordinates[2][0] - rigidNodesCoordinates[0][0]) * (rigidNodesCoordinates[1][1] - rigidNodesCoordinates[0][1]);
883  double a2 = 0; // slope x,y,z between 2 rigid nodes and a not rigid node
884  double b2 = 0;
885  double c2 = 0;
886  a2 = (rigidNodesCoordinates[1][1] - rigidNodesCoordinates[0][1]) * (notRigidNodeCoordinates[2] - rigidNodesCoordinates[0][2]) - (notRigidNodeCoordinates[1] - rigidNodesCoordinates[0][1]) * (rigidNodesCoordinates[1][2] - rigidNodesCoordinates[0][2]);
887  b2 = (rigidNodesCoordinates[1][2] - rigidNodesCoordinates[0][2]) * (notRigidNodeCoordinates[0] - rigidNodesCoordinates[0][0]) - (notRigidNodeCoordinates[2] - rigidNodesCoordinates[0][2]) * (rigidNodesCoordinates[1][0] - rigidNodesCoordinates[0][0]);
888  c2 = (rigidNodesCoordinates[1][0] - rigidNodesCoordinates[0][0]) * (notRigidNodeCoordinates[1] - rigidNodesCoordinates[0][1]) - (notRigidNodeCoordinates[0] - rigidNodesCoordinates[0][0]) * (rigidNodesCoordinates[1][1] - rigidNodesCoordinates[0][1]);
889  double a3 = 0; // slope x,y,z between 2 rigid nodes and a not rigid node
890  double b3 = 0;
891  double c3 = 0;
892  a3 = (rigidNodesCoordinates[1][1] - rigidNodesCoordinates[2][1]) * (notRigidNodeCoordinates[2] - rigidNodesCoordinates[2][2]) - (notRigidNodeCoordinates[1] - rigidNodesCoordinates[2][1]) * (rigidNodesCoordinates[1][2] - rigidNodesCoordinates[2][2]);
893  b3 = (rigidNodesCoordinates[1][2] - rigidNodesCoordinates[2][2]) * (notRigidNodeCoordinates[0] - rigidNodesCoordinates[2][0]) - (notRigidNodeCoordinates[2] - rigidNodesCoordinates[2][2]) * (rigidNodesCoordinates[1][0] - rigidNodesCoordinates[2][0]);
894  c3 = (rigidNodesCoordinates[1][0] - rigidNodesCoordinates[2][0]) * (notRigidNodeCoordinates[1] - rigidNodesCoordinates[2][1]) - (notRigidNodeCoordinates[0] - rigidNodesCoordinates[2][0]) * (rigidNodesCoordinates[1][1] - rigidNodesCoordinates[2][1]);
895  double a4 = 0; // slope x,y,z between 2 rigid nodes and a not rigid node
896  double b4 = 0;
897  double c4 = 0;
898  a4 = (rigidNodesCoordinates[0][1] - rigidNodesCoordinates[2][1]) * (notRigidNodeCoordinates[2] - rigidNodesCoordinates[2][2]) - (notRigidNodeCoordinates[1] - rigidNodesCoordinates[2][1]) * (rigidNodesCoordinates[0][2] - rigidNodesCoordinates[2][2]);
899  b4 = (rigidNodesCoordinates[0][2] - rigidNodesCoordinates[2][2]) * (notRigidNodeCoordinates[0] - rigidNodesCoordinates[2][0]) - (notRigidNodeCoordinates[2] - rigidNodesCoordinates[2][2]) * (rigidNodesCoordinates[0][0] - rigidNodesCoordinates[2][0]);
900  c4 = (rigidNodesCoordinates[0][0] - rigidNodesCoordinates[2][0]) * (notRigidNodeCoordinates[1] - rigidNodesCoordinates[2][1]) - (notRigidNodeCoordinates[0] - rigidNodesCoordinates[2][0]) * (rigidNodesCoordinates[0][1] - rigidNodesCoordinates[2][1]);
901 
902  // angle between the plane composed by rigid nodes only and the other plans. If the angle is small, the particle can pass through the wall
903  double cosAngle12 = (a1 * a2 + b1 * b2 + c1 * c2) / (std::sqrt(std::pow(a1, 2) + std::pow(b1, 2) + std::pow(c1, 2)) * std::sqrt(std::pow(a2, 2) + std::pow(b2, 2) + std::pow(c2, 2)));
904  double cosAngle13 = (a1 * a3 + b1 * b3 + c1 * c3) / (std::sqrt(std::pow(a1, 2) + std::pow(b1, 2) + std::pow(c1, 2)) * std::sqrt(std::pow(a3, 2) + std::pow(b3, 2) + std::pow(c3, 2)));
905  double cosAngle14 = (a1 * a4 + b1 * b4 + c1 * c4) / (std::sqrt(std::pow(a1, 2) + std::pow(b1, 2) + std::pow(c1, 2)) * std::sqrt(std::pow(a4, 2) + std::pow(b4, 2) + std::pow(c4, 2)));
906 
907  // angle between the normals of the rigid nodes. I want to avoid rigid elements at the corner
908  double cosAngleBetweenNormals01 = (rigidNodesNormals[0][0] * rigidNodesNormals[1][0] + rigidNodesNormals[0][1] * rigidNodesNormals[1][1]) /
909  (std::sqrt(std::pow(rigidNodesNormals[0][0], 2) + std::pow(rigidNodesNormals[0][1], 2)) *
910  std::sqrt(std::pow(rigidNodesNormals[1][0], 2) + std::pow(rigidNodesNormals[1][1], 2)));
911  double cosAngleBetweenNormals02 = (rigidNodesNormals[0][0] * rigidNodesNormals[2][0] + rigidNodesNormals[0][1] * rigidNodesNormals[2][1]) /
912  (std::sqrt(std::pow(rigidNodesNormals[0][0], 2) + std::pow(rigidNodesNormals[0][1], 2)) *
913  std::sqrt(std::pow(rigidNodesNormals[2][0], 2) + std::pow(rigidNodesNormals[2][1], 2)));
914  double cosAngleBetweenNormals12 = (rigidNodesNormals[1][0] * rigidNodesNormals[2][0] + rigidNodesNormals[1][1] * rigidNodesNormals[2][1]) /
915  (std::sqrt(std::pow(rigidNodesNormals[1][0], 2) + std::pow(rigidNodesNormals[1][1], 2)) *
916  std::sqrt(std::pow(rigidNodesNormals[2][0], 2) + std::pow(rigidNodesNormals[2][1], 2)));
917 
918  if ((std::abs(cosAngle12) > 0.995 || std::abs(cosAngle13) > 0.995 || std::abs(cosAngle14) > 0.995) && (cosAngleBetweenNormals01 > 0.99 && cosAngleBetweenNormals02 > 0.99 && cosAngleBetweenNormals12 > 0.99))
919  {
920  eElement[notRigidNodeId].Set(TO_ERASE);
921  inside_nodes_removed++;
922 
923  if (inletElement)
924  {
925  nodes_removed_inlet_zone++;
926  }
927 
928  const SizeType propertyIdNode = eElement[notRigidNodeId].FastGetSolutionStepValue(PROPERTY_ID);
929  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
930  {
931  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
932  }
933  }
934 
935  double pwdDistance = 0.0f;
936  for (std::size_t i = 0; i < 3; i++)
937  {
938  pwdDistance += std::pow(eElement[notRigidNodeId].Coordinates()[i] - WallBaricenter[i], 2);
939  }
940  wallNodeDistance = std::sqrt(pwdDistance);
941  }
942 
943  bool longDamBreak = false; // to attivate in case of long dam breaks to avoid separated elements in the water front
944  if (longDamBreak == true && freeSurfaceNodes > 2 && rigidNodes > 1)
945  {
946  for (SizeType i = 0; i < numNodes; i++)
947  {
948  if (eElement[i].Is(FREE_SURFACE) && eElement[i].IsNot(RIGID))
949  {
950 
951  GlobalPointersVector<Element> &neighb_elems = eElement[i].GetValue(NEIGHBOUR_ELEMENTS);
952  GlobalPointersVector<Node> &neighb_nodes = eElement[i].GetValue(NEIGHBOUR_NODES);
953 
954  if (neighb_elems.size() < 2)
955  {
956  eElement[i].Set(TO_ERASE);
957  inside_nodes_removed++;
958 
959  if (inletElement)
960  {
961  nodes_removed_inlet_zone++;
962  }
963 
964  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
965  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
966  {
967  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
968  }
969  }
970  else
971  {
972  if (neighb_nodes.size() < 10)
973  {
974  SizeType freeSurfaceNodesNeigh = 0;
975  for (SizeType j = 0; j < neighb_nodes.size(); j++)
976  {
977  if (neighb_nodes[j].Is(FREE_SURFACE) && neighb_nodes[j].IsNot(RIGID))
978  {
979  freeSurfaceNodesNeigh++;
980  }
981  if (neighb_nodes[j].IsNot(FREE_SURFACE) && neighb_nodes[j].IsNot(RIGID) && neighb_nodes[j].IsNot(TO_ERASE))
982  {
983  break;
984  }
985  if (j == (neighb_nodes.size() - 1) && freeSurfaceNodesNeigh < 2)
986  {
987  eElement[i].Set(TO_ERASE);
988  inside_nodes_removed++;
989 
990  if (inletElement)
991  {
992  nodes_removed_inlet_zone++;
993  }
994 
995  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
996  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
997  {
998  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
999  }
1000  }
1001  }
1002  }
1003  else
1004  {
1005  for (SizeType j = 0; j < neighb_nodes.size(); j++)
1006  {
1007  if (neighb_nodes[j].IsNot(RIGID))
1008  {
1009  break;
1010  }
1011  if (j == (neighb_nodes.size() - 1))
1012  {
1013  eElement[i].Set(TO_ERASE);
1014  inside_nodes_removed++;
1015 
1016  if (inletElement)
1017  {
1018  nodes_removed_inlet_zone++;
1019  }
1020 
1021  const SizeType propertyIdNode = eElement[i].FastGetSolutionStepValue(PROPERTY_ID);
1022  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
1023  {
1024  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
1025  }
1026  }
1027  }
1028  }
1029  }
1030  }
1031  }
1032  }
1033 
1034  array_1d<double, 6> Edges(6, 0.0);
1035  array_1d<SizeType, 6> FirstEdgeNode(6, 0);
1036  array_1d<SizeType, 6> SecondEdgeNode(6, 0);
1037  double wallLength = 0;
1038 
1039  // //////// to compute the length of the wall edge /////////
1040  array_1d<double, 3> CoorDifference = eElement[1].Coordinates() - eElement[0].Coordinates();
1041 
1042  double SquaredLength = CoorDifference[0] * CoorDifference[0] +
1043  CoorDifference[1] * CoorDifference[1] +
1044  CoorDifference[2] * CoorDifference[2];
1045  Edges[0] = std::sqrt(SquaredLength);
1046  double minimumLength = Edges[0] * 10.0;
1047  FirstEdgeNode[0] = 0;
1048  SecondEdgeNode[0] = 1;
1049  int erasableNode = -1;
1050  if (eElement[0].Is(RIGID) && eElement[1].Is(RIGID))
1051  {
1052  wallLength = Edges[0];
1053  }
1054  if ((eElement[0].Is(RIGID) && eElement[1].IsNot(RIGID)) ||
1055  (eElement[1].Is(RIGID) && eElement[0].IsNot(RIGID)))
1056  {
1057  minimumLength = Edges[0];
1058  if (eElement[0].IsNot(RIGID))
1059  {
1060  erasableNode = 0;
1061  }
1062  else
1063  {
1064  erasableNode = 1;
1065  }
1066  }
1067  SizeType counter = 0;
1068  for (SizeType i = 2; i < eElement.size(); i++)
1069  {
1070  for (SizeType j = 0; j < i; j++)
1071  {
1072  noalias(CoorDifference) = eElement[i].Coordinates() - eElement[j].Coordinates();
1073  // CoorDifference = Element[i].Coordinates() - Element[j].Coordinates();
1074  SquaredLength = CoorDifference[0] * CoorDifference[0] +
1075  CoorDifference[1] * CoorDifference[1] +
1076  CoorDifference[2] * CoorDifference[2];
1077  counter += 1;
1078  Edges[counter] = std::sqrt(SquaredLength);
1079  FirstEdgeNode[counter] = j;
1080  SecondEdgeNode[counter] = i;
1081  if (eElement[i].Is(RIGID) && eElement[j].Is(RIGID) && (wallLength == 0 || Edges[counter] > wallLength))
1082  {
1083  wallLength = Edges[counter];
1084  }
1085  if (((eElement[i].Is(RIGID) && eElement[j].IsNot(RIGID)) ||
1086  (eElement[j].Is(RIGID) && eElement[i].IsNot(RIGID))) &&
1087  eElement[i].IsNot(TO_ERASE) && eElement[j].IsNot(TO_ERASE) &&
1088  (Edges[counter] < minimumLength || minimumLength == 0))
1089  {
1090  minimumLength = Edges[counter];
1091  if (eElement[i].IsNot(RIGID))
1092  {
1093  erasableNode = i;
1094  }
1095  else
1096  {
1097  erasableNode = j;
1098  }
1099  }
1100  }
1101  }
1102 
1104  for (SizeType i = 0; i < eElement.size(); i++)
1105  {
1106  if (eElement[i].IsNot(RIGID) && eElement[i].IsNot(TO_ERASE) && eElement[i].IsNot(SOLID) && eElement[i].IsNot(ISOLATED))
1107  {
1109  if (eElement[i].Is(FREE_SURFACE))
1110  {
1111  NodeWeakPtrVectorType &neighb_nodes = eElement[i].GetValue(NEIGHBOUR_NODES);
1112  SizeType countRigid = 0;
1113  SizeType countFreeSurface = 0;
1114  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
1115  {
1116  if ((nn)->Is(RIGID) || (nn)->Is(SOLID))
1117  {
1118  countRigid++;
1119  }
1120  if ((nn)->Is(FREE_SURFACE) && (nn)->IsNot(RIGID) && (nn)->IsNot(SOLID))
1121  {
1122  countFreeSurface++;
1123  }
1124  }
1125  if ((countRigid + countFreeSurface) == neighb_nodes.size() && countRigid > 0)
1126  {
1127  safetyCoefficient3D = 0.25;
1128  }
1129  }
1130  }
1131  }
1132 
1133  if ((minimumLength < (0.35 * wallLength) || (wallNodeDistance < (0.25 * wallLength) && wallNodeDistance > 0)) && erasableNode > -1)
1134  {
1135  if (eElement[erasableNode].IsNot(RIGID) && eElement[erasableNode].IsNot(TO_ERASE) && eElement[erasableNode].IsNot(SOLID) && eElement[erasableNode].IsNot(ISOLATED) && eElement[erasableNode].IsNot(FREE_SURFACE))
1136  {
1137  eElement[erasableNode].Set(TO_ERASE);
1138  inside_nodes_removed++;
1139 
1140  if (inletElement)
1141  {
1142  nodes_removed_inlet_zone++;
1143  }
1144  const SizeType propertyIdNode = eElement[erasableNode].FastGetSolutionStepValue(PROPERTY_ID);
1145  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
1146  {
1147  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
1148  }
1149  }
1150  }
1151 
1152  // //////// to compare the non-wall length to wall edge length /////////
1153  for (SizeType i = 0; i < Edges.size(); i++)
1154  {
1155 
1156  if (mrRemesh.ExecutionOptions.Is(MesherUtilities::REFINE_WALL_CORNER))
1157  {
1158  const ProcessInfo &rCurrentProcessInfo = mrModelPart.GetProcessInfo();
1159  double currentTime = rCurrentProcessInfo[TIME];
1160  double deltaTime = rCurrentProcessInfo[DELTA_TIME];
1161  if (freeSurfaceNodes == 0 && currentTime > 2.0 * deltaTime)
1162  {
1163  safetyCoefficient3D *= 0.7;
1164  }
1165  if (currentTime < 2.0 * deltaTime)
1166  {
1167  safetyCoefficient3D *= 1.05;
1168  }
1169  }
1170  if (((eElement[FirstEdgeNode[i]].Is(RIGID) && eElement[SecondEdgeNode[i]].IsNot(RIGID)) ||
1171  (eElement[SecondEdgeNode[i]].Is(RIGID) && eElement[FirstEdgeNode[i]].IsNot(RIGID))) &&
1172  eElement[FirstEdgeNode[i]].IsNot(TO_ERASE) &&
1173  eElement[SecondEdgeNode[i]].IsNot(TO_ERASE) &&
1174  Edges[i] < safetyCoefficient3D * wallLength)
1175  {
1176  if (eElement[FirstEdgeNode[i]].IsNot(RIGID) && eElement[FirstEdgeNode[i]].IsNot(SOLID) && eElement[FirstEdgeNode[i]].IsNot(TO_ERASE) && eElement[FirstEdgeNode[i]].IsNot(ISOLATED))
1177  {
1178 
1179  eElement[FirstEdgeNode[i]].Set(TO_ERASE);
1180  inside_nodes_removed++;
1181 
1182  if (inletElement)
1183  {
1184  nodes_removed_inlet_zone++;
1185  }
1186 
1187  const SizeType propertyIdNode = eElement[FirstEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1188  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
1189  {
1190  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
1191  }
1192  }
1193  else if (eElement[SecondEdgeNode[i]].IsNot(RIGID) && eElement[SecondEdgeNode[i]].IsNot(SOLID) && eElement[SecondEdgeNode[i]].IsNot(TO_ERASE) && eElement[SecondEdgeNode[i]].IsNot(ISOLATED))
1194  {
1195  eElement[SecondEdgeNode[i]].Set(TO_ERASE);
1196  inside_nodes_removed++;
1197 
1198  if (inletElement)
1199  {
1200  nodes_removed_inlet_zone++;
1201  }
1202 
1203  const SizeType propertyIdNode = eElement[SecondEdgeNode[i]].FastGetSolutionStepValue(PROPERTY_ID);
1204  if (propertyIdNode != principalModelPartId) // this is to conserve the number of nodes of the smaller domain in case of a two-fluid analysis
1205  {
1206  mrRemesh.Info->BalancePrincipalSecondaryPartsNodes += -1;
1207  }
1208  }
1209  }
1210  }
1211 
1212  KRATOS_CATCH("")
1213  }
1214 
1215  bool CheckForMovingLayerNodes(Node &CheckedNode, const double wallLength)
1216  {
1217  KRATOS_TRY
1218  const SizeType dimension = mrModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
1219 
1220  NodeWeakPtrVectorType &neighb_nodes = CheckedNode.GetValue(NEIGHBOUR_NODES);
1221  bool eraseNode = true;
1222  double maxSquaredDistance = 0;
1223  NodeWeakPtrVectorType::iterator j = neighb_nodes.begin();
1224  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
1225  {
1226  if ((nn)->IsNot(RIGID) && (nn)->IsNot(SOLID))
1227  {
1228  // std::cout<<"neigh coordinates: "<<(nn)->X()<<" "<<(nn)->Y()<<std::endl;
1229  array_1d<double, 3> CoorNeighDifference = CheckedNode.Coordinates() - (nn)->Coordinates();
1230  double squaredDistance = CoorNeighDifference[0] * CoorNeighDifference[0] + CoorNeighDifference[1] * CoorNeighDifference[1];
1231  if (dimension == 3)
1232  {
1233  squaredDistance += CoorNeighDifference[2] * CoorNeighDifference[2];
1234  }
1235  if (squaredDistance > maxSquaredDistance)
1236  {
1237  // std::cout<<"(distances: "<<squaredDistance<<" vs "<<maxSquaredDistance<<")"<<std::endl;
1238  maxSquaredDistance = squaredDistance;
1239  j = nn;
1240  }
1241  }
1242  }
1243  // I have looked for the biggest edge for moving there the layer node
1244  double maxNeighDistance = std::sqrt(maxSquaredDistance);
1245  if (maxNeighDistance > wallLength && wallLength > 0)
1246  {
1247  for (NodeWeakPtrVectorType::iterator nn = neighb_nodes.begin(); nn != neighb_nodes.end(); nn++)
1248  {
1249  if (nn == j)
1250  {
1251 
1252  SizeType idMaster = CheckedNode.GetId();
1253  SizeType idSlave = (j)->GetId();
1254  InterpolateFromTwoNodes(idMaster, idMaster, idSlave);
1255  std::vector<double> NewCoordinates(3);
1256  NewCoordinates[0] = (CheckedNode.X() + (j)->X()) * 0.5;
1257  NewCoordinates[1] = (CheckedNode.Y() + (j)->Y()) * 0.5;
1258  CheckedNode.X() = NewCoordinates[0];
1259  CheckedNode.Y() = NewCoordinates[1];
1260  CheckedNode.X0() = NewCoordinates[0];
1261  CheckedNode.Y0() = NewCoordinates[1];
1262  CheckedNode.FastGetSolutionStepValue(DISPLACEMENT_X, 0) = 0;
1263  CheckedNode.FastGetSolutionStepValue(DISPLACEMENT_Y, 0) = 0;
1264  CheckedNode.FastGetSolutionStepValue(DISPLACEMENT_X, 1) = 0;
1265  CheckedNode.FastGetSolutionStepValue(DISPLACEMENT_Y, 1) = 0;
1266  if (dimension == 3)
1267  {
1268  NewCoordinates[2] = (CheckedNode.Z() + (j)->Z()) * 0.5;
1269  CheckedNode.Z() = NewCoordinates[2];
1270  CheckedNode.Z0() = NewCoordinates[2];
1271  CheckedNode.FastGetSolutionStepValue(DISPLACEMENT_Z, 0) = 0;
1272  CheckedNode.FastGetSolutionStepValue(DISPLACEMENT_Z, 1) = 0;
1273  }
1274  // std::cout<<"new coordinates: "<<CheckedNode.X()<<" "<<CheckedNode.Y()<<std::endl;
1275  // std::cout<<"(distances: "<<maxNeighDistance<<" vs "<<wallLength<<")"<<std::endl;
1276  }
1277  eraseNode = false;
1278  }
1279  }
1280  return eraseNode;
1281 
1282  KRATOS_CATCH("")
1283  }
1284 
1285  void InterpolateFromTwoNodes(SizeType idMaster, SizeType idSlave1, SizeType idSlave2)
1286  {
1287 
1288  KRATOS_TRY
1289 
1290  Node::Pointer MasterNode = mrModelPart.pGetNode(idMaster);
1291  Node::Pointer SlaveNode1 = mrModelPart.pGetNode(idSlave1);
1292  Node::Pointer SlaveNode2 = mrModelPart.pGetNode(idSlave2);
1293 
1294  VariablesList &rVariablesList = mrModelPart.GetNodalSolutionStepVariablesList();
1295 
1296  SizeType buffer_size = MasterNode->GetBufferSize();
1297 
1298  for (VariablesList::const_iterator i_variable = rVariablesList.begin(); i_variable != rVariablesList.end(); i_variable++)
1299  {
1300  // std::cout<<" name "<<i_variable->Name()<<std::endl;
1301  // std::cout<<" type "<<typeid(*i_variable).name()<<std::endl;
1302  std::string variable_name = i_variable->Name();
1303  if (KratosComponents<Variable<double>>::Has(variable_name))
1304  {
1305  // std::cout<<"double"<<std::endl;
1306  const Variable<double> &variable = KratosComponents<Variable<double>>::Get(variable_name);
1307  for (SizeType step = 0; step < buffer_size; step++)
1308  {
1309  // getting the data of the solution step
1310  double &node_data = MasterNode->FastGetSolutionStepValue(variable, step);
1311 
1312  double node0_data = SlaveNode1->FastGetSolutionStepValue(variable, step);
1313  double node1_data = SlaveNode2->FastGetSolutionStepValue(variable, step);
1314 
1315  node_data = (0.5 * node0_data + 0.5 * node1_data);
1316  }
1317  }
1318  else if (KratosComponents<Variable<array_1d<double, 3>>>::Has(variable_name))
1319  {
1320  // std::cout<<"array1d"<<std::endl;
1321  const Variable<array_1d<double, 3>> &variable = KratosComponents<Variable<array_1d<double, 3>>>::Get(variable_name);
1322  for (SizeType step = 0; step < buffer_size; step++)
1323  {
1324  // getting the data of the solution step
1325  array_1d<double, 3> &node_data = MasterNode->FastGetSolutionStepValue(variable, step);
1326  const array_1d<double, 3> &node0_data = SlaveNode1->FastGetSolutionStepValue(variable, step);
1327  const array_1d<double, 3> &node1_data = SlaveNode2->FastGetSolutionStepValue(variable, step);
1328  noalias(node_data) = (0.5 * node0_data + 0.5 * node1_data);
1329  }
1330  }
1331  else if (KratosComponents<Variable<int>>::Has(variable_name))
1332  {
1333  // std::cout<<"int"<<std::endl;
1334  // NO INTERPOLATION
1335  }
1336  else if (KratosComponents<Variable<bool>>::Has(variable_name))
1337  {
1338  // std::cout<<"bool"<<std::endl;
1339  // NO INTERPOLATION
1340  }
1341  else if (KratosComponents<Variable<Matrix>>::Has(variable_name))
1342  {
1343  // std::cout<<"Matrix"<<std::endl;
1344  const Variable<Matrix> &variable = KratosComponents<Variable<Matrix>>::Get(variable_name);
1345  for (SizeType step = 0; step < buffer_size; step++)
1346  {
1347  // getting the data of the solution step
1348  Matrix &node_data = MasterNode->FastGetSolutionStepValue(variable, step);
1349 
1350  Matrix &node0_data = SlaveNode1->FastGetSolutionStepValue(variable, step);
1351  Matrix &node1_data = SlaveNode2->FastGetSolutionStepValue(variable, step);
1352 
1353  if (node_data.size1() > 0 && node_data.size2())
1354  {
1355  if (node_data.size1() == node0_data.size1() && node_data.size2() == node0_data.size2() &&
1356  node_data.size1() == node1_data.size1() && node_data.size2() == node1_data.size2())
1357  {
1358  noalias(node_data) = (0.5 * node0_data + 0.5 * node1_data);
1359  }
1360  }
1361  }
1362  }
1363  else if (KratosComponents<Variable<Vector>>::Has(variable_name))
1364  {
1365  // std::cout<<"Vector"<<std::endl;
1366  const Variable<Vector> &variable = KratosComponents<Variable<Vector>>::Get(variable_name);
1367  for (SizeType step = 0; step < buffer_size; step++)
1368  {
1369  // getting the data of the solution step
1370  Vector &node_data = MasterNode->FastGetSolutionStepValue(variable, step);
1371 
1372  Vector &node0_data = SlaveNode1->FastGetSolutionStepValue(variable, step);
1373  Vector &node1_data = SlaveNode2->FastGetSolutionStepValue(variable, step);
1374 
1375  if (node_data.size() > 0)
1376  {
1377  if (node_data.size() == node0_data.size() &&
1378  node_data.size() == node1_data.size())
1379  {
1380  noalias(node_data) = (0.5 * node0_data + 0.5 * node1_data);
1381  // node_data = (0.5*node0_data + 0.5*node1_data);
1382  }
1383  }
1384  }
1385  }
1386  }
1387 
1388  KRATOS_CATCH("")
1389  }
1390 
1392  RemoveMeshNodesForFluidsProcess &operator=(RemoveMeshNodesForFluidsProcess const &rOther);
1393 
1395 
1397  // Process(Process const& rOther);
1398 
1400 
1401  }; // Class Process
1402 
1404 
1407 
1411 
1413  inline std::istream &operator>>(std::istream &rIStream,
1414  RemoveMeshNodesForFluidsProcess &rThis);
1415 
1417  inline std::ostream &operator<<(std::ostream &rOStream,
1418  const RemoveMeshNodesForFluidsProcess &rThis)
1419  {
1420  rThis.PrintInfo(rOStream);
1421  rOStream << std::endl;
1422  rThis.PrintData(rOStream);
1423 
1424  return rOStream;
1425  }
1427 
1428 } // namespace Kratos.
1429 
1430 #endif // KRATOS_REMOVE_MESH_NODES_FOR_FLUIDS_PROCESS_H_INCLUDED defined
Kratos::Bucket
Short class definition.
Definition: bucket.h:57
Kratos::Condition
Base class for all Conditions.
Definition: condition.h:59
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::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::PointsArrayType
PointerVector< TPointType > PointsArrayType
Definition: geometry.h:118
Kratos::GlobalPointersVector
This class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
Kratos::GlobalPointersVector::iterator
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: global_pointers_vector.h:79
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::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::NodesContainerType
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
Kratos::ModelPart::Nodes
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
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::RemoveMeshNodesForFluidsProcess
Remove Mesh Nodes Process for 2D and 3D cases.
Definition: remove_mesh_nodes_for_fluids_process.hpp:55
Kratos::RemoveMeshNodesForFluidsProcess::GeometryType
ConditionType::GeometryType GeometryType
Definition: remove_mesh_nodes_for_fluids_process.hpp:65
Kratos::RemoveMeshNodesForFluidsProcess::RemoveMeshNodesForFluidsProcess
RemoveMeshNodesForFluidsProcess(ModelPart &rModelPart, MesherUtilities::MeshingParameters &rRemeshingParameters, int EchoLevel)
Default constructor.
Definition: remove_mesh_nodes_for_fluids_process.hpp:77
Kratos::RemoveMeshNodesForFluidsProcess::ElementWeakPtrVectorType
GlobalPointersVector< Element > ElementWeakPtrVectorType
Definition: remove_mesh_nodes_for_fluids_process.hpp:70
Kratos::RemoveMeshNodesForFluidsProcess::Info
std::string Info() const override
Turn back information as a string.
Definition: remove_mesh_nodes_for_fluids_process.hpp:167
Kratos::RemoveMeshNodesForFluidsProcess::KdtreeType
Tree< KDTreePartition< BucketType > > KdtreeType
Definition: remove_mesh_nodes_for_fluids_process.hpp:67
Kratos::RemoveMeshNodesForFluidsProcess::PointsArrayType
ModelPart::MeshType::GeometryType::PointsArrayType PointsArrayType
Definition: remove_mesh_nodes_for_fluids_process.hpp:68
Kratos::RemoveMeshNodesForFluidsProcess::SizeType
std::size_t SizeType
Definition: remove_mesh_nodes_for_fluids_process.hpp:71
Kratos::RemoveMeshNodesForFluidsProcess::PrintInfo
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: remove_mesh_nodes_for_fluids_process.hpp:173
Kratos::RemoveMeshNodesForFluidsProcess::~RemoveMeshNodesForFluidsProcess
virtual ~RemoveMeshNodesForFluidsProcess()
Destructor.
Definition: remove_mesh_nodes_for_fluids_process.hpp:89
Kratos::RemoveMeshNodesForFluidsProcess::NodeWeakPtrVectorType
GlobalPointersVector< Node > NodeWeakPtrVectorType
Definition: remove_mesh_nodes_for_fluids_process.hpp:69
Kratos::RemoveMeshNodesForFluidsProcess::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(RemoveMeshNodesForFluidsProcess)
Pointer definition of Process.
Kratos::RemoveMeshNodesForFluidsProcess::ConditionType
ModelPart::ConditionType ConditionType
Definition: remove_mesh_nodes_for_fluids_process.hpp:63
Kratos::RemoveMeshNodesForFluidsProcess::Execute
void Execute() override
Execute method is used to execute the Process algorithms.
Definition: remove_mesh_nodes_for_fluids_process.hpp:106
Kratos::RemoveMeshNodesForFluidsProcess::BucketType
Bucket< 3, Node, std::vector< Node::Pointer >, Node::Pointer, std::vector< Node::Pointer >::iterator, std::vector< double >::iterator > BucketType
Definition: remove_mesh_nodes_for_fluids_process.hpp:66
Kratos::RemoveMeshNodesForFluidsProcess::PrintData
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: remove_mesh_nodes_for_fluids_process.hpp:179
Kratos::RemoveMeshNodesForFluidsProcess::PropertiesType
ModelPart::PropertiesType PropertiesType
Definition: remove_mesh_nodes_for_fluids_process.hpp:64
Kratos::RemoveMeshNodesForFluidsProcess::operator()
void operator()()
This operator is provided to call the process as a function and simply calls the Execute method.
Definition: remove_mesh_nodes_for_fluids_process.hpp:96
Kratos::Tree
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
Kratos::VariablesList::const_iterator
boost::indirect_iterator< VariablesContainerType::const_iterator > const_iterator
Definition: variables_list.h:77
Kratos::array_1d< double, 3 >
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
KRATOS_INFO
#define KRATOS_INFO(label)
Definition: logger.h:250
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
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::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::Vector
Internals::Matrix< double, AMatrix::dynamic, 1 > Vector
Definition: amatrix_interface.h:472
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::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
isotropic_damage_automatic_differentiation.dimension
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
mesh_to_mdpa_converter.radius
float radius
Definition: mesh_to_mdpa_converter.py:18
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
script_THERMAL_CORRECT.counter
int counter
Definition: script_THERMAL_CORRECT.py:218
tensors::i
integer i
Definition: TensorModule.f:17
pfem_fluid_dynamics_application_variables.h
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::RefiningBoxShiftedMaximumPoint
std::vector< array_1d< double, 3 > > RefiningBoxShiftedMaximumPoint
Definition: mesher_utilities.hpp:713
Kratos::MesherUtilities::MeshingParameters::RefiningBoxMeshSize
std::vector< double > RefiningBoxMeshSize
Definition: mesher_utilities.hpp:707
Kratos::MesherUtilities::MeshingParameters::BoundingBoxInitialTime
double BoundingBoxInitialTime
Definition: mesher_utilities.hpp:699
Kratos::MesherUtilities::MeshingParameters::RefiningBoxFinalTime
std::vector< double > RefiningBoxFinalTime
Definition: mesher_utilities.hpp:706
Kratos::MesherUtilities::MeshingParameters::UseBoundingBox
bool UseBoundingBox
Definition: mesher_utilities.hpp:698
Kratos::MesherUtilities::MeshingParameters::BoundingBoxFinalTime
double BoundingBoxFinalTime
Definition: mesher_utilities.hpp:700
Kratos::MesherUtilities::MeshingParameters::UseRefiningBox
std::vector< bool > UseRefiningBox
Definition: mesher_utilities.hpp:704
Kratos::MesherUtilities::MeshingParameters::BoundingBoxUpperPoint
array_1d< double, 3 > BoundingBoxUpperPoint
Definition: mesher_utilities.hpp:702
Kratos::MesherUtilities::MeshingParameters::BoundingBoxLowerPoint
array_1d< double, 3 > BoundingBoxLowerPoint
Definition: mesher_utilities.hpp:701
Kratos::MesherUtilities::MeshingParameters::ExecutionOptions
Flags ExecutionOptions
Definition: mesher_utilities.hpp:648
Kratos::MesherUtilities::MeshingParameters::RefiningBoxMinimumPoint
std::vector< array_1d< double, 3 > > RefiningBoxMinimumPoint
Definition: mesher_utilities.hpp:709
Kratos::MesherUtilities::MeshingParameters::RefiningBoxInitialTime
std::vector< double > RefiningBoxInitialTime
Definition: mesher_utilities.hpp:705
Kratos::MesherUtilities::MeshingParameters::RefiningBoxMaximumPoint
std::vector< array_1d< double, 3 > > RefiningBoxMaximumPoint
Definition: mesher_utilities.hpp:710
Kratos::MesherUtilities::MeshingParameters::RefiningBoxShiftedMinimumPoint
std::vector< array_1d< double, 3 > > RefiningBoxShiftedMinimumPoint
Definition: mesher_utilities.hpp:712
variables_list_data_value_container.h
Generated on Sat Feb 3 2024 02:08:33 for KratosMultiphysics by doxygen 1.9.1