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_pfem_refine.h
Go to the documentation of this file.
1 //
2 // Project Name: Kratos
3 // Last Modified by: $Author: rrossi $
4 // Date: $Date: 2009-01-22 17:13:57 $
5 // Revision: $Revision: 1.5 $
6 //
7 //
8 
9 
10 #if !defined(KRATOS_TRIGEN_PFEM_MODELER_H_INCLUDED )
11 #define KRATOS_TRIGEN_PFEM_MODELER_H_INCLUDED
12 
13 // System includes
14 #include <string>
15 #include <iostream>
16 #include <cstdlib>
17 
18 // External includes
19 #include "triangle.h"
20 
21 // Project includes
22 #include "includes/define.h"
23 #include "includes/model_part.h"
24 #include "utilities/timer.h"
30 
31 
32 
33 namespace Kratos
34 {
35 extern "C" {
36  void triangulate(char *, struct triangulateio *, struct triangulateio *,struct triangulateio *);
37  //void trifree();
38 }
39 
40 
43 
47 
51 
55 
59 
61 
64 {
65 public:
68 
71  typedef Node PointType;
72  typedef Node::Pointer PointPointerType;
73  typedef std::vector<PointType::Pointer> PointVector;
74  typedef PointVector::iterator PointIterator;
75  typedef std::vector<double> DistanceVector;
76  typedef std::vector<double>::iterator DistanceIterator;
79 
83 
86  mJ(ZeroMatrix(2,2)), //local jacobian
87  mJinv(ZeroMatrix(2,2)), //inverse jacobian
88  mC(ZeroVector(2)), //dimension = number of nodes
89  mRhs(ZeroVector(2)) {}
90  //mpNodeEraseProcess(NULL){} //dimension = number of nodes
91 
93  virtual ~TriGenPFEMModeler() = default;
94 
95 
99 
100 
104 
105 
106  //*******************************************************************************************
107  //*******************************************************************************************
109  ModelPart& ThisModelPart ,
110  Element const& rReferenceElement,
111  Condition const& rReferenceBoundaryCondition,
112  EntitiesEraseProcess<Node>& node_erase, bool rem_nodes = true, bool add_nodes=true,
113  double my_alpha = 1.4, double h_factor=0.5)
114  {
115 
116  KRATOS_TRY
117  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FREE_SURFACE)==false )
118  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_FREE_SURFACE---- variable!!!!!! ERROR","");
119  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_STRUCTURE)==false )
120  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_STRUCTURE---- variable!!!!!! ERROR","");
121  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_BOUNDARY)==false )
122  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_BOUNDARY---- variable!!!!!! ERROR","");
123  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FLUID)==false )
124  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_FLUID---- variable!!!!!! ERROR","");
125 
126  KRATOS_WATCH("Trigen PFEM Refining Mesher")
127  const auto inital_time = std::chrono::steady_clock::now();
128 
129 
130 //clearing elements
131 
132  ThisModelPart.Elements().clear();
133  ThisModelPart.Conditions().clear();
134 
136 
137 
138  int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
139 
140  // bucket types
141  //typedef Bucket<3, PointType, ModelPart::NodesContainerType, PointPointerType, PointIterator, DistanceIterator > BucketType;
142  //typedef Bins< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > StaticBins;
143  // bucket types
144 
145  //typedef Tree< StaticBins > Bin; //Binstree;
146  unsigned int bucket_size = 20;
147 
148  //performing the interpolation - all of the nodes in this list will be preserved
149  unsigned int max_results = 100;
150  //PointerVector<PointType> res(max_results);
151  //NodeIterator res(max_results);
152  PointVector res(max_results);
153  DistanceVector res_distances(max_results);
154  Node work_point(0,0.0,0.0,0.0);
155  //if the remove_node switch is activated, we check if the nodes got too close
156 
157  if (rem_nodes==true)
158  {
159  PointVector list_of_nodes;
160  list_of_nodes.reserve(ThisModelPart.Nodes().size());
161  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
162  {
163  (list_of_nodes).push_back(*(i_node.base()));
164  }
165 
166  KdtreeType nodes_tree1(list_of_nodes.begin(),list_of_nodes.end(), bucket_size);
167 
168  RemoveCloseNodes(ThisModelPart, nodes_tree1, node_erase, h_factor);
169  }
170 
172  // //
173  // Now we shall pass the Alpha Shape for the second time, having the "bad nodes" removed //
175  //creating the containers for the input and output
176  struct triangulateio in_mid;
177  struct triangulateio out_mid;
178  struct triangulateio vorout_mid;
179 
180  initialize_triangulateio(in_mid);
181  initialize_triangulateio(out_mid);
182  initialize_triangulateio(vorout_mid);
183 
184  //assigning the size of the input containers
185 
186  in_mid.numberofpoints = ThisModelPart.Nodes().size();
187  in_mid.pointlist = (REAL *) malloc(in_mid.numberofpoints * 2 * sizeof(REAL));
188 
189  //reorder the id's and give the coordinates to the mesher
190  ModelPart::NodesContainerType::iterator nodes_begin = ThisModelPart.NodesBegin();
191  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
192  {
193  int base = i*2;
194  //int base = ((nodes_begin + i)->Id() - 1 ) * 2;
195 
196  //from now on it is consecutive
197  (nodes_begin + i)->SetId(i+1);
198 // (nodes_begin + i)->Id() = i+1;
199 
200  in_mid.pointlist[base] = (nodes_begin + i)->X();
201  in_mid.pointlist[base+1] = (nodes_begin + i)->Y();
202 
203  auto& node_dofs = (nodes_begin + i)->GetDofs();
204  for(auto iii = node_dofs.begin(); iii != node_dofs.end(); iii++)
205  {
206  (**iii).SetEquationId(i+1);
207  }
208  }
209  //in_mid.numberoftriangles = ThisModelPart.Elements().size();
210  //in_mid.trianglelist = (int*) malloc(in_mid.numberoftriangles * 3 * sizeof(int));
211 
212  // "P" suppresses the output .poly file. Saves disk space, but you
213  //lose the ability to maintain constraining segments on later refinements of the mesh.
214  // "B" Suppresses boundary markers in the output .node, .poly, and .edge output files
215  // "n" outputs a list of triangles neighboring each triangle.
216  // "e" outputs edge list (i.e. all the "connectivities")
217  char options1[] = "Pne";
218  triangulate(options1, &in_mid, &out_mid, &vorout_mid);
219  //print out the mesh generation time
220  std::cout << "Mesh generation time = " << Timer::ElapsedSeconds(inital_time) << std::endl;
221  //number of newly generated triangles
222  unsigned int el_number=out_mid.numberoftriangles;
223 
224  //PASSING THE CREATED ELEMENTS TO THE ALPHA-SHAPE
225  std::vector<int> preserved_list1(el_number);
226  preserved_list1.resize(el_number);
227 
228  int number_of_preserved_elems= PassAlphaShape(ThisModelPart, preserved_list1, el_number, my_alpha, out_mid);
229 
230  //freeing memory
231 
232  clean_triangulateio(in_mid);
233  clean_triangulateio(vorout_mid);
234  KRATOS_WATCH("ln367");
235  //NOW WE SHALL PERFORM ADAPTIVE REMESHING, i.e. insert and remove nodes based upon mesh quality
236  // and prescribed sizes
237  struct triangulateio in2;
238  struct triangulateio out2;
239  struct triangulateio vorout2;
240 
241  initialize_triangulateio(in2);
242  initialize_triangulateio(out2);
243  initialize_triangulateio(vorout2);
244 
245 // in2.firstnumber = 1;
246  in2.numberofpoints = ThisModelPart.Nodes().size();
247  in2.pointlist = (REAL *) malloc(in2.numberofpoints * 2 * sizeof(REAL));
248 
249  //writing the input point list
250  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
251  {
252  int base = i*2;
253  in2.pointlist[base] = (nodes_begin + i)->X();
254  in2.pointlist[base+1] = (nodes_begin + i)->Y();
255  }
256  in2.numberoftriangles=number_of_preserved_elems;
257 
258  in2.trianglelist = (int *) malloc(in2.numberoftriangles * 3 * sizeof(int));
259  in2.trianglearealist = (REAL *) malloc(in2.numberoftriangles * sizeof(REAL));
260 // in2.trianglelist = new int[in2.numberoftriangles * 3];
261 // in2.trianglearealist = new double[in2.numberoftriangles];
262 
263  KRATOS_WATCH(el_number);
264  int counter = 0;
265  //here I will assign a huge number of NODAL_H to the free surface nodes, so that there no nodes will be added
266  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
267  {
268  if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
269  {
270 
271  double& val=i_node->FastGetSolutionStepValue(NODAL_H);
272  val*=2.0;
273  //i_node->FastGetSolutionStepValue(NODAL_H,1)=val;
274  //KRATOS_WATCH("AAAAAAAAAAAAAAAAAAAAAAAA!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
275  }
276  }
277 
278 
279  for (unsigned int el=0; el<el_number; el++)
280  {
281  if( static_cast<bool>(preserved_list1[el]) == true )
282  {
283  //saving the list of ONLY preserved triangles, the ones that passed alpha-shape check
284  int new_base = counter*3;
285  int old_base = el*3;
286  //copying in case it is preserved
287  in2.trianglelist[new_base] = out_mid.trianglelist[old_base];
288  in2.trianglelist[new_base+1] = out_mid.trianglelist[old_base+1];
289  in2.trianglelist[new_base+2] = out_mid.trianglelist[old_base+2];
290 
291  //calculate the prescribed h
292  double prescribed_h = (nodes_begin + out_mid.trianglelist[old_base]-1)->FastGetSolutionStepValue(NODAL_H);
293  prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+1]-1)->FastGetSolutionStepValue(NODAL_H);
294  prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+2]-1)->FastGetSolutionStepValue(NODAL_H);
295  prescribed_h *= 0.3333;
296  //if h is the height of a equilateral triangle, the area is 1/2*h*h
297  in2.trianglearealist[counter] = 0.5*(1.5*prescribed_h*1.5*prescribed_h);
298  counter += 1;
299  }
300 
301  }
302  //now I set back the nodal_h
303  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
304  {
305  if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
306  {
307  double& nodal_h=i_node->FastGetSolutionStepValue(NODAL_H);
308  nodal_h/=2.0;
309  }
310  }
311 
312 
313  clean_triangulateio(out_mid);
314  KRATOS_WATCH("ln420");
315  //here we generate a new mesh adding/removing nodes, by initializing "q"-quality mesh and "a"-area constraint switches
316  //
317  // MOST IMPORTANT IS "r" switch, that refines previously generated mesh!!!!!!!!!!(that is the one given inside in2)
318  //char mesh_regen_opts[] = "YYJaqrn";
319  //char mesh_regen_opts[] = "YJq1.4arn";
320  if (add_nodes==true)
321  {
322  char mesh_regen_opts[] = "YJq1.4arn";
323  triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
324  KRATOS_WATCH("Adaptive remeshing executed")
325  }
326  else
327  {
328  char mesh_regen_opts[] = "YJrn";
329  triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
330  KRATOS_WATCH("Non-Adaptive remeshing executed")
331  }
332 
333  //and now we shall find out where the new nodes belong to
334  //defintions for spatial search
335  typedef Node PointType;
336  typedef Node::Pointer PointPointerType;
337  typedef std::vector<PointType::Pointer> PointVector;
338  //typedef std::vector<PointType::Pointer>::iterator PointIterator;
339  typedef PointVector::iterator PointIterator;
340  typedef std::vector<double> DistanceVector;
341  typedef std::vector<double>::iterator DistanceIterator;
342  KRATOS_WATCH("ln449");
343 
345 
346  typedef Tree< KDTreePartition<BucketType> > kd_tree; //Kdtree;
347 
348  //int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
349 
350  //creating an auxiliary list for the new nodes
351  PointVector list_of_new_nodes;
352 
353  //node to get the DOFs from
354  Node::DofsContainerType& reference_dofs = (ThisModelPart.NodesBegin())->GetDofs();
355 
356  double z = 0.0;
357  int n_points_before_refinement = in2.numberofpoints;
358  //if points were added, we add them as nodes to the ModelPart
359  if (out2.numberofpoints > n_points_before_refinement )
360  {
361  for(int i = n_points_before_refinement; i<out2.numberofpoints; i++)
362  {
363  int id=i+1;
364  int base = i*2;
365  double& x= out2.pointlist[base];
366  double& y= out2.pointlist[base+1];
367 
368  Node::Pointer pnode = ThisModelPart.CreateNewNode(id,x,y,z);
369 
370  pnode->SetBufferSize(ThisModelPart.NodesBegin()->GetBufferSize() );
371 
372  list_of_new_nodes.push_back( pnode );
373 
374  //std::cout << "new node id = " << pnode->Id() << std::endl;
375  //generating the dofs
376  for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
377  {
378  Node::DofType& rDof = **iii;
379  Node::DofType::Pointer p_new_dof = pnode->pAddDof( rDof );
380 
381  (p_new_dof)->FreeDof();
382 // (p_new_dof)->EquationId() = -1;
383 
384  }
385 
386  }
387  }
388 
389  std::cout << "During refinement we added " << out2.numberofpoints-n_points_before_refinement<< " nodes " <<std::endl;
390  //unsigned int bucket_size = 20;
391  //performing the interpolation - all of the nodes in this list will be preserved
392  //unsigned int max_results = 50;
393  //PointVector results(max_results);
394  //DistanceVector results_distances(max_results);
396 
398 
399  //int number_of_preserved_elems=0;
400 
401  int point_base;
402  //WHAT ARE THOSE????
403 // Node work_point(0,0.0,0.0,0.0);
404  unsigned int MaximumNumberOfResults = list_of_new_nodes.size();
405  PointVector Results(MaximumNumberOfResults);
406  DistanceVector ResultsDistances(MaximumNumberOfResults);
407 
408  step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
409 
410  if(out2.numberofpoints-n_points_before_refinement > 0) //if we added points
411  {
412 
413  kd_tree nodes_tree2(list_of_new_nodes.begin(),list_of_new_nodes.end(),bucket_size);
414  nodes_begin = ThisModelPart.NodesBegin();
415 
416  for(int el = 0; el< in2.numberoftriangles; el++)
417  {
418  int base = el * 3;
419  //coordinates
420  point_base = (in2.trianglelist[base] - 1)*2;
421  x1[0] = in2.pointlist[point_base];
422  x1[1] = in2.pointlist[point_base+1];
423 
424  point_base = (in2.trianglelist[base+1] - 1)*2;
425  x2[0] = in2.pointlist[point_base];
426  x2[1] = in2.pointlist[point_base+1];
427 
428  point_base = (in2.trianglelist[base+2] - 1)*2;
429  x3[0] = in2.pointlist[point_base];
430  x3[1] = in2.pointlist[point_base+1];
431 
432 
433  //find the center and "radius" of the element
434  double xc, yc, radius;
435  CalculateCenterAndSearchRadius( x1[0], x1[1],
436  x2[0], x2[1],
437  x3[0], x3[1],
438  xc,yc,radius);
439 
440  //find all of the new nodes within the radius
441  int number_of_points_in_radius;
442  work_point.X() = xc;
443  work_point.Y() = yc;
444  work_point.Z() = 0.0;
445 
446  number_of_points_in_radius = nodes_tree2.SearchInRadius(work_point, radius*1.01, Results.begin(),
447  ResultsDistances.begin(), MaximumNumberOfResults);
448 
449  Triangle2D3<Node > geom(
450  *( (nodes_begin + in2.trianglelist[base]-1).base() ),
451  *( (nodes_begin + in2.trianglelist[base+1]-1).base() ),
452  *( (nodes_begin + in2.trianglelist[base+2]-1).base() )
453  );
454 
455  //check if inside and eventually interpolate
456  for( auto it_found = Results.begin(); it_found != Results.begin() + number_of_points_in_radius; it_found++)
457  {
458  bool is_inside = false;
459  is_inside = CalculatePosition(x1[0], x1[1],
460  x2[0], x2[1],
461  x3[0], x3[1],
462  (*it_found)->X(),(*it_found)->Y(),N);
463 
464 
465  if(is_inside == true)
466  {
467  Interpolate( geom, N, step_data_size, *(it_found ) );
468 
469  }
470  }
471  }
472  }
473 
474  ThisModelPart.Elements().clear();
475 
476  //set the coordinates to the original value
477  for(auto & list_of_new_node : list_of_new_nodes)
478  {
479  const array_1d<double,3>& disp = list_of_new_node->FastGetSolutionStepValue(DISPLACEMENT);
480  list_of_new_node->X0() = list_of_new_node->X() - disp[0];
481  list_of_new_node->Y0() = list_of_new_node->Y() - disp[1];
482  list_of_new_node->Z0() = 0.0;
483  }
484  //cleaning unnecessary data
485  //in2.deinitialize();
486  //in2.initialize();
487  //free( in2.pointlist );
488 
489  //add preserved elements to the kratos
490  Properties::Pointer properties = ThisModelPart.GetMesh().pGetProperties(1);
491  nodes_begin = ThisModelPart.NodesBegin();
492  (ThisModelPart.Elements()).reserve(out2.numberoftriangles);
493 
494  for(int iii = 0; iii< out2.numberoftriangles; iii++)
495  {
496  int id = iii + 1;
497  int base = iii * 3;
498  Triangle2D3<Node > geom(
499  *( (nodes_begin + out2.trianglelist[base]-1).base() ),
500  *( (nodes_begin + out2.trianglelist[base+1]-1).base() ),
501  *( (nodes_begin + out2.trianglelist[base+2]-1).base() )
502  );
503 
504 
505 #ifdef _DEBUG
506  ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
507  if( *(ModelNodes).find( out2.trianglelist[base]).base() == *(ThisModelPart.Nodes().end()).base() )
508  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
509  if( *(ModelNodes).find( out2.trianglelist[base+1]).base() == *(ThisModelPart.Nodes().end()).base() )
510  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
511  if( *(ModelNodes).find( out2.trianglelist[base+2]).base() == *(ThisModelPart.Nodes().end()).base() )
512  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
513 
514 #endif
515 
516  Element::Pointer p_element = rReferenceElement.Create(id, geom, properties);
517  (ThisModelPart.Elements()).push_back(p_element);
518 
519  }
520  ThisModelPart.Elements().Sort();
521 
522  //filling the neighbour list
523  ModelPart::ElementsContainerType::const_iterator el_begin = ThisModelPart.ElementsBegin();
524  for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.ElementsBegin();
525  iii != ThisModelPart.ElementsEnd(); iii++)
526  {
527  //Geometry< Node >& geom = iii->GetGeometry();
528  int base = ( iii->Id() - 1 )*3;
529 
530  (iii->GetValue(NEIGHBOUR_ELEMENTS)).resize(3);
531  GlobalPointersVector< Element >& neighb = iii->GetValue(NEIGHBOUR_ELEMENTS);
532  for(int i = 0; i<3; i++)
533  {
534  int index = out2.neighborlist[base+i];
535  if(index > 0)
536  neighb(i) = GlobalPointer<Element>(&*(el_begin + index-1));
537  else
538  neighb(i) = Element::WeakPointer();
539  }
540  }
541  //identifying boundary, creating skin
542  IdentifyBoundary(ThisModelPart, rReferenceBoundaryCondition, properties, out2);
543 
544  KRATOS_WATCH("ln749");
545 
546  clean_triangulateio(in2);
547  KRATOS_WATCH("ln752");
548  clean_triangulateio(out2);
549  KRATOS_WATCH("ln754");
550  clean_triangulateio(vorout2);
551  KRATOS_WATCH("Finished remeshing with Trigen_PFEM_Refine")
552 
553  KRATOS_CATCH("")
554  }
555 
556 
560 
561 
565 
566 
570 
572  virtual std::string Info() const
573  {
574  return "";
575  }
576 
578  virtual void PrintInfo(std::ostream& rOStream) const {}
579 
581  virtual void PrintData(std::ostream& rOStream) const {}
582 
583 
587 
588 
590 
591 protected:
594 
595 
599  //KdtreeType mKdtree;
600 
604 
605 
609 
610 
614 
615 
619 
620 
624 
625 
627 
628 private:
631 
632 
636  boost::numeric::ublas::bounded_matrix<double,2,2> mJ; //local jacobian
637  boost::numeric::ublas::bounded_matrix<double,2,2> mJinv; //inverse jacobian
638  array_1d<double,2> mC; //center pos
639  array_1d<double,2> mRhs; //center pos
640  //EntitiesEraseProcess<Node>* mpNodeEraseProcess;
641 
642 
646  void RemoveCloseNodes(ModelPart& ThisModelPart, KdtreeType& nodes_tree1, EntitiesEraseProcess<Node>& node_erase, double& h_factor)
647  {
648  //unsigned int bucket_size = 20;
649 
650  //performing the interpolation - all of the nodes in this list will be preserved
651  unsigned int max_results = 100;
652  //PointerVector<PointType> res(max_results);
653  //NodeIterator res(max_results);
654  PointVector res(max_results);
655  DistanceVector res_distances(max_results);
656  Node work_point(0,0.0,0.0,0.0);
657 
658  unsigned int n_points_in_radius;
659  //radius means the distance, closer than which no node shall be allowd. if closer -> mark for erasing
660  double radius;
661 
662  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in != ThisModelPart.NodesEnd(); in++)
663  {
664  radius=h_factor*in->FastGetSolutionStepValue(NODAL_H);
665 
666  work_point[0]=in->X();
667  work_point[1]=in->Y();
668  work_point[2]=in->Z();
669 
670  n_points_in_radius = nodes_tree1.SearchInRadius(work_point, radius, res.begin(),res_distances.begin(), max_results);
671  if (n_points_in_radius>1)
672  {
673  if (in->FastGetSolutionStepValue(IS_BOUNDARY)==0.0 && in->FastGetSolutionStepValue(IS_STRUCTURE)==0.0)
674  {
675  //look if we are already erasing any of the other nodes
676  double erased_nodes = 0;
677  for(auto i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
678  erased_nodes += in->Is(TO_ERASE);
679 
680  if( erased_nodes < 1) //we cancel the node if no other nodes are being erased
681  in->Set(TO_ERASE,true);
682 
683  }
684  else if ( (in)->FastGetSolutionStepValue(IS_STRUCTURE)!=1.0) //boundary nodes will be removed if they get REALLY close to another boundary node (0.2 * h_factor)
685  {
686  //here we loop over the neighbouring nodes and if there are nodes
687  //with IS_BOUNDARY=1 which are closer than 0.2*nodal_h from our we remove the node we are considering
688  unsigned int k = 0;
689  unsigned int counter = 0;
690  for(auto i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
691  {
692  if ( (*i)->FastGetSolutionStepValue(IS_BOUNDARY,1)==1.0 && res_distances[k] < 0.2*radius && res_distances[k] > 0.0 )
693  {
694  // KRATOS_WATCH( res_distances[k] );
695  counter += 1;
696  }
697  k++;
698  }
699  if(counter > 0)
700  in->Set(TO_ERASE,true);
701  }
702  }
703 
704  }
705 
706 
707  node_erase.Execute();
708 
709  KRATOS_WATCH("Number of nodes after erasing")
710  KRATOS_WATCH(ThisModelPart.Nodes().size())
711  }
712 
713 
714  int PassAlphaShape(ModelPart& ThisModelPart, std::vector<int>& preserved_list1, unsigned int & el_number, double& my_alpha, struct triangulateio& out_mid)
715  {
716  //NOTE THAT preserved_list1 will be overwritten, only the elements that passed alpha-shaoe check will enter it
717 
718  //prepairing for alpha shape passing : creating necessary arrays
719  //list of preserved elements is created: at max el_number can be preserved (all elements)
720 
721 
722  array_1d<double,3> x1,x2,x3,xc;
723 
724  //int number_of_preserved_elems=0;
725  int number_of_preserved_elems=0;
726  int point_base;
727  //loop for passing alpha shape
728  for(unsigned int el = 0; el< el_number; el++)
729  {
730  int base = el * 3;
731 
732  //coordinates
733  point_base = (out_mid.trianglelist[base] - 1)*2;
734  x1[0] = out_mid.pointlist[point_base];
735  x1[1] = out_mid.pointlist[point_base+1];
736 
737  point_base = (out_mid.trianglelist[base+1] - 1)*2;
738  x2[0] = out_mid.pointlist[point_base];
739  x2[1] = out_mid.pointlist[point_base+1];
740 
741  point_base = (out_mid.trianglelist[base+2] - 1)*2;
742  x3[0] = out_mid.pointlist[point_base];
743  x3[1] = out_mid.pointlist[point_base+1];
744 
745  //here we shall temporarily save the elements and pass them afterwards to the alpha shape
746  Geometry<Node > temp;
747 
748  temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base]).base() ) );
749  temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+1]).base() ) );
750  temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+2]).base() ) );
751 
752  int number_of_structure_nodes = int( temp[0].FastGetSolutionStepValue(IS_STRUCTURE) );
753  number_of_structure_nodes += int( temp[1].FastGetSolutionStepValue(IS_STRUCTURE) );
754  number_of_structure_nodes += int( temp[2].FastGetSolutionStepValue(IS_STRUCTURE) );
755 
756  //check the number of nodes of boundary
757  int nfs = int( temp[0].FastGetSolutionStepValue(IS_FREE_SURFACE) );
758  nfs += int( temp[1].FastGetSolutionStepValue(IS_FREE_SURFACE) );
759  nfs += int( temp[2].FastGetSolutionStepValue(IS_FREE_SURFACE) );
760 
761  //check the number of nodes of boundary
762  int nfluid = int( temp[0].FastGetSolutionStepValue(IS_FLUID) );
763  nfluid += int( temp[1].FastGetSolutionStepValue(IS_FLUID) );
764  nfluid += int( temp[2].FastGetSolutionStepValue(IS_FLUID) );
765 
766  //check the number of nodes of boundary
767  // int nboundary = int( temp[0].FastGetSolutionStepValue(IS_BOUNDARY) );
768  // nboundary += int( temp[1].FastGetSolutionStepValue(IS_BOUNDARY) );
769  // nboundary += int( temp[2].FastGetSolutionStepValue(IS_BOUNDARY) );
770  //first check that we are working with fluid elements, otherwise throw an error
771  //if (nfluid!=3)
772  // KRATOS_THROW_ERROR(std::logic_error,"THATS NOT FLUID or NOT TRIANGLE!!!!!! ERROR","");
773  //otherwise perform alpha shape check
774 
775 
776  if(number_of_structure_nodes!=3) //if it is = 3 it is a completely fixed element -> do not add it
777  {
778 
779  if (nfs != 0 || nfluid != 3) //in this case it is close to the surface so i should use alpha shape
780  {
781 
782  if( AlphaShape(my_alpha,temp) && number_of_structure_nodes!=3) //if alpha shape says to preserve
783  {
784 // if(nboundary==3 && number_of_structure_nodes > 1 && nfs > 0) //if it is = 3 pressure problems -> do not add it
785 // {
786 // preserved_list1[el] = false;
787 // }
788 // else
789 // {
790  preserved_list1[el] = true;
791  number_of_preserved_elems += 1;
792 // }
793  }
794  }
795  else //internal triangle --- should be ALWAYS preserved
796  {
797  double bigger_alpha = my_alpha*10.0;
798  if( AlphaShape(bigger_alpha,temp) && number_of_structure_nodes!=3)
799  {
800 // if(nboundary==3 && number_of_structure_nodes > 1 && nfs > 0) //if it is = 3 pressure problems -> do not add it
801 // {
802 // preserved_list1[el] = false;
803 // }
804 // else
805 // {
806  preserved_list1[el] = true;
807  number_of_preserved_elems += 1;
808 // }
809  }
810  }
811  }
812  else
813  preserved_list1[el] = false;
814 
815  }
816  return number_of_preserved_elems;
817  }
818 
819 
820 
821 
822 
823  void IdentifyBoundary(ModelPart& ThisModelPart, Condition const& rReferenceBoundaryCondition, Properties::Pointer& properties, struct triangulateio& out2)
824  {
825 
826  //reset the boundary flag
827  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in!=ThisModelPart.NodesEnd(); in++)
828  {
829  in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
830  }
831  //filling the elemental neighbours list (from now on the elements list can not change)
832  ModelPart::ElementsContainerType::iterator elements_end = ThisModelPart.Elements().end();
833 
834  ThisModelPart.Elements().Unique();
835 
836  //now the boundary faces
837  for(ModelPart::ElementsContainerType::iterator iii = ThisModelPart.ElementsBegin(); iii != ThisModelPart.ElementsEnd(); iii++)
838  {
839  int base = ( iii->Id() - 1 )*3;
840 
841  ModelPart::ElementsContainerType::iterator el_neighb;
842  /*each face is opposite to the corresponding node number so
843  0 ----- 1 2
844  1 ----- 2 0
845  2 ----- 0 1
846  */
847 
849  //
850  //********************************************************************
851  //first face
852  el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base] );
853  if( el_neighb == elements_end )
854  {
855  //std::cout << "node0" << std::endl;
856  //if no neighnour is present => the face is free surface
857  iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
858  iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
859 
860  //Generate condition
862  temp1.reserve(2);
863  temp1.push_back(iii->GetGeometry()(1));
864  temp1.push_back(iii->GetGeometry()(2));
865 
866  Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
867  int id = (iii->Id()-1)*3;
868  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
869  ThisModelPart.Conditions().push_back(p_cond);
870 
871  }
872 
873  //********************************************************************
874  //second face
875  el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+1] );
876  //if( el != ThisModelPart.Elements().end() )
877  if( el_neighb == elements_end )
878  {
879  //if no neighnour is present => the face is free surface
880  iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
881  iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
882 
883  //Generate condition
885  temp1.reserve(2);
886  temp1.push_back(iii->GetGeometry()(2));
887  temp1.push_back(iii->GetGeometry()(0));
888 
889  Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
890  int id = (iii->Id()-1)*3+1;
891  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
892  ThisModelPart.Conditions().push_back(p_cond);
893 
894 
895  }
896 
897  //********************************************************************
898  //third face
899  el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+2] );
900  if( el_neighb == elements_end )
901  {
902  //if no neighnour is present => the face is free surface
903  iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
904  iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
905 
906 // Generate condition
908  temp1.reserve(2);
909  temp1.push_back(iii->GetGeometry()(0));
910  temp1.push_back(iii->GetGeometry()(1));
911 
912  Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
913  int id = (iii->Id()-1)*3+2;
914  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
915  ThisModelPart.Conditions().push_back(p_cond);
916 
917  }
918 
919 
920  }
921  }
922 
923 
924  //returns false if it should be removed
925  bool AlphaShape(double alpha_param, Geometry<Node >& pgeom)
926  //bool AlphaShape(double alpha_param, Triangle2D<Node >& pgeom)
927  {
928  KRATOS_TRY
929 
930 
931  double x0 = pgeom[0].X();
932  double x1 = pgeom[1].X();
933  double x2 = pgeom[2].X();
934 
935  double y0 = pgeom[0].Y();
936  double y1 = pgeom[1].Y();
937  double y2 = pgeom[2].Y();
938 
939  mJ(0,0)=2.0*(x1-x0);
940  mJ(0,1)=2.0*(y1-y0);
941  mJ(1,0)=2.0*(x2-x0);
942  mJ(1,1)=2.0*(y2-y0);
943 
944 
945  double detJ = mJ(0,0)*mJ(1,1)-mJ(0,1)*mJ(1,0);
946 
947  mJinv(0,0) = mJ(1,1);
948  mJinv(0,1) = -mJ(0,1);
949  mJinv(1,0) = -mJ(1,0);
950  mJinv(1,1) = mJ(0,0);
951 
952  bounded_matrix<double,2,2> check;
953 
954 //calculate average h
955  double h;
956  h = pgeom[0].FastGetSolutionStepValue(NODAL_H);
957  h += pgeom[1].FastGetSolutionStepValue(NODAL_H);
958  h += pgeom[2].FastGetSolutionStepValue(NODAL_H);
959  h *= 0.333333333;
960 
961 
962  if(detJ < 5e-3*h*h)
963  {
964  //std::cout << "detJ = " << detJ << std::endl;
966  pgeom[0].GetSolutionStepValue(IS_BOUNDARY) = 1;
967  pgeom[1].GetSolutionStepValue(IS_BOUNDARY) = 1;
968  pgeom[2].GetSolutionStepValue(IS_BOUNDARY) = 1;
969  return false;
970  }
971 
972  else
973  {
974 
975  double x0_2 = x0*x0 + y0*y0;
976  double x1_2 = x1*x1 + y1*y1;
977  double x2_2 = x2*x2 + y2*y2;
978 
979  //finalizing the calculation of the inverted matrix
980  //std::cout<<"MATR INV"<<MatrixInverse(mJ)<<std::endl;
981  mJinv /= detJ;
982  //calculating the RHS
983  mRhs[0] = (x1_2 - x0_2);
984  mRhs[1] = (x2_2 - x0_2);
985 
986  //calculate position of the center
987  noalias(mC) = prod(mJinv,mRhs);
988 
989  double radius = sqrt(pow(mC[0]-x0,2)+pow(mC[1]-y0,2));
990 
991 
992  if (radius < h*alpha_param)
993  {
994  return true;
995  }
996  else
997  {
998  return false;
999  }
1000  }
1001 
1002 
1003  KRATOS_CATCH("")
1004  }
1005  //AUXILLIARY FCTNS
1006  inline void CalculateCenterAndSearchRadius(const double x0, const double y0,
1007  const double x1, const double y1,
1008  const double x2, const double y2,
1009  double& xc, double& yc, double& R
1010  )
1011  {
1012  xc = 0.3333333333333333333*(x0+x1+x2);
1013  yc = 0.3333333333333333333*(y0+y1+y2);
1014 
1015  double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0);
1016  double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1);
1017  double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2);
1018 
1019  R = R1;
1020  if(R2 > R) R = R2;
1021  if(R3 > R) R = R3;
1022 
1023  R = sqrt(R);
1024  }
1025 
1026  inline double CalculateVol( const double x0, const double y0,
1027  const double x1, const double y1,
1028  const double x2, const double y2
1029  )
1030  {
1031  return 0.5*( (x1-x0)*(y2-y0)- (y1-y0)*(x2-x0) );
1032  }
1033 
1034  inline bool CalculatePosition( const double x0, const double y0,
1035  const double x1, const double y1,
1036  const double x2, const double y2,
1037  const double xc, const double yc,
1038  array_1d<double,3>& N
1039  )
1040  {
1041  double area = CalculateVol(x0,y0,x1,y1,x2,y2);
1042 
1043  if(area < 0.000000000001)
1044  {
1045  KRATOS_THROW_ERROR(std::logic_error,"element with zero area found","");
1046  }
1047 
1048 
1049 
1050  N[0] = CalculateVol(x1,y1,x2,y2,xc,yc) / area;
1051  N[1] = CalculateVol(x2,y2,x0,y0,xc,yc) / area;
1052  N[2] = CalculateVol(x0,y0,x1,y1,xc,yc) / area;
1053 
1054  /* N[0] = CalculateVol(x0,y0,x1,y1,xc,yc) * inv_area;
1055  N[1] = CalculateVol(x1,y1,x2,y2,xc,yc) * inv_area;
1056  N[2] = CalculateVol(x2,y2,x0,y0,xc,yc) * inv_area;*/
1057  double tol = 1e-4;
1058  double upper_limit = 1.0+tol;
1059  double lower_limit = -tol;
1060 
1061  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
1062  //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
1063  return true;
1064 
1065  return false;
1066  }
1067 
1068  void Interpolate( Triangle2D3<Node >& geom, const array_1d<double,3>& N,
1069  unsigned int step_data_size,
1070  Node::Pointer pnode)
1071  {
1072  unsigned int buffer_size = pnode->GetBufferSize();
1073  //KRATOS_WATCH("Buffer size")
1074  //KRATOS_WATCH(buffer_size)
1075 
1076  for(unsigned int step = 0; step<buffer_size; step++)
1077  {
1078 
1079  //getting the data of the solution step
1080  double* step_data = (pnode)->SolutionStepData().Data(step);
1081 
1082 
1083  double* node0_data = geom[0].SolutionStepData().Data(step);
1084  double* node1_data = geom[1].SolutionStepData().Data(step);
1085  double* node2_data = geom[2].SolutionStepData().Data(step);
1086 
1087  //copying this data in the position of the vector we are interested in
1088  for(unsigned int j= 0; j<step_data_size; j++)
1089  {
1090 
1091  step_data[j] = N[0]*node0_data[j] + N[1]*node1_data[j] + N[2]*node2_data[j];
1092 
1093 
1094  }
1095  }
1096  if (N[0]==0.0 && N[1]==0.0 && N[2]==0.0)
1097  KRATOS_THROW_ERROR(std::logic_error,"SOMETHING's wrong with the added nodes!!!!!! ERROR","");
1098 
1099  //if ( pnode->FastGetSolutionStepValue(BULK_MODULUS)==0.0)
1100  // KRATOS_THROW_ERROR(std::logic_error,"SOMETHING's wrong with the added nodes!!!!!! ERROR","");
1101 
1102  //now we assure that the flag variables are set coorect!! since we add nodes inside of the fluid volume only
1103  //we manually reset the IS_BOUNDARY, IS_FLUID, IS_STRUCTURE, IS_FREE_SURFACE values in a right way
1104  //not to have values, like 0.33 0.66 resulting if we would have been interpolating them in the same way
1105  //as the normal variables, like Velocity etc
1106 
1107 
1108  pnode->FastGetSolutionStepValue(IS_BOUNDARY)=0.0;
1109  pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1110  pnode->FastGetSolutionStepValue(IS_INTERFACE)=0.0;
1111  pnode->Set(TO_ERASE,false);
1112  pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=0.0;
1113  pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1114  }
1115 
1116  void initialize_triangulateio( triangulateio& tr )
1117  {
1118  tr.pointlist = (REAL*) nullptr;
1119  tr.pointattributelist = (REAL*) nullptr;
1120  tr.pointmarkerlist = (int*) nullptr;
1121  tr.numberofpoints = 0;
1122  tr.numberofpointattributes = 0;
1123  tr.trianglelist = (int*) nullptr;
1124  tr.triangleattributelist = (REAL*) nullptr;
1125  tr.trianglearealist = (REAL*) nullptr;
1126  tr.neighborlist = (int*) nullptr;
1127  tr.numberoftriangles = 0;
1128  tr.numberofcorners = 3;
1129  tr.numberoftriangleattributes = 0;
1130  tr.segmentlist = (int*) nullptr;
1131  tr.segmentmarkerlist = (int*) nullptr;
1132  tr.numberofsegments = 0;
1133  tr.holelist = (REAL*) nullptr;
1134  tr.numberofholes = 0;
1135  tr.regionlist = (REAL*) nullptr;
1136  tr.numberofregions = 0;
1137  tr.edgelist = (int*) nullptr;
1138  tr.edgemarkerlist = (int*) nullptr;
1139  tr.normlist = (REAL*) nullptr;
1140  tr.numberofedges = 0;
1141  };
1142 
1143  void clean_triangulateio( triangulateio& tr )
1144  {
1145  if(tr.pointlist != nullptr) free(tr.pointlist );
1146  if(tr.pointattributelist != nullptr) free(tr.pointattributelist );
1147  if(tr.pointmarkerlist != nullptr) free(tr.pointmarkerlist );
1148  if(tr.trianglelist != nullptr) free(tr.trianglelist );
1149  if(tr.triangleattributelist != nullptr) free(tr.triangleattributelist );
1150  if(tr.trianglearealist != nullptr) free(tr.trianglearealist );
1151  if(tr.neighborlist != nullptr) free(tr.neighborlist );
1152  if(tr.segmentlist != nullptr) free(tr.segmentlist );
1153  if(tr.segmentmarkerlist != nullptr) free(tr.segmentmarkerlist );
1154  if(tr.holelist != nullptr) free(tr.holelist );
1155  if(tr.regionlist != nullptr) free(tr.regionlist );
1156  if(tr.edgelist != nullptr) free(tr.edgelist );
1157  if(tr.edgemarkerlist != nullptr) free(tr.edgemarkerlist );
1158  if(tr.normlist != nullptr) free(tr.normlist );
1159  };
1163 
1164 
1168 
1169 
1173 
1174 
1178 
1181 
1182 
1184 
1185 }; // Class TriGenPFEMModeler
1186 
1188 
1191 
1192 
1196 
1197 
1199 inline std::istream& operator >> (std::istream& rIStream,
1200  TriGenPFEMModeler& rThis);
1201 
1203 inline std::ostream& operator << (std::ostream& rOStream,
1204  const TriGenPFEMModeler& rThis)
1205 {
1206  rThis.PrintInfo(rOStream);
1207  rOStream << std::endl;
1208  rThis.PrintData(rOStream);
1209 
1210  return rOStream;
1211 }
1213 
1214 
1215 } // namespace Kratos.
1216 
1217 #endif // KRATOS_TRIGEN_PFEM_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.
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
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
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_pfem_refine.h:64
Node::Pointer PointPointerType
Definition: trigen_pfem_refine.h:72
virtual void PrintInfo(std::ostream &rOStream) const
Print information about this object.
Definition: trigen_pfem_refine.h:578
Node PointType
Definition: trigen_pfem_refine.h:71
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: trigen_pfem_refine.h:581
Tree< KDTreePartition< BucketType > > KdtreeType
Definition: trigen_pfem_refine.h:78
PointVector::iterator PointIterator
Definition: trigen_pfem_refine.h:74
TriGenPFEMModeler()
Default constructor.
Definition: trigen_pfem_refine.h:85
Bucket< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType
Definition: trigen_pfem_refine.h:77
std::vector< double > DistanceVector
Definition: trigen_pfem_refine.h:75
virtual std::string Info() const
Turn back information as a string.
Definition: trigen_pfem_refine.h:572
KRATOS_CLASS_POINTER_DEFINITION(TriGenPFEMModeler)
Pointer definition of TriGenModeler.
virtual ~TriGenPFEMModeler()=default
Destructor.
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_pfem_refine.h:108
std::vector< PointType::Pointer > PointVector
Definition: trigen_pfem_refine.h:73
std::vector< double >::iterator DistanceIterator
Definition: trigen_pfem_refine.h:76
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
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