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.
trigen_glass_forming.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: Riccardo Rossi
11 //
12 
13 #if !defined(KRATOS_TRIGEN_GLASS_MODELER_H_INCLUDED )
14 #define KRATOS_TRIGEN_GLASS_MODELER_H_INCLUDED
15 
16 // System includes
17 
18 // External includes
19 #include "triangle.h"
20 
21 // Project includes
22 #include "includes/define.h"
23 #include "includes/model_part.h"
26 #include "utilities/timer.h"
29 
30 namespace Kratos
31 {
32 extern "C" {
33  void triangulate(char *, struct triangulateio *, struct triangulateio *,struct triangulateio *);
34  //void trifree();
35 }
36 
37 
40 
44 
48 
52 
56 
58 
61 {
62 public:
65 
68  typedef Node PointType;
69  typedef Node::Pointer PointPointerType;
70  typedef std::vector<PointType::Pointer> PointVector;
71  typedef PointVector::iterator PointIterator;
72  typedef std::vector<double> DistanceVector;
73  typedef std::vector<double>::iterator DistanceIterator;
76 
80 
83  mJ(ZeroMatrix(2,2)), //local jacobian
84  mJinv(ZeroMatrix(2,2)), //inverse jacobian
85  mC(ZeroVector(2)), //dimension = number of nodes
86  mRhs(ZeroVector(2)) {}
87  //mpNodeEraseProcess(NULL){} //dimension = number of nodes
88 
90  virtual ~TriGenGLASSModeler() {}
91 
92 
96 
97 
101 
102 
103  //*******************************************************************************************
104  //*******************************************************************************************
106  ModelPart& ThisModelPart ,
107  Element const& rReferenceElement,
108  Condition const& rReferenceBoundaryCondition,
109  EntitiesEraseProcess<Node>& node_erase, bool rem_nodes = true, bool add_nodes=true,
110  double my_alpha = 1.4, double h_factor=0.5)
111  {
112 
113  KRATOS_TRY
114  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FREE_SURFACE)==false )
115  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_FREE_SURFACE---- variable!!!!!! ERROR","");
116  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_STRUCTURE)==false )
117  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_STRUCTURE---- variable!!!!!! ERROR","");
118  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_BOUNDARY)==false )
119  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_BOUNDARY---- variable!!!!!! ERROR","");
120  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FLUID)==false )
121  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_FLUID---- variable!!!!!! ERROR","");
122 
123  KRATOS_WATCH("Trigen GLASS Refining Mesher")
124  const auto inital_time = std::chrono::steady_clock::now();
125 
126 
127 //clearing elements
128  int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
129 
130  ThisModelPart.Elements().clear();
131  //REFINING EDGES FOR GLASS SIMULATION PROCESSES
132  //NOTE: INBLOW SURFACE IS MARKED WITH FLAG_VARIABLE=1 and OUTER SURFACE IS MARKED WITH IS_INTERFACE=1.. BOTH ARE IS_FREE_SURFACE
133  int id = (ThisModelPart.Nodes().end() - 1)->Id() + 1;
134  for(ModelPart::ConditionsContainerType::iterator ic = ThisModelPart.ConditionsBegin() ; ic != ThisModelPart.ConditionsEnd() ; ic++)
135  {
136 
137  if (ic->GetGeometry().size()==2)
138  {
139  unsigned int n_fs=ic->GetGeometry()[0].FastGetSolutionStepValue(IS_FREE_SURFACE);
140  n_fs+=ic->GetGeometry()[1].FastGetSolutionStepValue(IS_FREE_SURFACE);
141 
142  //THIS REFINES THE NODES OF INTERBAL ELEMENTS OF THE SURFACE WHERE THE INBLOW IS: FLAG_VAR=1
143  if (n_fs==ic->GetGeometry().size())
144  {
145  double x0=ic->GetGeometry()[0].X(); double y0=ic->GetGeometry()[0].Y();
146  double x1=ic->GetGeometry()[1].X(); double y1=ic->GetGeometry()[1].Y();
147 
148 
149  double edge=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1));
150 
151  double nodal_h=ic->GetGeometry()[0].FastGetSolutionStepValue(NODAL_H);
152  nodal_h+=ic->GetGeometry()[1].FastGetSolutionStepValue(NODAL_H);
153 
154  nodal_h*=0.5;
155  //if the edge of the facet (condition) is too long, we split it into two by adding a node in the middle
156 
157 
158  Node::DofsContainerType& reference_dofs = (ThisModelPart.NodesBegin())->GetDofs();
159 
160  double factor=1.25;
161  if (edge>factor*nodal_h)
162  {
163  id++;
164 
165  double x=0.5*(x0+x1); double y=0.5*(y0+y1); double z=0.0;
166  Node::Pointer pnode = ThisModelPart.CreateNewNode(id,x,y,z);
167 
168  //putting the new node also in an auxiliary list
169  //KRATOS_WATCH("adding nodes AT THE EDGE==================================")
170 
171 
172  //std::cout << "new node id = " << pnode->Id() << std::endl;
173  //generating the dofs
174  for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
175  {
176  Node::DofType& rDof = **iii;
177  Node::DofType::Pointer p_new_dof = pnode->pAddDof( rDof );
178 
179  (p_new_dof)->FreeDof();
180  }
181  Geometry<Node >& geom = ic->GetGeometry();
182  InterpolateOnEdge(geom, step_data_size, pnode);
183  const array_1d<double,3>& disp = pnode->FastGetSolutionStepValue(DISPLACEMENT);
184  pnode->X0() = pnode->X() - disp[0];
185  pnode->Y0() = pnode->Y() - disp[1];
186  pnode->Z0() = 0.0;
187  //KRATOS_WATCH("Added node at the EDGE")
188  }
189 
190 
191 
192  }
193  }
194  }
195 
196 
197  ThisModelPart.Conditions().clear();
198 
200 
201 
202  //int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
203 
204  // bucket types
205  //typedef Bucket<3, PointType, ModelPart::NodesContainerType, PointPointerType, PointIterator, DistanceIterator > BucketType;
206  //typedef Bins< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > StaticBins;
207  // bucket types
208 
209  //typedef Tree< StaticBins > Bin; //Binstree;
210  unsigned int bucket_size = 20;
211 
212  //performing the interpolation - all of the nodes in this list will be preserved
213  unsigned int max_results = 100;
214  //PointerVector<PointType> res(max_results);
215  //NodeIterator res(max_results);
216  PointVector res(max_results);
217  DistanceVector res_distances(max_results);
218  Node work_point(0,0.0,0.0,0.0);
219  //if the remove_node switch is activated, we check if the nodes got too close
220 
221  if (rem_nodes==true)
222  {
223  PointVector list_of_nodes;
224  list_of_nodes.reserve(ThisModelPart.Nodes().size());
225  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
226  {
227  (list_of_nodes).push_back(*(i_node.base()));
228  }
229 
230  KdtreeType nodes_tree1(list_of_nodes.begin(),list_of_nodes.end(), bucket_size);
231 
232  RemoveCloseNodes(ThisModelPart, nodes_tree1, node_erase, h_factor);
233  }
234 
236  // //
237  // Now we shall pass the Alpha Shape for the second time, having the "bad nodes" removed //
239  //creating the containers for the input and output
240  struct triangulateio in_mid;
241  struct triangulateio out_mid;
242  struct triangulateio vorout_mid;
243 
244  initialize_triangulateio(in_mid);
245  initialize_triangulateio(out_mid);
246  initialize_triangulateio(vorout_mid);
247 
248  //assigning the size of the input containers
249 
250  in_mid.numberofpoints = ThisModelPart.Nodes().size();
251  in_mid.pointlist = (REAL *) malloc(in_mid.numberofpoints * 2 * sizeof(REAL));
252 
253  //reorder the id's and give the coordinates to the mesher
254  ModelPart::NodesContainerType::iterator nodes_begin = ThisModelPart.NodesBegin();
255  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
256  {
257  int base = i*2;
258  //int base = ((nodes_begin + i)->Id() - 1 ) * 2;
259 
260  //from now on it is consecutive
261  (nodes_begin + i)->SetId(i+1);
262 // (nodes_begin + i)->Id() = i+1;
263 
264  in_mid.pointlist[base] = (nodes_begin + i)->X();
265  in_mid.pointlist[base+1] = (nodes_begin + i)->Y();
266 
267  auto& node_dofs = (nodes_begin + i)->GetDofs();
268  for(auto iii = node_dofs.begin(); iii != node_dofs.end(); iii++)
269  {
270  (**iii).SetEquationId(i+1);
271  }
272  }
273  //in_mid.numberoftriangles = ThisModelPart.Elements().size();
274  //in_mid.trianglelist = (int*) malloc(in_mid.numberoftriangles * 3 * sizeof(int));
275 
276  // "P" suppresses the output .poly file. Saves disk space, but you
277  //lose the ability to maintain constraining segments on later refinements of the mesh.
278  // "B" Suppresses boundary markers in the output .node, .poly, and .edge output files
279  // "n" outputs a list of triangles neighboring each triangle.
280  // "e" outputs edge list (i.e. all the "connectivities")
281  char options1[] = "Pne";
282  triangulate(options1, &in_mid, &out_mid, &vorout_mid);
283  //print out the mesh generation time
284  std::cout << "mesh generation time = " << Timer::ElapsedSeconds(inital_time) << std::endl;
285  //number of newly generated triangles
286  unsigned int el_number=out_mid.numberoftriangles;
287 
288  //PASSING THE CREATED ELEMENTS TO THE ALPHA-SHAPE
289  std::vector<int> preserved_list1(el_number);
290  preserved_list1.resize(el_number);
291 
292  int number_of_preserved_elems= PassAlphaShape(ThisModelPart, preserved_list1, el_number, my_alpha, out_mid);
293 
294  //freeing memory
295 
296  clean_triangulateio(in_mid);
297  clean_triangulateio(vorout_mid);
298  //KRATOS_WATCH("ln367");
299  //NOW WE SHALL PERFORM ADAPTIVE REMESHING, i.e. insert and remove nodes based upon mesh quality
300  // and prescribed sizes
301  struct triangulateio in2;
302  struct triangulateio out2;
303  struct triangulateio vorout2;
304 
305  initialize_triangulateio(in2);
306  initialize_triangulateio(out2);
307  initialize_triangulateio(vorout2);
308 
309 // in2.firstnumber = 1;
310  in2.numberofpoints = ThisModelPart.Nodes().size();
311  in2.pointlist = (REAL *) malloc(in2.numberofpoints * 2 * sizeof(REAL));
312 
313  //writing the input point list
314  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
315  {
316  int base = i*2;
317  in2.pointlist[base] = (nodes_begin + i)->X();
318  in2.pointlist[base+1] = (nodes_begin + i)->Y();
319  }
320  in2.numberoftriangles=number_of_preserved_elems;
321 
322  in2.trianglelist = (int *) malloc(in2.numberoftriangles * 3 * sizeof(int));
323  in2.trianglearealist = (REAL *) malloc(in2.numberoftriangles * sizeof(REAL));
324 // in2.trianglelist = new int[in2.numberoftriangles * 3];
325 // in2.trianglearealist = new double[in2.numberoftriangles];
326 
327  //KRATOS_WATCH(el_number);
328  int counter = 0;
329  //here I will assign a huge number of NODAL_H to the free surface nodes, so that there no nodes will be added
330  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
331  {
332  if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
333  {
334 
335  double& val=i_node->FastGetSolutionStepValue(NODAL_H);
336  val*=2.0;
337  //i_node->FastGetSolutionStepValue(NODAL_H,1)=val;
338  //KRATOS_WATCH("AAAAAAAAAAAAAAAAAAAAAAAA!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
339  }
340  }
341 
342 
343  for (unsigned int el=0; el<el_number; el++)
344  {
345  if( static_cast<bool>(preserved_list1[el]) == true )
346  {
347  //saving the list of ONLY preserved triangles, the ones that passed alpha-shape check
348  int new_base = counter*3;
349  int old_base = el*3;
350  //copying in case it is preserved
351  in2.trianglelist[new_base] = out_mid.trianglelist[old_base];
352  in2.trianglelist[new_base+1] = out_mid.trianglelist[old_base+1];
353  in2.trianglelist[new_base+2] = out_mid.trianglelist[old_base+2];
354 
355  //calculate the prescribed h
356  double prescribed_h = (nodes_begin + out_mid.trianglelist[old_base]-1)->FastGetSolutionStepValue(NODAL_H);
357  prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+1]-1)->FastGetSolutionStepValue(NODAL_H);
358  prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+2]-1)->FastGetSolutionStepValue(NODAL_H);
359  prescribed_h *= 0.3333;
360  //if h is the height of a equilateral triangle, the area is 1/2*h*h
361  in2.trianglearealist[counter] = 0.5*(1.5*prescribed_h*1.5*prescribed_h);
362  counter += 1;
363  }
364 
365  }
366  //now I set back the nodal_h
367  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
368  {
369  if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
370  {
371  double& nodal_h=i_node->FastGetSolutionStepValue(NODAL_H);
372  nodal_h/=2.0;
373  }
374  }
375 
376 
377  clean_triangulateio(out_mid);
378  //KRATOS_WATCH("ln420");
379  //here we generate a new mesh adding/removing nodes, by initializing "q"-quality mesh and "a"-area constraint switches
380  //
381  // MOST IMPORTANT IS "r" switch, that refines previously generated mesh!!!!!!!!!!(that is the one given inside in2)
382  //char mesh_regen_opts[] = "YYJaqrn";
383  //char mesh_regen_opts[] = "YJq1.4arn";
384  if (add_nodes==true)
385  {
386  char mesh_regen_opts[] = "YJq1.4arn";
387  triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
388  KRATOS_WATCH("Adaptive remeshing executed")
389  }
390  else
391  {
392  char mesh_regen_opts[] = "YJrn";
393  triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
394  KRATOS_WATCH("Non-Adaptive remeshing executed")
395  }
396 
397  //and now we shall find out where the new nodes belong to
398  //defintions for spatial search
399  typedef Node PointType;
400  typedef Node::Pointer PointPointerType;
401  typedef std::vector<PointType::Pointer> PointVector;
402  //typedef std::vector<PointType::Pointer>::iterator PointIterator;
403  typedef PointVector::iterator PointIterator;
404  typedef std::vector<double> DistanceVector;
405  typedef std::vector<double>::iterator DistanceIterator;
406  //KRATOS_WATCH("ln449");
407 
409 
410  typedef Tree< KDTreePartition<BucketType> > kd_tree; //Kdtree;
411 
412  //int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
413 
414  //creating an auxiliary list for the new nodes
415  PointVector list_of_new_nodes;
416 
417  //node to get the DOFs from
418  Node::DofsContainerType& reference_dofs = (ThisModelPart.NodesBegin())->GetDofs();
419 
420  double z = 0.0;
421  int n_points_before_refinement = in2.numberofpoints;
422  //if points were added, we add them as nodes to the ModelPart
423  if (out2.numberofpoints > n_points_before_refinement )
424  {
425  for(int i = n_points_before_refinement; i<out2.numberofpoints; i++)
426  {
427  int id=i+1;
428  int base = i*2;
429  double& x= out2.pointlist[base];
430  double& y= out2.pointlist[base+1];
431 
432  Node::Pointer pnode = ThisModelPart.CreateNewNode(id,x,y,z);
433 
434  pnode->SetBufferSize(ThisModelPart.NodesBegin()->GetBufferSize() );
435 
436  list_of_new_nodes.push_back( pnode );
437 
438  //std::cout << "new node id = " << pnode->Id() << std::endl;
439  //generating the dofs
440  for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
441  {
442  Node::DofType& rDof = **iii;
443  Node::DofType::Pointer p_new_dof = pnode->pAddDof( rDof );
444 
445  (p_new_dof)->FreeDof();
446 // (p_new_dof)->EquationId() = -1;
447 
448  }
449 
450  }
451  }
452 
453  std::cout << "During refinement we added " << out2.numberofpoints-n_points_before_refinement<< " nodes " <<std::endl;
454  //unsigned int bucket_size = 20;
455  //performing the interpolation - all of the nodes in this list will be preserved
456  //unsigned int max_results = 50;
457  //PointVector results(max_results);
458  //DistanceVector results_distances(max_results);
460 
462 
463  //int number_of_preserved_elems=0;
464 
465  int point_base;
466  //WHAT ARE THOSE????
467 // Node work_point(0,0.0,0.0,0.0);
468  unsigned int MaximumNumberOfResults = list_of_new_nodes.size();
469  PointVector Results(MaximumNumberOfResults);
470  DistanceVector ResultsDistances(MaximumNumberOfResults);
471 
472  step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
473 
474  if(out2.numberofpoints-n_points_before_refinement > 0) //if we added points
475  {
476 
477  kd_tree nodes_tree2(list_of_new_nodes.begin(),list_of_new_nodes.end(),bucket_size);
478  nodes_begin = ThisModelPart.NodesBegin();
479 
480  for(int el = 0; el< in2.numberoftriangles; el++)
481  {
482  int base = el * 3;
483  //coordinates
484  point_base = (in2.trianglelist[base] - 1)*2;
485  x1[0] = in2.pointlist[point_base];
486  x1[1] = in2.pointlist[point_base+1];
487 
488  point_base = (in2.trianglelist[base+1] - 1)*2;
489  x2[0] = in2.pointlist[point_base];
490  x2[1] = in2.pointlist[point_base+1];
491 
492  point_base = (in2.trianglelist[base+2] - 1)*2;
493  x3[0] = in2.pointlist[point_base];
494  x3[1] = in2.pointlist[point_base+1];
495 
496 
497  //find the center and "radius" of the element
498  double xc, yc, radius;
499  CalculateCenterAndSearchRadius( x1[0], x1[1],
500  x2[0], x2[1],
501  x3[0], x3[1],
502  xc,yc,radius);
503 
504  //find all of the new nodes within the radius
505  int number_of_points_in_radius;
506  work_point.X() = xc;
507  work_point.Y() = yc;
508  work_point.Z() = 0.0;
509 
510  number_of_points_in_radius = nodes_tree2.SearchInRadius(work_point, radius*1.01, Results.begin(),
511  ResultsDistances.begin(), MaximumNumberOfResults);
512 
513  Triangle2D3<Node > geom(
514  *( (nodes_begin + in2.trianglelist[base]-1).base() ),
515  *( (nodes_begin + in2.trianglelist[base+1]-1).base() ),
516  *( (nodes_begin + in2.trianglelist[base+2]-1).base() )
517  );
518 
519  //check if inside and eventually interpolate
520  for( PointIterator it_found = Results.begin(); it_found != Results.begin() + number_of_points_in_radius; it_found++)
521  {
522  bool is_inside = false;
523  is_inside = CalculatePosition(x1[0], x1[1],
524  x2[0], x2[1],
525  x3[0], x3[1],
526  (*it_found)->X(),(*it_found)->Y(),N);
527 
528 
529  if(is_inside == true)
530  {
531  Interpolate( geom, N, step_data_size, *(it_found ) );
532 
533  }
534  }
535  }
536  }
537 
538  ThisModelPart.Elements().clear();
539 
540  //set the coordinates to the original value
541  for( PointVector::iterator it = list_of_new_nodes.begin(); it!=list_of_new_nodes.end(); it++)
542  {
543  const array_1d<double,3>& disp = (*it)->FastGetSolutionStepValue(DISPLACEMENT);
544  (*it)->X0() = (*it)->X() - disp[0];
545  (*it)->Y0() = (*it)->Y() - disp[1];
546  (*it)->Z0() = 0.0;
547  }
548  //cleaning unnecessary data
549  //in2.deinitialize();
550  //in2.initialize();
551  //free( in2.pointlist );
552 
553  //add preserved elements to the kratos
554  Properties::Pointer properties = ThisModelPart.GetMesh().pGetProperties(1);
555  nodes_begin = ThisModelPart.NodesBegin();
556  (ThisModelPart.Elements()).reserve(out2.numberoftriangles);
557 
558  for(int iii = 0; iii< out2.numberoftriangles; iii++)
559  {
560  int id = iii + 1;
561  int base = iii * 3;
562  Triangle2D3<Node > geom(
563  *( (nodes_begin + out2.trianglelist[base]-1).base() ),
564  *( (nodes_begin + out2.trianglelist[base+1]-1).base() ),
565  *( (nodes_begin + out2.trianglelist[base+2]-1).base() )
566  );
567 
568 
569 #ifdef _DEBUG
570  ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
571  if( *(ModelNodes).find( out2.trianglelist[base]).base() == *(ThisModelPart.Nodes().end()).base() )
572  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
573  if( *(ModelNodes).find( out2.trianglelist[base+1]).base() == *(ThisModelPart.Nodes().end()).base() )
574  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
575  if( *(ModelNodes).find( out2.trianglelist[base+2]).base() == *(ThisModelPart.Nodes().end()).base() )
576  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
577 
578 #endif
579 
580  Element::Pointer p_element = rReferenceElement.Create(id, geom, properties);
581  (ThisModelPart.Elements()).push_back(p_element);
582 
583  }
584  ThisModelPart.Elements().Sort();
585 
586  //filling the neighbour list
587  ModelPart::ElementsContainerType::const_iterator el_begin = ThisModelPart.ElementsBegin();
588  for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.ElementsBegin();
589  iii != ThisModelPart.ElementsEnd(); iii++)
590  {
591  //Geometry< Node >& geom = iii->GetGeometry();
592  int base = ( iii->Id() - 1 )*3;
593 
594  (iii->GetValue(NEIGHBOUR_ELEMENTS)).resize(3);
595  GlobalPointersVector< Element >& neighb = iii->GetValue(NEIGHBOUR_ELEMENTS);
596  for(int i = 0; i<3; i++)
597  {
598  int index = out2.neighborlist[base+i];
599  if(index > 0)
600  neighb(i) = GlobalPointer<Element>(&*(el_begin + index-1));
601  else
602  neighb(i) = Element::WeakPointer();
603  }
604  }
605  //identifying boundary, creating skin
606  IdentifyBoundary(ThisModelPart, rReferenceBoundaryCondition, properties, out2);
607 
608  //KRATOS_WATCH("ln749");
609 
610  clean_triangulateio(in2);
611  //KRATOS_WATCH("ln752");
612  clean_triangulateio(out2);
613  //KRATOS_WATCH("ln754");
614  clean_triangulateio(vorout2);
615  KRATOS_WATCH("Finished remeshing with Trigen_GLASS_Refine")
616 
617  KRATOS_CATCH("")
618  }
619 
620 
624 
625 
629 
630 
634 
636  virtual std::string Info() const
637  {
638  return "";
639  }
640 
642  virtual void PrintInfo(std::ostream& rOStream) const {}
643 
645  virtual void PrintData(std::ostream& rOStream) const {}
646 
647 
651 
652 
654 
655 protected:
658 
659 
663  //KdtreeType mKdtree;
664 
668 
669 
673 
674 
678 
679 
683 
684 
688 
689 
691 
692 private:
695 
696 
700  boost::numeric::ublas::bounded_matrix<double,2,2> mJ; //local jacobian
701  boost::numeric::ublas::bounded_matrix<double,2,2> mJinv; //inverse jacobian
702  array_1d<double,2> mC; //center pos
703  array_1d<double,2> mRhs; //center pos
704  //EntitiesEraseProcess<Node>* mpNodeEraseProcess;
705 
706 
710  void RemoveCloseNodes(ModelPart& ThisModelPart, KdtreeType& nodes_tree1, EntitiesEraseProcess<Node>& node_erase, double& h_factor)
711  {
712  //unsigned int bucket_size = 20;
713 
714  //performing the interpolation - all of the nodes in this list will be preserved
715  unsigned int max_results = 100;
716  //PointerVector<PointType> res(max_results);
717  //NodeIterator res(max_results);
718  PointVector res(max_results);
719  DistanceVector res_distances(max_results);
720  Node work_point(0,0.0,0.0,0.0);
721 
722  unsigned int n_points_in_radius;
723  //radius means the distance, closer than which no node shall be allowd. if closer -> mark for erasing
724  double radius;
725 
726  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in != ThisModelPart.NodesEnd(); in++)
727  {
728  radius=h_factor*in->FastGetSolutionStepValue(NODAL_H);
729 
730  work_point[0]=in->X();
731  work_point[1]=in->Y();
732  work_point[2]=in->Z();
733 
734  n_points_in_radius = nodes_tree1.SearchInRadius(work_point, radius, res.begin(),res_distances.begin(), max_results);
735  if (n_points_in_radius>1)
736  {
737  if (in->FastGetSolutionStepValue(IS_BOUNDARY)==0.0 && in->FastGetSolutionStepValue(IS_STRUCTURE)==0.0)
738  {
739  //look if we are already erasing any of the other nodes
740  double erased_nodes = 0;
741  for(PointIterator i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
742  erased_nodes += in->Is(TO_ERASE);
743 
744  if( erased_nodes < 1) //we cancel the node if no other nodes are being erased
745  in->Set(TO_ERASE,true);
746 
747  }
748  //GLASS: WE DISTNIGUISH THE INBLOW BY FLAG=1
749  else if ( (in)->FastGetSolutionStepValue(IS_STRUCTURE)!=1.0 && (in)->FastGetSolutionStepValue(FLAG_VARIABLE)!=1.0) //boundary nodes will be removed if they get REALLY close to another boundary node (0.2 * h_factor)
750  {
751  //here we loop over the neighbouring nodes and if there are nodes
752  //with IS_BOUNDARY=1 which are closer than 0.2*nodal_h from our we remove the node we are considering
753  unsigned int k = 0;
754  unsigned int counter = 0;
755  for(PointIterator i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
756  {
757  if ( (*i)->FastGetSolutionStepValue(IS_BOUNDARY,1)==1.0 && res_distances[k] < 0.2*radius && res_distances[k] > 0.0 )
758  {
759  // KRATOS_WATCH( res_distances[k] );
760  counter += 1;
761  }
762  k++;
763  }
764  if(counter > 0)
765  in->Set(TO_ERASE,true);
766  }
767  }
768 
769  }
770 
771 
772  node_erase.Execute();
773 
774  KRATOS_WATCH("Number of nodes after erasing")
775  KRATOS_WATCH(ThisModelPart.Nodes().size())
776  }
777 
778 
779  int PassAlphaShape(ModelPart& ThisModelPart, std::vector<int>& preserved_list1, unsigned int & el_number, double& my_alpha, struct triangulateio& out_mid)
780  {
781  //NOTE THAT preserved_list1 will be overwritten, only the elements that passed alpha-shaoe check will enter it
782 
783  //prepairing for alpha shape passing : creating necessary arrays
784  //list of preserved elements is created: at max el_number can be preserved (all elements)
785 
786 
787  array_1d<double,3> x1,x2,x3,xc;
788 
789  //int number_of_preserved_elems=0;
790  int number_of_preserved_elems=0;
791  int point_base;
792  //loop for passing alpha shape
793  for(unsigned int el = 0; el< el_number; el++)
794  {
795  int base = el * 3;
796 
797  //coordinates
798  point_base = (out_mid.trianglelist[base] - 1)*2;
799  x1[0] = out_mid.pointlist[point_base];
800  x1[1] = out_mid.pointlist[point_base+1];
801 
802  point_base = (out_mid.trianglelist[base+1] - 1)*2;
803  x2[0] = out_mid.pointlist[point_base];
804  x2[1] = out_mid.pointlist[point_base+1];
805 
806  point_base = (out_mid.trianglelist[base+2] - 1)*2;
807  x3[0] = out_mid.pointlist[point_base];
808  x3[1] = out_mid.pointlist[point_base+1];
809 
810  //here we shall temporarily save the elements and pass them afterwards to the alpha shape
811  Geometry<Node > temp;
812 
813  temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base]).base() ) );
814  temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+1]).base() ) );
815  temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+2]).base() ) );
816 
817  int number_of_structure_nodes = int( temp[0].FastGetSolutionStepValue(IS_STRUCTURE) );
818  number_of_structure_nodes += int( temp[1].FastGetSolutionStepValue(IS_STRUCTURE) );
819  number_of_structure_nodes += int( temp[2].FastGetSolutionStepValue(IS_STRUCTURE) );
820 
821  //check the number of nodes of boundary
822  int nfs = int( temp[0].FastGetSolutionStepValue(IS_FREE_SURFACE) );
823  nfs += int( temp[1].FastGetSolutionStepValue(IS_FREE_SURFACE) );
824  nfs += int( temp[2].FastGetSolutionStepValue(IS_FREE_SURFACE) );
825 
826  //check the number of nodes of boundary
827  int nfluid = int( temp[0].FastGetSolutionStepValue(IS_FLUID) );
828  nfluid += int( temp[1].FastGetSolutionStepValue(IS_FLUID) );
829  nfluid += int( temp[2].FastGetSolutionStepValue(IS_FLUID) );
830 
831  //check the number of nodes of boundary
832  // int nboundary = int( temp[0].FastGetSolutionStepValue(IS_BOUNDARY) );
833  // nboundary += int( temp[1].FastGetSolutionStepValue(IS_BOUNDARY) );
834  // nboundary += int( temp[2].FastGetSolutionStepValue(IS_BOUNDARY) );
835 
836  //check the number of nodes of interface - this is how we mark the outer surface in the GLASS FORMING MODELING
837  // int n_int = int( temp[0].FastGetSolutionStepValue(IS_INTERFACE) );
838  // n_int += int( temp[1].FastGetSolutionStepValue(IS_INTERFACE) );
839  // n_int += int( temp[2].FastGetSolutionStepValue(IS_INTERFACE) );
840  //first check that we are working with fluid elements, otherwise throw an error
841  //if (nfluid!=3)
842  // KRATOS_THROW_ERROR(std::logic_error,"THATS NOT FLUID or NOT TRIANGLE!!!!!! ERROR","");
843  //otherwise perform alpha shape check
844 
845 
846  if(number_of_structure_nodes!=3) //if it is = 3 it is a completely fixed element -> do not add it
847  {
848 
849 // if ((nfs != 0 || nfluid != 3) && n_int==0) //in this case it is close to the surface so i should use alpha shape
850  if (nfs != 0 || nfluid != 3) //in this case it is close to the surface so i should use alpha shape
851  {
852 
853  if( AlphaShape(my_alpha,temp)) //if alpha shape says to preserve
854  {
855 // if(nboundary==3 && number_of_structure_nodes > 1 && nfs > 0) //if it is = 3 pressure problems -> do not add it
856 // {
857 // preserved_list1[el] = false;
858 // }
859 // else
860 // {
861  preserved_list1[el] = true;
862  number_of_preserved_elems += 1;
863 // }
864  }
865  }
866  //FOR GLASS FORMING: we keep the contact elements as much as possible
867  /*
868  else if (n_int!=0 && number_of_structure_nodes!=0)
869  {
870  double bigger_alpha = my_alpha*2.0;
871  if( AlphaShape(bigger_alpha,temp) && number_of_structure_nodes!=3)
872  {
873  //KRATOS_WATCH ("PRESRVING INTERFACE ELEMENT")
874  preserved_list1[el] = true;
875  number_of_preserved_elems += 1;
876 //
877  }
878 
879 
880  }
881  */
882  else //internal triangle --- should be ALWAYS preserved
883  {
884  double bigger_alpha = my_alpha*10.0;
885  if( AlphaShape(bigger_alpha,temp) && number_of_structure_nodes!=3)
886  {
887 
888  preserved_list1[el] = true;
889  number_of_preserved_elems += 1;
890 //
891  }
892  }
893  }
894 
895  else
896  preserved_list1[el] = false;
897 
898  }
899  return number_of_preserved_elems;
900  }
901 
902 
903 
904 
905 
906  void IdentifyBoundary(ModelPart& ThisModelPart, Condition const& rReferenceBoundaryCondition, Properties::Pointer& properties, struct triangulateio& out2)
907  {
908 
909  //reset the boundary flag
910  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in!=ThisModelPart.NodesEnd(); in++)
911  {
912  in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
913  }
914  //filling the elemental neighbours list (from now on the elements list can not change)
915  ModelPart::ElementsContainerType::iterator elements_end = ThisModelPart.Elements().end();
916 
917  ThisModelPart.Elements().Unique();
918 
919  //now the boundary faces
920  for(ModelPart::ElementsContainerType::iterator iii = ThisModelPart.ElementsBegin(); iii != ThisModelPart.ElementsEnd(); iii++)
921  {
922  int base = ( iii->Id() - 1 )*3;
923 
924  ModelPart::ElementsContainerType::iterator el_neighb;
925  /*each face is opposite to the corresponding node number so
926  0 ----- 1 2
927  1 ----- 2 0
928  2 ----- 0 1
929  */
930 
932  //
933  //********************************************************************
934  //first face
935  el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base] );
936  if( el_neighb == elements_end )
937  {
938  //std::cout << "node0" << std::endl;
939  //if no neighnour is present => the face is free surface
940  iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
941  iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
942 
943  //Generate condition
945  temp1.reserve(2);
946  temp1.push_back(iii->GetGeometry()(1));
947  temp1.push_back(iii->GetGeometry()(2));
948 
949  Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
950  int id = (iii->Id()-1)*3;
951  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
952  ThisModelPart.Conditions().push_back(p_cond);
953 
954  }
955 
956  //********************************************************************
957  //second face
958  el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+1] );
959  //if( el != ThisModelPart.Elements().end() )
960  if( el_neighb == elements_end )
961  {
962  //if no neighnour is present => the face is free surface
963  iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
964  iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
965 
966  //Generate condition
968  temp1.reserve(2);
969  temp1.push_back(iii->GetGeometry()(2));
970  temp1.push_back(iii->GetGeometry()(0));
971 
972  Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
973  int id = (iii->Id()-1)*3+1;
974  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
975  ThisModelPart.Conditions().push_back(p_cond);
976 
977 
978  }
979 
980  //********************************************************************
981  //third face
982  el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+2] );
983  if( el_neighb == elements_end )
984  {
985  //if no neighnour is present => the face is free surface
986  iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
987  iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
988 
989 // Generate condition
991  temp1.reserve(2);
992  temp1.push_back(iii->GetGeometry()(0));
993  temp1.push_back(iii->GetGeometry()(1));
994 
995  Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
996  int id = (iii->Id()-1)*3+2;
997  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
998  ThisModelPart.Conditions().push_back(p_cond);
999 
1000  }
1001 
1002 
1003  }
1004  }
1005 
1006 
1007  //returns false if it should be removed
1008  bool AlphaShape(double alpha_param, Geometry<Node >& pgeom)
1009  //bool AlphaShape(double alpha_param, Triangle2D<Node >& pgeom)
1010  {
1011  KRATOS_TRY
1012 
1013 
1014  double x0 = pgeom[0].X();
1015  double x1 = pgeom[1].X();
1016  double x2 = pgeom[2].X();
1017 
1018  double y0 = pgeom[0].Y();
1019  double y1 = pgeom[1].Y();
1020  double y2 = pgeom[2].Y();
1021 
1022  mJ(0,0)=2.0*(x1-x0);
1023  mJ(0,1)=2.0*(y1-y0);
1024  mJ(1,0)=2.0*(x2-x0);
1025  mJ(1,1)=2.0*(y2-y0);
1026 
1027 
1028  double detJ = mJ(0,0)*mJ(1,1)-mJ(0,1)*mJ(1,0);
1029 
1030  mJinv(0,0) = mJ(1,1);
1031  mJinv(0,1) = -mJ(0,1);
1032  mJinv(1,0) = -mJ(1,0);
1033  mJinv(1,1) = mJ(0,0);
1034 
1035  bounded_matrix<double,2,2> check;
1036 
1037 //calculate average h
1038  double h;
1039  h = pgeom[0].FastGetSolutionStepValue(NODAL_H);
1040  h += pgeom[1].FastGetSolutionStepValue(NODAL_H);
1041  h += pgeom[2].FastGetSolutionStepValue(NODAL_H);
1042  h *= 0.333333333;
1043 
1044 
1045  if(detJ < 5e-3*h*h)
1046  {
1047  //std::cout << "detJ = " << detJ << std::endl;
1049  pgeom[0].GetSolutionStepValue(IS_BOUNDARY) = 1;
1050  pgeom[1].GetSolutionStepValue(IS_BOUNDARY) = 1;
1051  pgeom[2].GetSolutionStepValue(IS_BOUNDARY) = 1;
1052  return false;
1053  }
1054 
1055  else
1056  {
1057 
1058  double x0_2 = x0*x0 + y0*y0;
1059  double x1_2 = x1*x1 + y1*y1;
1060  double x2_2 = x2*x2 + y2*y2;
1061 
1062  //finalizing the calculation of the inverted matrix
1063  //std::cout<<"MATR INV"<<MatrixInverse(mJ)<<std::endl;
1064  mJinv /= detJ;
1065  //calculating the RHS
1066  mRhs[0] = (x1_2 - x0_2);
1067  mRhs[1] = (x2_2 - x0_2);
1068 
1069  //calculate position of the center
1070  noalias(mC) = prod(mJinv,mRhs);
1071 
1072  double radius = sqrt(pow(mC[0]-x0,2)+pow(mC[1]-y0,2));
1073 
1074 
1075  if (radius < h*alpha_param)
1076  {
1077  return true;
1078  }
1079  else
1080  {
1081  return false;
1082  }
1083  }
1084 
1085 
1086  KRATOS_CATCH("")
1087  }
1088  //AUXILLIARY FCTNS
1089  inline void CalculateCenterAndSearchRadius(const double x0, const double y0,
1090  const double x1, const double y1,
1091  const double x2, const double y2,
1092  double& xc, double& yc, double& R
1093  )
1094  {
1095  xc = 0.3333333333333333333*(x0+x1+x2);
1096  yc = 0.3333333333333333333*(y0+y1+y2);
1097 
1098  double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0);
1099  double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1);
1100  double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2);
1101 
1102  R = R1;
1103  if(R2 > R) R = R2;
1104  if(R3 > R) R = R3;
1105 
1106  R = sqrt(R);
1107  }
1108 
1109  inline double CalculateVol( const double x0, const double y0,
1110  const double x1, const double y1,
1111  const double x2, const double y2
1112  )
1113  {
1114  return 0.5*( (x1-x0)*(y2-y0)- (y1-y0)*(x2-x0) );
1115  }
1116 
1117  inline bool CalculatePosition( const double x0, const double y0,
1118  const double x1, const double y1,
1119  const double x2, const double y2,
1120  const double xc, const double yc,
1121  array_1d<double,3>& N
1122  )
1123  {
1124  double area = CalculateVol(x0,y0,x1,y1,x2,y2);
1125 
1126  if(area < 0.000000000001)
1127  {
1128  KRATOS_THROW_ERROR(std::logic_error,"element with zero area found","");
1129  }
1130 
1131 
1132 
1133  N[0] = CalculateVol(x1,y1,x2,y2,xc,yc) / area;
1134  N[1] = CalculateVol(x2,y2,x0,y0,xc,yc) / area;
1135  N[2] = CalculateVol(x0,y0,x1,y1,xc,yc) / area;
1136 
1137  /* N[0] = CalculateVol(x0,y0,x1,y1,xc,yc) * inv_area;
1138  N[1] = CalculateVol(x1,y1,x2,y2,xc,yc) * inv_area;
1139  N[2] = CalculateVol(x2,y2,x0,y0,xc,yc) * inv_area;*/
1140  double tol = 1e-4;
1141  double upper_limit = 1.0+tol;
1142  double lower_limit = -tol;
1143 
1144  if(N[0] >= lower_limit && N[1] >= lower_limit && N[2] >= lower_limit && N[0] <= upper_limit && N[1] <= upper_limit && N[2] <= upper_limit) //if the xc yc is inside the triangle
1145  //if(N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0) //if the xc yc is inside the triangle return true
1146  return true;
1147 
1148  return false;
1149  }
1150 
1151  void Interpolate( Triangle2D3<Node >& geom, const array_1d<double,3>& N,
1152  unsigned int step_data_size,
1153  Node::Pointer pnode)
1154  {
1155  unsigned int buffer_size = pnode->GetBufferSize();
1156  //KRATOS_WATCH("Buffer size")
1157  //KRATOS_WATCH(buffer_size)
1158 
1159  for(unsigned int step = 0; step<buffer_size; step++)
1160  {
1161 
1162  //getting the data of the solution step
1163  double* step_data = (pnode)->SolutionStepData().Data(step);
1164 
1165 
1166  double* node0_data = geom[0].SolutionStepData().Data(step);
1167  double* node1_data = geom[1].SolutionStepData().Data(step);
1168  double* node2_data = geom[2].SolutionStepData().Data(step);
1169 
1170  //copying this data in the position of the vector we are interested in
1171  for(unsigned int j= 0; j<step_data_size; j++)
1172  {
1173 
1174  step_data[j] = N[0]*node0_data[j] + N[1]*node1_data[j] + N[2]*node2_data[j];
1175 
1176 
1177  }
1178  }
1179  if (N[0]==0.0 && N[1]==0.0 && N[2]==0.0)
1180  KRATOS_THROW_ERROR(std::logic_error,"SOMETHING's wrong with the added nodes!!!!!! ERROR","");
1181 
1182  //if ( pnode->FastGetSolutionStepValue(BULK_MODULUS)==0.0)
1183  // KRATOS_THROW_ERROR(std::logic_error,"SOMETHING's wrong with the added nodes!!!!!! ERROR","");
1184 
1185  //now we assure that the flag variables are set coorect!! since we add nodes inside of the fluid volume only
1186  //we manually reset the IS_BOUNDARY, IS_FLUID, IS_STRUCTURE, IS_FREE_SURFACE values in a right way
1187  //not to have values, like 0.33 0.66 resulting if we would have been interpolating them in the same way
1188  //as the normal variables, like Velocity etc
1189 
1190 
1191  pnode->FastGetSolutionStepValue(IS_BOUNDARY)=0.0;
1192  pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1193  pnode->FastGetSolutionStepValue(IS_INTERFACE)=0.0;
1194  pnode->Set(TO_ERASE,false);
1195  pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=0.0;
1196  pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1197  }
1198 
1199  void InterpolateOnEdge( Geometry<Node >& geom, unsigned int step_data_size, Node::Pointer pnode)
1200  {
1201  unsigned int buffer_size = pnode->GetBufferSize();
1202  //KRATOS_WATCH(buffer_size)
1203 
1204  for(unsigned int step = 0; step<buffer_size; step++)
1205  {
1206 
1207  //getting the data of the solution step
1208  double* step_data = (pnode)->SolutionStepData().Data(step);
1209 
1210 
1211  double* node0_data = geom[0].SolutionStepData().Data(step);
1212  double* node1_data = geom[1].SolutionStepData().Data(step);
1213 
1214  //copying this data in the position of the vector we are interested in
1215  for(unsigned int j= 0; j<step_data_size; j++)
1216  {
1217 
1218  step_data[j] = 0.5*(node0_data[j] + node1_data[j]);
1219 
1220 
1221  }
1222  }
1223  //now we assure that the flag variables are set coorect!! since we add nodes inside of the fluid volume only
1224  //we manually reset the IS_BOUNDARY, IS_FLUID, IS_STRUCTURE, IS_FREE_SURFACE values in a right way
1225  //not to have values, like 0.33 0.66 resulting if we would have been interpolating them in the same way
1226  //as the normal variables, like Velocity etc
1227 
1228  //pnode->FastGetSolutionStepValue(IS_BOUNDARY)=1.0;
1229  //pnode->FastGetSolutionStepValue(FLAG_VARIABLE)=1.0;
1230 
1231  pnode->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET)=0.0;
1232  pnode->FastGetSolutionStepValue(IS_BOUNDARY)=1.0;
1233  pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=1.0;
1234  pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1235  pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1236 
1237 
1238  pnode->Set(TO_ERASE,false);
1239 
1240  if (pnode->FastGetSolutionStepValue(IS_INTERFACE)>0.99)
1241  pnode->FastGetSolutionStepValue(IS_INTERFACE)=1.0;
1242  else
1243  pnode->FastGetSolutionStepValue(IS_INTERFACE)=0.0;
1244 
1245 
1246 
1247  }
1248 
1249 
1250 
1251  void initialize_triangulateio( triangulateio& tr )
1252  {
1253  tr.pointlist = (REAL*) NULL;
1254  tr.pointattributelist = (REAL*) NULL;
1255  tr.pointmarkerlist = (int*) NULL;
1256  tr.numberofpoints = 0;
1257  tr.numberofpointattributes = 0;
1258  tr.trianglelist = (int*) NULL;
1259  tr.triangleattributelist = (REAL*) NULL;
1260  tr.trianglearealist = (REAL*) NULL;
1261  tr.neighborlist = (int*) NULL;
1262  tr.numberoftriangles = 0;
1263  tr.numberofcorners = 3;
1264  tr.numberoftriangleattributes = 0;
1265  tr.segmentlist = (int*) NULL;
1266  tr.segmentmarkerlist = (int*) NULL;
1267  tr.numberofsegments = 0;
1268  tr.holelist = (REAL*) NULL;
1269  tr.numberofholes = 0;
1270  tr.regionlist = (REAL*) NULL;
1271  tr.numberofregions = 0;
1272  tr.edgelist = (int*) NULL;
1273  tr.edgemarkerlist = (int*) NULL;
1274  tr.normlist = (REAL*) NULL;
1275  tr.numberofedges = 0;
1276  };
1277 
1278  void clean_triangulateio( triangulateio& tr )
1279  {
1280  if(tr.pointlist != NULL) free(tr.pointlist );
1281  if(tr.pointattributelist != NULL) free(tr.pointattributelist );
1282  if(tr.pointmarkerlist != NULL) free(tr.pointmarkerlist );
1283  if(tr.trianglelist != NULL) free(tr.trianglelist );
1284  if(tr.triangleattributelist != NULL) free(tr.triangleattributelist );
1285  if(tr.trianglearealist != NULL) free(tr.trianglearealist );
1286  if(tr.neighborlist != NULL) free(tr.neighborlist );
1287  if(tr.segmentlist != NULL) free(tr.segmentlist );
1288  if(tr.segmentmarkerlist != NULL) free(tr.segmentmarkerlist );
1289  if(tr.holelist != NULL) free(tr.holelist );
1290  if(tr.regionlist != NULL) free(tr.regionlist );
1291  if(tr.edgelist != NULL) free(tr.edgelist );
1292  if(tr.edgemarkerlist != NULL) free(tr.edgemarkerlist );
1293  if(tr.normlist != NULL) free(tr.normlist );
1294  };
1298 
1299 
1303 
1304 
1308 
1309 
1313 
1316 
1317 
1319 
1320 }; // Class TriGenGLASSModeler
1321 
1323 
1326 
1327 
1331 
1332 
1334 inline std::istream& operator >> (std::istream& rIStream,
1335  TriGenGLASSModeler& rThis);
1336 
1338 inline std::ostream& operator << (std::ostream& rOStream,
1339  const TriGenGLASSModeler& rThis)
1340 {
1341  rThis.PrintInfo(rOStream);
1342  rOStream << std::endl;
1343  rThis.PrintData(rOStream);
1344 
1345  return rOStream;
1346 }
1348 
1349 
1350 } // namespace Kratos.
1351 
1352 #endif // KRATOS_TRIGEN_GLASS_MODELER_H_INCLUDED defined
PeriodicInterfaceProcess & operator=(const PeriodicInterfaceProcess &)=delete
Short class definition.
Definition: bucket.h:57
Base class for all Conditions.
Definition: condition.h:59
Geometry< NodeType >::PointsArrayType NodesArrayType
definition of nodes container type, redefined from GeometryType
Definition: condition.h:86
Dof represents a degree of freedom (DoF).
Definition: dof.h:86
Base class for all Elements.
Definition: element.h:60
virtual Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const
It creates a new element pointer.
Definition: element.h:202
This process removes the entities from a model part with the flag TO_ERASE.
Definition: entity_erase_process.h:70
void Execute() override
Execute method is used to execute the Process algorithms.
Geometry base class.
Definition: geometry.h:71
This class is a wrapper for a pointer to a data that is located in a different rank.
Definition: global_pointer.h:44
Definition: amatrix_interface.h:497
Definition: amatrix_interface.h:530
PropertiesType::Pointer pGetProperties(IndexType PropertiesId)
Definition: mesh.h:394
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
NodeType::Pointer CreateNewNode(int Id, double x, double y, double z, VariablesList::Pointer pNewVariablesList, IndexType ThisIndex=0)
Definition: model_part.cpp:270
ConditionIterator ConditionsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1361
SizeType GetNodalSolutionStepDataSize()
Definition: model_part.h:571
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
ConditionsContainerType & Conditions(IndexType ThisIndex=0)
Definition: model_part.h:1381
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
ElementIterator ElementsEnd(IndexType ThisIndex=0)
Definition: model_part.h:1179
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
MeshType & GetMesh(IndexType ThisIndex=0)
Definition: model_part.h:1791
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
ConditionIterator ConditionsEnd(IndexType ThisIndex=0)
Definition: model_part.h:1371
This class defines the node.
Definition: node.h:65
std::vector< std::unique_ptr< Dof< double > >> DofsContainerType
The DoF container type definition.
Definition: node.h:92
double Y() const
Definition: point.h:187
double Z() const
Definition: point.h:193
double X() const
Definition: point.h:181
void reserve(size_type dim)
Definition: pointer_vector.h:319
static double ElapsedSeconds(const std::chrono::steady_clock::time_point StartTime)
This method returns the resulting time.
Definition: timer.h:179
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
Short class definition.
Definition: trigen_glass_forming.h:61
std::vector< PointType::Pointer > PointVector
Definition: trigen_glass_forming.h:70
PointVector::iterator PointIterator
Definition: trigen_glass_forming.h:71
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: trigen_glass_forming.h:645
Tree< KDTreePartition< BucketType > > KdtreeType
Definition: trigen_glass_forming.h:75
void ReGenerateMesh(ModelPart &ThisModelPart, Element const &rReferenceElement, Condition const &rReferenceBoundaryCondition, EntitiesEraseProcess< Node > &node_erase, bool rem_nodes=true, bool add_nodes=true, double my_alpha=1.4, double h_factor=0.5)
Definition: trigen_glass_forming.h:105
TriGenGLASSModeler()
Default constructor.
Definition: trigen_glass_forming.h:82
virtual std::string Info() const
Turn back information as a string.
Definition: trigen_glass_forming.h:636
std::vector< double > DistanceVector
Definition: trigen_glass_forming.h:72
std::vector< double >::iterator DistanceIterator
Definition: trigen_glass_forming.h:73
KRATOS_CLASS_POINTER_DEFINITION(TriGenGLASSModeler)
Pointer definition of TriGenModeler.
virtual void PrintInfo(std::ostream &rOStream) const
Print information about this object.
Definition: trigen_glass_forming.h:642
Bucket< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType
Definition: trigen_glass_forming.h:74
Node PointType
Definition: trigen_glass_forming.h:68
Node::Pointer PointPointerType
Definition: trigen_glass_forming.h:69
virtual ~TriGenGLASSModeler()
Destructor.
Definition: trigen_glass_forming.h:90
A three node 2D triangle geometry with linear shape functions.
Definition: triangle_2d_3.h:74
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_WATCH(variable)
Definition: define.h:806
#define KRATOS_TRY
Definition: define.h:109
#define REAL
Definition: delaunator_utilities.cpp:16
Kratos::ModelPart ModelPart
Definition: kratos_wrapper.h:31
z
Definition: GenerateWind.py:163
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
bool CalculatePosition(Geometry< Node > &geom, const double xc, const double yc, const double zc, array_1d< double, 3 > &N)
Definition: transfer_utility.h:343
double CalculateVol(const double x0, const double y0, const double x1, const double y1, const double x2, const double y2)
Definition: transfer_utility.h:424
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
void triangulate(char *, struct triangulateio *, struct triangulateio *, struct triangulateio *)
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
PointVector::iterator PointIterator
Definition: internal_variables_interpolation_process.h:54
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
std::vector< double > DistanceVector
Definition: internal_variables_interpolation_process.h:55
int step
Definition: face_heat.py:88
y
Other simbols definition.
Definition: generate_axisymmetric_navier_stokes_element.py:54
res
Definition: generate_convection_diffusion_explicit_element.py:211
h
Definition: generate_droplet_dynamics.py:91
x2
Definition: generate_frictional_mortar_condition.py:122
x1
Definition: generate_frictional_mortar_condition.py:121
int tol
Definition: hinsberg_optimization.py:138
R
Definition: isotropic_damage_automatic_differentiation.py:172
float radius
Definition: mesh_to_mdpa_converter.py:18
def Interpolate(variable, entity, sf_values, historical_value)
Definition: point_output_process.py:231
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
el
Definition: read_stl.py:25
float temp
Definition: rotating_cone.py:85
float xc
Definition: rotating_cone.py:77
float yc
Definition: rotating_cone.py:78
int counter
Definition: script_THERMAL_CORRECT.py:218
add_nodes
Definition: script.py:96
N
Definition: sensitivityMatrix.py:29
x
Definition: sensitivityMatrix.py:49
namespace
Definition: array_1d.h:793
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31
float factor
for node in (self.combined_model_part).Nodes: pold = node.GetSolutionStepValue(PRESSURE,...
Definition: ulf_PGLASS.py:254