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.
sliding_edge_process.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Klaus B. Sautter
11 //
12 //
13 
14 #ifndef SLIDING_EDGE_PROCESS_H
15 #define SLIDING_EDGE_PROCESS_H
16 
17 // System includes
18 #include <string>
19 #include <iostream>
20 
21 // External includes
22 
23 // Project includes
24 #include "includes/node.h"
25 #include "includes/define.h"
26 #include "processes/process.h"
27 #include "utilities/math_utils.h"
30 #include "includes/model_part.h"
31 
32 namespace Kratos
33 {
34 
35 
37  : public Process
38 {
39  public:
40 
41  typedef Node NodeType;
42  typedef Node ::Pointer NodeTypePointer;
43  typedef std::vector<NodeTypePointer> NodeVector;
45  typedef std::vector<NodeTypePointer>::iterator NodeIterator;
46  typedef std::vector<double> DoubleVector;
47  typedef DoubleVector::iterator DoubleVectorIterator;
48  typedef std::size_t SizeType;
50  typedef ModelPart::MasterSlaveConstraintType::Pointer ConstraintPointer;
51 
52  // Type definitions for tree-search
55 
56 
59 
62  Parameters InputParameters):mrModelPart(rModelPart),mParameters(InputParameters)
63  {
64  KRATOS_TRY;
65  Parameters default_parameters = Parameters(R"(
66  {
67  "constraint_name" : "LinearMasterSlaveConstraint",
68  "master_sub_model_part_name" : "master_connect",
69  "slave_sub_model_part_name" : "slave_connect",
70  "variable_names" : ["DISPLACEMENT_Y","DISPLACEMENT_Z"],
71  "reform_every_step" : true,
72  "debug_info" : true,
73  "angled_initial_line" : false,
74  "follow_line" : false
75  })" );
76  default_parameters.ValidateAndAssignDefaults(InputParameters);
77 
78 
79  KRATOS_ERROR_IF(mParameters["angled_initial_line"].GetBool() && mParameters["follow_line"].GetBool())
80  << "either 'angled_initial_line' or 'follow_line' or both 'false'" << std::endl;
81  KRATOS_CATCH("")
82  }
83 
84 
86  {
87  KRATOS_TRY;
88  if (this->GetmIsInitialized())
89  {
90  if (mParameters["reform_every_step"].GetBool())
91  {
92  //this->CoupleModelParts();
93  this->CoupleModelPartsCustom();
94  }
95  }
96  //else this->CoupleModelParts();
97  else this->CoupleModelPartsCustom();
98  KRATOS_CATCH("");
99  }
100 
102  {
103  KRATOS_TRY;
104  if (mParameters["reform_every_step"].GetBool())
105  {
106  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
107  master_model_part.RemoveMasterSlaveConstraintsFromAllLevels(TO_ERASE);
108  }
109  KRATOS_CATCH("")
110  }
111 
113  {
114  KRATOS_TRY;
115  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
116  ModelPart &slave_model_part = mrModelPart.GetSubModelPart(mParameters["slave_sub_model_part_name"].GetString());
117  const int bucket_size = mParameters["bucket_size"].GetInt();
118  NodesArrayType &r_nodes_master = master_model_part.Nodes();
119  NodesArrayType &r_nodes_slave = slave_model_part.Nodes();
120 
121  NodeVector master_node_list(r_nodes_master.size());
122  this->CreateListOfNodesOfMasterSubModelPart(master_node_list);
123 
124  KDTree::Pointer search_tree =
125  Kratos::shared_ptr<KDTree>(new KDTree(master_node_list.begin(), master_node_list.end(), bucket_size));
126 
127  const int max_number_of_neighbors = 2;
128 
129  for(NodeType& node_i : r_nodes_slave)
130  {
131  double neighbor_search_radius = mParameters["neighbor_search_radius"].GetDouble();
132  SizeType number_of_neighbors = 0;
133  NodeVector neighbor_nodes( max_number_of_neighbors );
134  DoubleVector resulting_squared_distances( max_number_of_neighbors );
135 
136  neighbor_nodes.clear();
137  resulting_squared_distances.clear();
138  //1.) find nodal neighbors
139  number_of_neighbors = search_tree->SearchInRadius( node_i,
140  neighbor_search_radius,
141  neighbor_nodes.begin(),
142  resulting_squared_distances.begin(),
143  max_number_of_neighbors );
144 
145  if (mParameters["must_find_neighbor"].GetBool())
146  {
147  if (number_of_neighbors<1)
148  {
149  KRATOS_ERROR << "found no neighbor for slave node " << node_i.Id() << " " << node_i.Coordinates() << std::endl;
150  }
151  }
152 
153  if (number_of_neighbors>0)
154  {
155  if(mParameters["debug_info"].GetBool()) std::cout << "nr.ne.: " << number_of_neighbors << std::endl;
156  DoubleVector list_of_weights( number_of_neighbors, 0.0 );
157 
158  this->CalculateNodalWeights(resulting_squared_distances,list_of_weights,number_of_neighbors);
159 
160  if(mParameters["angled_initial_line"].GetBool())
161  {this->CoupleSlaveToNeighborMasterNodesInitialAngledLine(node_i,neighbor_nodes,list_of_weights,number_of_neighbors);}
162  else
163  {this->CoupleSlaveToNeighborMasterNodes(node_i,neighbor_nodes,list_of_weights,number_of_neighbors);}
164  }
165 
166  }
167  this->SetmIsInitialized(true);
168  KRATOS_CATCH("");
169  }
170 
171  void CreateListOfNodesOfMasterSubModelPart(NodeVector& rMasterNodeList)
172  {
173  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
174  NodesArrayType &r_nodes = master_model_part.Nodes();
175 
176  rMasterNodeList.resize(r_nodes.size());
177  auto i_begin = master_model_part.NodesBegin();
178 
179  for (SizeType i(0);i<r_nodes.size();++i)
180  {
181  NodeTypePointer pnode = *((i_begin+i).base());
182  rMasterNodeList[i] = pnode;
183  }
184  }
185 
186  void CalculateNodalWeights(const DoubleVector& rResultingSquaredDistances,
187  DoubleVector& rNodalNeighborWeights, const SizeType& rNumberOfNeighbors)
188  {
189  const double numerical_limit = std::numeric_limits<double>::epsilon();
190  double total_nodal_distance = 0.00;
191 
192  if (rNumberOfNeighbors == 1) rNodalNeighborWeights[0] = 1.00;
193  else
194  {
195  if (std::abs(rResultingSquaredDistances[0])<numerical_limit) rNodalNeighborWeights[0] = 1.00;
196  else if (std::abs(rResultingSquaredDistances[1])<numerical_limit) rNodalNeighborWeights[1] = 1.00;
197  else
198  {
199  for (SizeType i=0;i<rNumberOfNeighbors;++i) total_nodal_distance+=std::sqrt(rResultingSquaredDistances[i]);
200  for (SizeType i=0;i<rNumberOfNeighbors;++i)
201  {
202  // change order because smaller distance gets bigger weight!
203  double current_weight = std::sqrt(rResultingSquaredDistances[rNumberOfNeighbors-(i+1)])/total_nodal_distance;
204  (current_weight>numerical_limit) ? (rNodalNeighborWeights[i]=current_weight) : (rNodalNeighborWeights[i]=0.00);
205  }
206  }
207  }
208  }
209 
211  const NodeVector& rNeighborNodes, const DoubleVector& rNodalNeighborWeights,
212  const SizeType& rNumberOfNeighbors)
213  {
214  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
215  ModelPart &slave_model_part = mrModelPart.GetSubModelPart(mParameters["slave_sub_model_part_name"].GetString());
216  NodesArrayType &r_nodes_master = master_model_part.Nodes();
217  NodesArrayType &r_nodes_slave = slave_model_part.Nodes();
218  const double numerical_limit = std::numeric_limits<double>::epsilon();
219 
220 
221  if (rNumberOfNeighbors<2)
222  {
223  for(SizeType dof_iterator=0;dof_iterator<mParameters["variable_names"].size();++dof_iterator)
224  {
225  const auto& current_dof =
226  KratosComponents<VariableType>::Get(mParameters["variable_names"][dof_iterator].GetString());
227 
228  for(SizeType master_iterator =0;master_iterator<rNumberOfNeighbors;++master_iterator)
229  {
230 
232 
233  ConstraintPointer p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(mParameters["constraint_set_name"].GetString(),current_id,
234  r_nodes_master[rNeighborNodes[master_iterator]->Id()],current_dof,r_nodes_slave[rCurrentSlaveNode.Id()],
235  current_dof,rNodalNeighborWeights[master_iterator],0.0);
236 
237  p_current_constraint->Set(TO_ERASE);
238 
239 
240  if(mParameters["debug_info"].GetBool()){
241  std::cout << rNeighborNodes[master_iterator]->Id()
242  << "-----" << rCurrentSlaveNode.Id() << "-----" << rNodalNeighborWeights[master_iterator]
243  << " ::: " << mParameters["variable_names"][dof_iterator].GetString() << std::endl;
244  }
245 
246  } // each master node
247  } // each dof
248  }
249  else if (rNumberOfNeighbors==2)
250  {
251  const auto node_i = rNeighborNodes[0];
252  const auto node_j = rNeighborNodes[1];
253 
254  const auto& check_dof_z = KratosComponents<VariableType>::Get("DISPLACEMENT_Z");
255  const auto& check_dof_y = KratosComponents<VariableType>::Get("DISPLACEMENT_Y");
256  const auto& check_dof_x = KratosComponents<VariableType>::Get("DISPLACEMENT_X");
257 
258  const double dx = (node_j->FastGetSolutionStepValue(DISPLACEMENT_X)+node_j->X0()) -
259  (node_i->FastGetSolutionStepValue(DISPLACEMENT_X)+node_i->X0());
260  const double dy = (node_j->FastGetSolutionStepValue(DISPLACEMENT_Y)+node_j->Y0()) -
261  (node_i->FastGetSolutionStepValue(DISPLACEMENT_Y)+node_i->Y0());
262  const double dz = (node_j->FastGetSolutionStepValue(DISPLACEMENT_Z)+node_j->Z0()) -
263  (node_i->FastGetSolutionStepValue(DISPLACEMENT_Z)+node_i->Z0());
264 
265  // realize the coupling
266  if (std::abs(dx)>numerical_limit)
267  {
269  ConstraintPointer p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(
270  mParameters["constraint_set_name"].GetString(),current_id,
271  r_nodes_slave[rCurrentSlaveNode.Id()],
272  check_dof_x,r_nodes_slave[rCurrentSlaveNode.Id()],
273  check_dof_y,dy/dx,0.0);
274  p_current_constraint->Set(TO_ERASE);
275 
276  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
277  p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(
278  mParameters["constraint_set_name"].GetString(),current_id,
279  r_nodes_slave[rCurrentSlaveNode.Id()],
280  check_dof_x,r_nodes_slave[rCurrentSlaveNode.Id()],
281  check_dof_z,dz/dx,0.0);
282  p_current_constraint->Set(TO_ERASE);
283  }
284  else if (std::abs(dy)>numerical_limit)
285  {
287  ConstraintPointer p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(
288  mParameters["constraint_set_name"].GetString(),current_id,
289  r_nodes_slave[rCurrentSlaveNode.Id()],
290  check_dof_y,r_nodes_slave[rCurrentSlaveNode.Id()],
291  check_dof_x,dx/dy,0.0);
292  p_current_constraint->Set(TO_ERASE);
293 
294  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
295  p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(
296  mParameters["constraint_set_name"].GetString(),current_id,
297  r_nodes_slave[rCurrentSlaveNode.Id()],
298  check_dof_y,r_nodes_slave[rCurrentSlaveNode.Id()],
299  check_dof_z,dz/dy,0.0);
300  p_current_constraint->Set(TO_ERASE);
301  }
302  else if (std::abs(dz)>numerical_limit)
303  {
305  ConstraintPointer p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(
306  mParameters["constraint_set_name"].GetString(),current_id,
307  r_nodes_slave[rCurrentSlaveNode.Id()],
308  check_dof_z,r_nodes_slave[rCurrentSlaveNode.Id()],
309  check_dof_x,dx/dz,0.0);
310  p_current_constraint->Set(TO_ERASE);
311 
312  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
313  p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(
314  mParameters["constraint_set_name"].GetString(),current_id,
315  r_nodes_slave[rCurrentSlaveNode.Id()],
316  check_dof_z,r_nodes_slave[rCurrentSlaveNode.Id()],
317  check_dof_y,dy/dz,0.0);
318  p_current_constraint->Set(TO_ERASE);
319  }
320  else KRATOS_ERROR << "sliding edge is a point for slave node " << rCurrentSlaveNode.Id() << std::endl;
321 
322  }
323  else KRATOS_ERROR << "maximal 2 neighbors allowed for slave node " << rCurrentSlaveNode.Id() << std::endl;
324  }
325 
326  void CoupleSlaveToNeighborMasterNodes(const NodeType& rCurrentSlaveNode,
327  const NodeVector& rNeighborNodes, const DoubleVector& rNodalNeighborWeights,
328  const SizeType& rNumberOfNeighbors)
329  {
330  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
331  ModelPart &slave_model_part = mrModelPart.GetSubModelPart(mParameters["slave_sub_model_part_name"].GetString());
332  NodesArrayType &r_nodes_master = master_model_part.Nodes();
333  NodesArrayType &r_nodes_slave = slave_model_part.Nodes();
334  //const double numerical_limit = std::numeric_limits<double>::epsilon();
335 
336  for(SizeType dof_iterator=0;dof_iterator<mParameters["variable_names"].size();++dof_iterator)
337  {
338  const auto& current_dof =
339  KratosComponents<VariableType>::Get(mParameters["variable_names"][dof_iterator].GetString());
340 
341  for(SizeType master_iterator =0;master_iterator<rNumberOfNeighbors;++master_iterator)
342  {
343 
345 
346  ConstraintPointer p_current_constraint = master_model_part.CreateNewMasterSlaveConstraint(mParameters["constraint_set_name"].GetString(),current_id,
347  r_nodes_master[rNeighborNodes[master_iterator]->Id()],current_dof,r_nodes_slave[rCurrentSlaveNode.Id()],
348  current_dof,rNodalNeighborWeights[master_iterator],0.0);
349 
350  p_current_constraint->Set(TO_ERASE);
351 
352 
353  if(mParameters["debug_info"].GetBool()){
354  std::cout << rNeighborNodes[master_iterator]->Id()
355  << "-----" << rCurrentSlaveNode.Id() << "-----" << rNodalNeighborWeights[master_iterator]
356  << " ::: " << mParameters["variable_names"][dof_iterator].GetString() << std::endl;
357  }
358 
359  } // each master node
360  } // each dof
361  }
362 
363  void SetmIsInitialized(const bool& check) {this->mIsInitialized = check;}
364  bool GetmIsInitialized() const {return this->mIsInitialized;}
365 
366 
367 
368  std::vector<std::size_t> FindNearestNeighboursCustom(const NodeType& node_i, std::vector<double>& weights)
369  {
370  KRATOS_TRY;
371  //const double numerical_limit = std::numeric_limits<double>::epsilon();
372  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
373  NodesArrayType &r_nodes_master = master_model_part.Nodes();
374 
375  double distance = 1e12; // better: std::numeric_limits<double>::max()
376  double distance_i = 0.0;
377  std::vector<std::size_t> neighbour_ids = {0,0};
378  weights.clear();
379  weights = {0.0,0.0};
380 
381  for(NodeType& node_j : r_nodes_master)
382  {
383  distance_i = GetNodalDistance(node_i,node_j);
384  if (distance_i<distance)
385  {
386  distance=distance_i;
387  neighbour_ids[0] = node_j.Id();
388  }
389  }
390  weights[0] = distance;
391  distance = 1e12; // better: std::numeric_limits<double>::max()
392 
393  for(NodeType& node_j : r_nodes_master)
394  {
395  if (node_j.Id()==neighbour_ids[0]) continue;
396 
397  distance_i = GetNodalDistance(node_i,node_j);
398  if (distance_i<distance)
399  {
400  distance=distance_i;
401  neighbour_ids[1] = node_j.Id();
402  }
403  }
404  weights[1] = distance;
405  const double total_distance = weights[0]+weights[1];
406  weights[0] = 1.0 - (weights[0]/total_distance);
407  weights[1] = 1.0 - (weights[1]/total_distance);
408 
409 
410  if (mParameters["debug_info"].GetBool())
411  {
412  KRATOS_WATCH(node_i.Id())
413  KRATOS_WATCH(neighbour_ids)
414  KRATOS_WATCH(weights)
415  std::cout << "_________________________" << std::endl;
416  }
417  return neighbour_ids;
418  KRATOS_CATCH("")
419  }
420 
421  double GetNodalDistance(const NodeType& node_i, const NodeType& node_j)
422  {
423  const array_1d<double, 3> delta_pos =
424  node_j.GetInitialPosition().Coordinates() -
425  node_i.GetInitialPosition().Coordinates() +
426  node_j.FastGetSolutionStepValue(DISPLACEMENT) -
427  node_i.FastGetSolutionStepValue(DISPLACEMENT);
428 
429  const double l = MathUtils<double>::Norm3(delta_pos);
430  return l;
431  }
432 
434  {
435  KRATOS_TRY;
436  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
437  ModelPart &slave_model_part = mrModelPart.GetSubModelPart(mParameters["slave_sub_model_part_name"].GetString());
438  //NodesArrayType &r_nodes_master = master_model_part.Nodes();
439  NodesArrayType &r_nodes_slave = slave_model_part.Nodes();
440 
441  const int max_number_of_neighbours = 2;
442 
443  for(NodeType& node_i : r_nodes_slave)
444  {
445  std::vector<double> nodal_weights = {0.0,0.0};
446  std::vector<std::size_t> neighbour_node_ids = FindNearestNeighboursCustom(node_i,nodal_weights);
447  NodeVector neighbour_nodes( max_number_of_neighbours );
448  neighbour_nodes[0] = master_model_part.pGetNode(neighbour_node_ids[0]);
449  neighbour_nodes[1] = master_model_part.pGetNode(neighbour_node_ids[1]);
450 
451 
452  // couples displacement of slave to itself to follow rigid angled line
453  if(mParameters["angled_initial_line"].GetBool())
454  {this->CoupleSlaveToNeighborMasterNodesInitialAngledLine(node_i,neighbour_nodes,nodal_weights,max_number_of_neighbours);}
455  // linearized line equation
456  else if(mParameters["follow_line"].GetBool())
457  {this->CoupleSlaveToNeighborMasterNodesLineEquation(node_i,neighbour_nodes,nodal_weights,max_number_of_neighbours);}
458  //basic version
459  else
460  {this->CoupleSlaveToNeighborMasterNodes(node_i,neighbour_nodes,nodal_weights,max_number_of_neighbours);}
461  }
462 
463  this->SetmIsInitialized(true);
464  KRATOS_CATCH("");
465  }
466 
467  void CoupleSlaveToNeighborMasterNodesLineEquation(const NodeType& rCurrentSlaveNode,
468  const NodeVector& rNeighborNodes, const DoubleVector& rNodalNeighborWeights,
469  const SizeType& rNumberOfNeighbors)
470  {
471  ModelPart &master_model_part = mrModelPart.GetSubModelPart(mParameters["master_sub_model_part_name"].GetString());
472  ModelPart &slave_model_part = mrModelPart.GetSubModelPart(mParameters["slave_sub_model_part_name"].GetString());
473  NodesArrayType &r_nodes_master = master_model_part.Nodes();
474  NodesArrayType &r_nodes_slave = slave_model_part.Nodes();
475  const double numerical_limit = std::numeric_limits<double>::epsilon();
476 
477 
478  const auto node_i = rNeighborNodes[0];
479  const auto node_j = rNeighborNodes[1];
480 
481  const auto& check_dof_z = KratosComponents<VariableType>::Get("DISPLACEMENT_Z");
482  const auto& check_dof_y = KratosComponents<VariableType>::Get("DISPLACEMENT_Y");
483  const auto& check_dof_x = KratosComponents<VariableType>::Get("DISPLACEMENT_X");
484 
485  const double ua = node_i->FastGetSolutionStepValue(DISPLACEMENT_X);
486  const double va = node_i->FastGetSolutionStepValue(DISPLACEMENT_Y);
487  const double wa = node_i->FastGetSolutionStepValue(DISPLACEMENT_Z);
488 
489  const double ub = node_j->FastGetSolutionStepValue(DISPLACEMENT_X);
490  const double vb = node_j->FastGetSolutionStepValue(DISPLACEMENT_Y);
491  const double wb = node_j->FastGetSolutionStepValue(DISPLACEMENT_Z);
492 
493  /* const double us = rCurrentSlaveNode.FastGetSolutionStepValue(DISPLACEMENT_X);
494  const double vs = rCurrentSlaveNode.FastGetSolutionStepValue(DISPLACEMENT_Y);
495  const double ws = rCurrentSlaveNode.FastGetSolutionStepValue(DISPLACEMENT_Z); */
496 
497  const double xa = node_i->X0();
498  const double ya = node_i->Y0();
499  const double za = node_i->Z0();
500 
501  const double xb = node_j->X0();
502  const double yb = node_j->Y0();
503  const double zb = node_j->Z0();
504 
505  const double xs = rCurrentSlaveNode.X0();
506  const double ys = rCurrentSlaveNode.Y0();
507  const double zs = rCurrentSlaveNode.Z0();
508 
509  const double dx = xb+ub-xa-ua;
510  const double dy = yb+vb-ya-va;
511  const double dz = zb+wb-za-wa;
512 
515  //if explicit time int we might have to couple vel/acc here !
518 
519 
520  // realize the coupling
521  if (std::abs(dx)>numerical_limit)
522  {
523  double constant = (ya-ys) - (xa*(yb-ya)/dx) + (xs*(yb-ya)/dx);
524 
525  double weight_i = 1.0 + (xa/dx) - (xs/dx);
527  ConstraintPointer p_current_constraint_1 = master_model_part.CreateNewMasterSlaveConstraint(
528  mParameters["constraint_set_name"].GetString(),current_id,
529  r_nodes_master[node_i->Id()],
530  check_dof_y,
531  r_nodes_slave[rCurrentSlaveNode.Id()],
532  check_dof_y,
533  weight_i,
534  constant);
535  p_current_constraint_1->Set(TO_ERASE);
536 
537  weight_i = -1.0*dy/dx;
538  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
539  ConstraintPointer p_current_constraint_2 = master_model_part.CreateNewMasterSlaveConstraint(
540  mParameters["constraint_set_name"].GetString(),current_id,
541  r_nodes_master[node_i->Id()],
542  check_dof_x,
543  r_nodes_slave[rCurrentSlaveNode.Id()],
544  check_dof_y,
545  weight_i,
546  0.0);
547  p_current_constraint_2->Set(TO_ERASE);
548 
549 
550  weight_i = (xs/dx) - (xa/dx);
551  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
552  ConstraintPointer p_current_constraint_3 = master_model_part.CreateNewMasterSlaveConstraint(
553  mParameters["constraint_set_name"].GetString(),current_id,
554  r_nodes_master[node_j->Id()],
555  check_dof_y,
556  r_nodes_slave[rCurrentSlaveNode.Id()],
557  check_dof_y,
558  weight_i,
559  0.0);
560  p_current_constraint_3->Set(TO_ERASE);
561 
562 
563  weight_i = 1.0*dy/dx;
564  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
565  ConstraintPointer p_current_constraint_4 = master_model_part.CreateNewMasterSlaveConstraint(
566  mParameters["constraint_set_name"].GetString(),current_id,
567  r_nodes_slave[rCurrentSlaveNode.Id()],
568  check_dof_x,
569  r_nodes_slave[rCurrentSlaveNode.Id()],
570  check_dof_y,
571  weight_i,
572  0.0);
573  p_current_constraint_4->Set(TO_ERASE);
574 
575 
576 
577 
578 
579  constant = (za-zs) - (xa*(zb-za)/dx) + (xs*(zb-za)/dx);
580 
581  weight_i = 1.0 + (xa/dx) - (xs/dx);
582  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
583  ConstraintPointer p_current_constraint_5 = master_model_part.CreateNewMasterSlaveConstraint(
584  mParameters["constraint_set_name"].GetString(),current_id,
585  r_nodes_master[node_i->Id()],
586  check_dof_z,
587  r_nodes_slave[rCurrentSlaveNode.Id()],
588  check_dof_z,
589  weight_i,
590  constant);
591  p_current_constraint_5->Set(TO_ERASE);
592 
593  weight_i = -1.0*dz/dx;
594  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
595  ConstraintPointer p_current_constraint_6 = master_model_part.CreateNewMasterSlaveConstraint(
596  mParameters["constraint_set_name"].GetString(),current_id,
597  r_nodes_master[node_i->Id()],
598  check_dof_x,
599  r_nodes_slave[rCurrentSlaveNode.Id()],
600  check_dof_z,
601  weight_i,
602  0.0);
603  p_current_constraint_6->Set(TO_ERASE);
604 
605 
606  weight_i = (xs/dx) - (xa/dx);
607  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
608  ConstraintPointer p_current_constraint_7 = master_model_part.CreateNewMasterSlaveConstraint(
609  mParameters["constraint_set_name"].GetString(),current_id,
610  r_nodes_master[node_j->Id()],
611  check_dof_z,
612  r_nodes_slave[rCurrentSlaveNode.Id()],
613  check_dof_z,
614  weight_i,
615  0.0);
616  p_current_constraint_7->Set(TO_ERASE);
617 
618 
619  weight_i = 1.0*dz/dx;
620  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
621  ConstraintPointer p_current_constraint_8 = master_model_part.CreateNewMasterSlaveConstraint(
622  mParameters["constraint_set_name"].GetString(),current_id,
623  r_nodes_slave[rCurrentSlaveNode.Id()],
624  check_dof_x,
625  r_nodes_slave[rCurrentSlaveNode.Id()],
626  check_dof_z,
627  weight_i,
628  0.0);
629  p_current_constraint_8->Set(TO_ERASE);
630 
631  }
632  else if (std::abs(dy)>numerical_limit)
633  {
634  double constant = (xa-xs) - (ya*(xb-xa)/dy) + (ys*(xb-xa)/dy);
635 
636  double weight_i = 1.0 + (ya/dy) - (ys/dy);
638  ConstraintPointer p_current_constraint_1 = master_model_part.CreateNewMasterSlaveConstraint(
639  mParameters["constraint_set_name"].GetString(),current_id,
640  r_nodes_master[node_i->Id()],
641  check_dof_x,
642  r_nodes_slave[rCurrentSlaveNode.Id()],
643  check_dof_x,
644  weight_i,
645  constant);
646  p_current_constraint_1->Set(TO_ERASE);
647 
648  weight_i = -1.0*dx/dy;
649  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
650  ConstraintPointer p_current_constraint_2 = master_model_part.CreateNewMasterSlaveConstraint(
651  mParameters["constraint_set_name"].GetString(),current_id,
652  r_nodes_master[node_i->Id()],
653  check_dof_y,
654  r_nodes_slave[rCurrentSlaveNode.Id()],
655  check_dof_x,
656  weight_i,
657  0.0);
658  p_current_constraint_2->Set(TO_ERASE);
659 
660 
661  weight_i = (ys/dy) - (ya/dy);
662  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
663  ConstraintPointer p_current_constraint_3 = master_model_part.CreateNewMasterSlaveConstraint(
664  mParameters["constraint_set_name"].GetString(),current_id,
665  r_nodes_master[node_j->Id()],
666  check_dof_x,
667  r_nodes_slave[rCurrentSlaveNode.Id()],
668  check_dof_x,
669  weight_i,
670  0.0);
671  p_current_constraint_3->Set(TO_ERASE);
672 
673 
674  weight_i = 1.0*dx/dy;
675  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
676  ConstraintPointer p_current_constraint_4 = master_model_part.CreateNewMasterSlaveConstraint(
677  mParameters["constraint_set_name"].GetString(),current_id,
678  r_nodes_slave[rCurrentSlaveNode.Id()],
679  check_dof_y,
680  r_nodes_slave[rCurrentSlaveNode.Id()],
681  check_dof_x,
682  weight_i,
683  0.0);
684  p_current_constraint_4->Set(TO_ERASE);
685 
686 
687  constant = (za-zs) - (ya*(zb-za)/dy) + (ys*(zb-za)/dy);
688 
689  weight_i = 1.0 + (ya/dy) - (ys/dy);
690  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
691  ConstraintPointer p_current_constraint_5 = master_model_part.CreateNewMasterSlaveConstraint(
692  mParameters["constraint_set_name"].GetString(),current_id,
693  r_nodes_master[node_i->Id()],
694  check_dof_z,
695  r_nodes_slave[rCurrentSlaveNode.Id()],
696  check_dof_z,
697  weight_i,
698  constant);
699  p_current_constraint_5->Set(TO_ERASE);
700 
701  weight_i = -1.0*dz/dy;
702  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
703  ConstraintPointer p_current_constraint_6 = master_model_part.CreateNewMasterSlaveConstraint(
704  mParameters["constraint_set_name"].GetString(),current_id,
705  r_nodes_master[node_i->Id()],
706  check_dof_y,
707  r_nodes_slave[rCurrentSlaveNode.Id()],
708  check_dof_z,
709  weight_i,
710  0.0);
711  p_current_constraint_6->Set(TO_ERASE);
712 
713 
714  weight_i = (ys/dy) - (ya/dy);
715  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
716  ConstraintPointer p_current_constraint_7 = master_model_part.CreateNewMasterSlaveConstraint(
717  mParameters["constraint_set_name"].GetString(),current_id,
718  r_nodes_master[node_j->Id()],
719  check_dof_z,
720  r_nodes_slave[rCurrentSlaveNode.Id()],
721  check_dof_z,
722  weight_i,
723  0.0);
724  p_current_constraint_7->Set(TO_ERASE);
725 
726 
727  weight_i = 1.0*dz/dy;
728  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
729  ConstraintPointer p_current_constraint_8 = master_model_part.CreateNewMasterSlaveConstraint(
730  mParameters["constraint_set_name"].GetString(),current_id,
731  r_nodes_slave[rCurrentSlaveNode.Id()],
732  check_dof_y,
733  r_nodes_slave[rCurrentSlaveNode.Id()],
734  check_dof_z,
735  weight_i,
736  0.0);
737  p_current_constraint_8->Set(TO_ERASE);
738  }
739  else if (std::abs(dz)>numerical_limit)
740  {
741  double constant = (xa-xs) - (za*(xb-xa)/dz) + (zs*(xb-xa)/dz);
742 
743  double weight_i = 1.0 + (za/dz) - (zs/dz);
745  ConstraintPointer p_current_constraint_1 = master_model_part.CreateNewMasterSlaveConstraint(
746  mParameters["constraint_set_name"].GetString(),current_id,
747  r_nodes_master[node_i->Id()],
748  check_dof_x,
749  r_nodes_slave[rCurrentSlaveNode.Id()],
750  check_dof_x,
751  weight_i,
752  constant);
753  p_current_constraint_1->Set(TO_ERASE);
754 
755  weight_i = -1.0*dx/dz;
756  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
757  ConstraintPointer p_current_constraint_2 = master_model_part.CreateNewMasterSlaveConstraint(
758  mParameters["constraint_set_name"].GetString(),current_id,
759  r_nodes_master[node_i->Id()],
760  check_dof_z,
761  r_nodes_slave[rCurrentSlaveNode.Id()],
762  check_dof_x,
763  weight_i,
764  0.0);
765  p_current_constraint_2->Set(TO_ERASE);
766 
767 
768  weight_i = (zs/dz) - (za/dz);
769  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
770  ConstraintPointer p_current_constraint_3 = master_model_part.CreateNewMasterSlaveConstraint(
771  mParameters["constraint_set_name"].GetString(),current_id,
772  r_nodes_master[node_j->Id()],
773  check_dof_x,
774  r_nodes_slave[rCurrentSlaveNode.Id()],
775  check_dof_x,
776  weight_i,
777  0.0);
778  p_current_constraint_3->Set(TO_ERASE);
779 
780 
781  weight_i = 1.0*dx/dz;
782  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
783  ConstraintPointer p_current_constraint_4 = master_model_part.CreateNewMasterSlaveConstraint(
784  mParameters["constraint_set_name"].GetString(),current_id,
785  r_nodes_slave[rCurrentSlaveNode.Id()],
786  check_dof_z,
787  r_nodes_slave[rCurrentSlaveNode.Id()],
788  check_dof_x,
789  weight_i,
790  0.0);
791  p_current_constraint_4->Set(TO_ERASE);
792 
793 
794 
795 
796 
797  constant = (ya-ys) - (za*(yb-ya)/dz) + (zs*(yb-ya)/dz);
798 
799  weight_i = 1.0 + (za/dz) - (zs/dz);
800  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
801  ConstraintPointer p_current_constraint_5 = master_model_part.CreateNewMasterSlaveConstraint(
802  mParameters["constraint_set_name"].GetString(),current_id,
803  r_nodes_master[node_i->Id()],
804  check_dof_y,
805  r_nodes_slave[rCurrentSlaveNode.Id()],
806  check_dof_y,
807  weight_i,
808  constant);
809  p_current_constraint_5->Set(TO_ERASE);
810 
811  weight_i = -1.0*dy/dz;
812  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
813  ConstraintPointer p_current_constraint_6 = master_model_part.CreateNewMasterSlaveConstraint(
814  mParameters["constraint_set_name"].GetString(),current_id,
815  r_nodes_master[node_i->Id()],
816  check_dof_z,
817  r_nodes_slave[rCurrentSlaveNode.Id()],
818  check_dof_y,
819  weight_i,
820  0.0);
821  p_current_constraint_6->Set(TO_ERASE);
822 
823 
824  weight_i = (zs/dz) - (za/dz);
825  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
826  ConstraintPointer p_current_constraint_7 = master_model_part.CreateNewMasterSlaveConstraint(
827  mParameters["constraint_set_name"].GetString(),current_id,
828  r_nodes_master[node_j->Id()],
829  check_dof_y,
830  r_nodes_slave[rCurrentSlaveNode.Id()],
831  check_dof_y,
832  weight_i,
833  0.0);
834  p_current_constraint_7->Set(TO_ERASE);
835 
836 
837  weight_i = 1.0*dy/dz;
838  current_id = mrModelPart.NumberOfMasterSlaveConstraints()+1;
839  ConstraintPointer p_current_constraint_8 = master_model_part.CreateNewMasterSlaveConstraint(
840  mParameters["constraint_set_name"].GetString(),current_id,
841  r_nodes_slave[rCurrentSlaveNode.Id()],
842  check_dof_z,
843  r_nodes_slave[rCurrentSlaveNode.Id()],
844  check_dof_y,
845  weight_i,
846  0.0);
847  p_current_constraint_8->Set(TO_ERASE);
848  }
849  else KRATOS_ERROR << "sliding edge is a point for slave node " << rCurrentSlaveNode.Id() << std::endl;
850  }
851 
852 
853 
854  protected:
855 
856  private:
857 
858  ModelPart& mrModelPart;
859  Parameters mParameters;
860  bool mIsInitialized = false;
861 
862 }; // Class
863 
864 }; // namespace
865 
866 #endif
Short class definition.
Definition: bucket.h:57
static const TComponentType & Get(const std::string &rName)
Retrieves a component with the specified name.
Definition: kratos_components.h:114
static double Norm3(const TVectorType &a)
Calculates the norm of vector "a" which is assumed to be of size 3.
Definition: math_utils.h:691
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
void RemoveMasterSlaveConstraintsFromAllLevels(Flags IdentifierFlag=TO_ERASE)
It erases all constraints identified by "IdentifierFlag" by removing the pointer.
Definition: model_part.cpp:1411
std::size_t IndexType
Pointer definition of ModelPart.
Definition: model_part.h:105
SizeType NumberOfMasterSlaveConstraints(IndexType ThisIndex=0) const
Definition: model_part.h:649
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
NodeType::Pointer pGetNode(IndexType NodeId, IndexType ThisIndex=0)
Definition: model_part.h:421
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
MasterSlaveConstraint::Pointer CreateNewMasterSlaveConstraint(const std::string &ConstraintName, IndexType Id, DofsVectorType &rMasterDofsVector, DofsVectorType &rSlaveDofsVector, const MatrixType &RelationMatrix, const VectorType &ConstantVector, IndexType ThisIndex=0)
Creates a new master-slave constraint in the current modelpart.
Definition: model_part.cpp:1234
ModelPart & GetSubModelPart(std::string const &SubModelPartName)
Definition: model_part.cpp:2029
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
This class defines the node.
Definition: node.h:65
TVariableType::Type & FastGetSolutionStepValue(const TVariableType &rThisVariable)
Definition: node.h:435
IndexType Id() const
Definition: node.h:262
double & Z0()
Definition: node.h:576
double & X0()
Definition: node.h:568
double & Y0()
Definition: node.h:572
const PointType & GetInitialPosition() const
Definition: node.h:559
This class provides to Kratos a data structure for I/O based on the standard of JSON.
Definition: kratos_parameters.h:59
double GetDouble() const
This method returns the double contained in the current Parameter.
Definition: kratos_parameters.cpp:657
int GetInt() const
This method returns the integer contained in the current Parameter.
Definition: kratos_parameters.cpp:666
void ValidateAndAssignDefaults(const Parameters &rDefaultParameters)
This function is designed to verify that the parameters under testing match the form prescribed by th...
Definition: kratos_parameters.cpp:1306
SizeType size() const
This method returns the total size of the current array parameter.
Definition: kratos_parameters.cpp:993
std::string GetString() const
This method returns the string contained in the current Parameter.
Definition: kratos_parameters.cpp:684
CoordinatesArrayType const & Coordinates() const
Definition: point.h:215
The base class for all processes in Kratos.
Definition: process.h:49
Definition: sliding_edge_process.h:38
DoubleVector::iterator DoubleVectorIterator
Definition: sliding_edge_process.h:47
ModelPart::MasterSlaveConstraintType::Pointer ConstraintPointer
Definition: sliding_edge_process.h:50
void CoupleSlaveToNeighborMasterNodesLineEquation(const NodeType &rCurrentSlaveNode, const NodeVector &rNeighborNodes, const DoubleVector &rNodalNeighborWeights, const SizeType &rNumberOfNeighbors)
Definition: sliding_edge_process.h:467
void CalculateNodalWeights(const DoubleVector &rResultingSquaredDistances, DoubleVector &rNodalNeighborWeights, const SizeType &rNumberOfNeighbors)
Definition: sliding_edge_process.h:186
void CoupleModelParts()
Definition: sliding_edge_process.h:112
void CoupleModelPartsCustom()
Definition: sliding_edge_process.h:433
double GetNodalDistance(const NodeType &node_i, const NodeType &node_j)
Definition: sliding_edge_process.h:421
void SetmIsInitialized(const bool &check)
Definition: sliding_edge_process.h:363
std::vector< NodeTypePointer > NodeVector
Definition: sliding_edge_process.h:43
std::size_t SizeType
Definition: sliding_edge_process.h:48
void ExecuteInitializeSolutionStep() override
This function will be executed at every time step BEFORE performing the solve phase.
Definition: sliding_edge_process.h:85
bool GetmIsInitialized() const
Definition: sliding_edge_process.h:364
std::vector< std::size_t > FindNearestNeighboursCustom(const NodeType &node_i, std::vector< double > &weights)
Definition: sliding_edge_process.h:368
SlidingEdgeProcess(ModelPart &rModelPart, Parameters InputParameters)
Constructor.
Definition: sliding_edge_process.h:61
void CoupleSlaveToNeighborMasterNodes(const NodeType &rCurrentSlaveNode, const NodeVector &rNeighborNodes, const DoubleVector &rNodalNeighborWeights, const SizeType &rNumberOfNeighbors)
Definition: sliding_edge_process.h:326
Bucket< 3, NodeType, NodeVector, NodeTypePointer, NodeIterator, DoubleVectorIterator > BucketType
Definition: sliding_edge_process.h:53
Tree< KDTreePartition< BucketType > > KDTree
Definition: sliding_edge_process.h:54
void CoupleSlaveToNeighborMasterNodesInitialAngledLine(const NodeType &rCurrentSlaveNode, const NodeVector &rNeighborNodes, const DoubleVector &rNodalNeighborWeights, const SizeType &rNumberOfNeighbors)
Definition: sliding_edge_process.h:210
ModelPart::NodesContainerType NodesArrayType
Definition: sliding_edge_process.h:44
void ExecuteFinalizeSolutionStep() override
This function will be executed at every time step AFTER performing the solve phase.
Definition: sliding_edge_process.h:101
Node ::Pointer NodeTypePointer
Definition: sliding_edge_process.h:42
void CreateListOfNodesOfMasterSubModelPart(NodeVector &rMasterNodeList)
Definition: sliding_edge_process.h:171
Variable< double > VariableType
Definition: sliding_edge_process.h:49
std::vector< double > DoubleVector
Definition: sliding_edge_process.h:46
std::vector< NodeTypePointer >::iterator NodeIterator
Definition: sliding_edge_process.h:45
KRATOS_CLASS_POINTER_DEFINITION(SlidingEdgeProcess)
Pointer definition of ApplyMultipointConstraintsProcess.
Node NodeType
Definition: sliding_edge_process.h:41
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_WATCH(variable)
Definition: define.h:806
#define KRATOS_TRY
Definition: define.h:109
#define KRATOS_ERROR
Definition: exception.h:161
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
std::shared_ptr< T > shared_ptr
Definition: smart_pointers.h:27
int current_id
Output settings end ####.
Definition: script.py:194
xs
Definition: sp_statistics.py:65
ys
Definition: sp_statistics.py:66
integer i
Definition: TensorModule.f:17
integer l
Definition: TensorModule.f:17