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.
tetgen_pfem_refine_vms.h
Go to the documentation of this file.
1 //
2 // Project Name: Kratos
3 // Last Modified by: $Author: anonymous $
4 // Date: $Date: 2009-01-15 14:50:34 $
5 // Revision: $Revision: 1.8 $
6 //
7 
8 
9 
10 
11 #if !defined(KRATOS_TETGEN_PFEM_MODELER_VMS_H_INCLUDED )
12 #define KRATOS_TETGEN_PFEM_MODELER_VMS_H_INCLUDED
13 
14 
15 
16 // System includes
17 #include <string>
18 #include <iostream>
19 #include <cstdlib>
20 #include <boost/timer.hpp>
21 
22 
23 
24 #include "tetgen.h" // Defined tetgenio, tetrahedralize().
25 
26 // Project includes
27 #include "includes/define.h"
28 #include "utilities/geometry_utilities.h"
29 #include "includes/model_part.h"
34 
36 //#include "containers/bucket.h"
37 //#include "containers/kd_tree.h"
38 //#include "external_includes/trigen_refine.h"
39 #include "tetgen_pfem_refine.h"
40 
41 namespace Kratos
42 {
43 
44 
47 
51 
55 
59 
63 
65 
68  {
69  public:
72 
75 
79 
82  mJ(ZeroMatrix(3,3)), //local jacobian
83  mJinv(ZeroMatrix(3,3)), //inverse jacobian
84  mc(ZeroVector(3)), //dimension = number of nodes
85  mRhs(ZeroVector(3)){} //dimension = number of nodes
86 
88  ~TetGenPfemModelerVms() override= default;
89 
90 
94 
95 
99 
100 
101  //*******************************************************************************************
102  //*******************************************************************************************
104  ModelPart& ThisModelPart ,
105  Element const& rReferenceElement,
106  Condition const& rReferenceBoundaryCondition,
107  EntitiesEraseProcess<Node>& node_erase, bool rem_nodes = true, bool add_nodes=true,
108  double alpha_param = 1.4, double h_factor=0.5)
109  {
110 
111  KRATOS_TRY
112  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FREE_SURFACE)==false )
113  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_FREE_SURFACE---- variable!!!!!! ERROR","");
114  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_STRUCTURE)==false )
115  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_STRUCTURE---- variable!!!!!! ERROR","");
116  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_BOUNDARY)==false )
117  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_BOUNDARY---- variable!!!!!! ERROR","");
118  if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FLUID)==false )
119  KRATOS_THROW_ERROR(std::logic_error,"Add ----IS_FLUID---- variable!!!!!! ERROR","");
120 
121  KRATOS_WATCH(" ENTERED TETGENMESHSUITE PFEM of Meshing Application")
122 
123  //clearing elements
124 
125  ThisModelPart.Elements().clear();
126  ThisModelPart.Conditions().clear();
127 
128  boost::timer auxiliary;
130  typedef Node PointType;
131  typedef Node::Pointer PointPointerType;
132  //typedef PointerVector<PointType> PointVector;
133  typedef std::vector<PointType::Pointer> PointVector;
134  typedef PointVector::iterator PointIterator;
135  typedef std::vector<double> DistanceVector;
136  typedef std::vector<double>::iterator DistanceIterator;
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
145 
146 
147  //*************
148  // DynamicBins;
149  typedef Tree< KDTreePartition<BucketType> > kd_tree; //Kdtree;
150  //typedef Tree< StaticBins > Bin; //Binstree;
151  unsigned int bucket_size = 20;
152 
153  //performing the interpolation - all of the nodes in this list will be preserved
154  unsigned int max_results = 50;
155  //PointerVector<PointType> res(max_results);
156  //NodeIterator res(max_results);
157  PointVector res(max_results);
158  DistanceVector res_distances(max_results);
159  Node work_point(0,0.0,0.0,0.0);
160  KRATOS_WATCH(h_factor)
161  //if the remove_node switch is activated, we check if the nodes got too close
162  if (rem_nodes==true)
163  {
164 
165  PointVector list_of_nodes;
166  list_of_nodes.reserve(ThisModelPart.Nodes().size());
167  for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
168  {
169  (list_of_nodes).push_back(*(i_node.base()));
170  //(list_of_nodes).push_back(i_node.base());
171  }
172 
173  kd_tree nodes_tree1(list_of_nodes.begin(),list_of_nodes.end(), bucket_size);
174  //std::cout<<nodes_tree2<<std::endl;
175 
176  unsigned int n_points_in_radius;
177  //radius means the distance, closer than which no node shall be allowd. if closer -> mark for erasing
178  double radius;
179 
180  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in != ThisModelPart.NodesEnd(); in++)
181  {
182  radius=h_factor*in->FastGetSolutionStepValue(NODAL_H);
183 
184  work_point[0]=in->X();
185  work_point[1]=in->Y();
186  work_point[2]=in->Z();
187 
188  n_points_in_radius = nodes_tree1.SearchInRadius(work_point, radius, res.begin(),res_distances.begin(), max_results);
189  if (n_points_in_radius>1)
190  {
191  if (in->FastGetSolutionStepValue(IS_BOUNDARY)==0.0 && in->FastGetSolutionStepValue(IS_STRUCTURE)==0.0)
192  {
193  //look if we are already erasing any of the other nodes
194  double erased_nodes = 0;
195  for(auto i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
196  erased_nodes += in->Is(TO_ERASE);
197 
198  if( erased_nodes < 1) //we cancel the node if no other nodes are being erased
199  in->Set(TO_ERASE,true);
200  }
201  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)
202  {
203  //here we loop over the neighbouring nodes and if there are nodes
204  //with IS_BOUNDARY=1 which are closer than 0.2*nodal_h from our we remove the node we are considering
205  unsigned int k = 0;
206  unsigned int counter = 0;
207  for(auto i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
208  {
209  if ( (*i)->FastGetSolutionStepValue(IS_BOUNDARY,1)==1.0 && res_distances[k] < 0.2*radius && res_distances[k] > 0.0 )
210  {
211  // KRATOS_WATCH( res_distances[k] );
212  counter += 1;
213  }
214  k++;
215  }
216  if(counter > 0)
217  in->Set(TO_ERASE,true);
218  }
219  }
220 
221  }
222 
223 
224  node_erase.Execute();
225  }
226 
230 
231 
232  tetgenio in, out, in2, outnew;
233 
234  // All indices start from 1.
235  in.firstnumber = 1;
236 
237  in.numberofpoints = ThisModelPart.Nodes().size();
238  in.pointlist = new REAL[in.numberofpoints * 3];
239 
240  //writing the point coordinates in a vector
241  ModelPart::NodesContainerType::iterator nodes_begin = ThisModelPart.NodesBegin();
242 
243  //reorder node Ids
244  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
245  {
246  (nodes_begin + i)->SetId(i+1);
247 // (nodes_begin + i)->Id() = i+1;
248  }
249 
250  //give the corrdinates to the mesher
251  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
252  {
253  int base = i*3;
254  in.pointlist[base] = (nodes_begin + i)->X();
255  in.pointlist[base+1] = (nodes_begin + i)->Y();
256  in.pointlist[base+2] = (nodes_begin + i)->Z();
257  }
258  char tetgen_options[] = "SQJ";
259 
260 
261  tetrahedralize(tetgen_options, &in, &out); //with option to remove slivers
262 
263 
264  double first_part_time = auxiliary.elapsed();
265  std::cout << "mesh generation time = " << first_part_time << std::endl;
266 
267  //generate Kratos Tetrahedra3D4
268  int el_number = out.numberoftetrahedra;
269 
270  boost::timer alpha_shape_time;
271  //calculate r , center for all of the elements
272  std::vector<int> preserved_list(el_number);
273  array_1d<double,3> x1,x2,x3,x4, xc;
274  int point_base;
275  int number_of_preserved_elems = 0;
276  for(int el = 0; el< el_number; el++)
277  {
278  int base = el * 4;
279 
280  //coordinates
281  point_base = (out.tetrahedronlist[base] - 1)*3;
282  x1[0] = out.pointlist[point_base];
283  x1[1] = out.pointlist[point_base+1];
284  x1[2] = out.pointlist[point_base+2];
285 
286  point_base = (out.tetrahedronlist[base+1] - 1)*3;
287  x2[0] = out.pointlist[point_base];
288  x2[1] = out.pointlist[point_base+1];
289  x2[2] = out.pointlist[point_base+2];
290 
291  point_base = (out.tetrahedronlist[base+2] - 1)*3;
292  x3[0] = out.pointlist[point_base];
293  x3[1] = out.pointlist[point_base+1];
294  x3[2] = out.pointlist[point_base+2];
295 
296  point_base = (out.tetrahedronlist[base+3] - 1)*3;
297  x4[0] = out.pointlist[point_base];
298  x4[1] = out.pointlist[point_base+1];
299  x4[2] = out.pointlist[point_base+2];
300 
301  //calculate the geometrical data
302  //degenerate elements are given a very high radius
303  double geometrical_hmin, geometrical_hmax;
304  double radius;
305  double vol;
306  CalculateElementData(x1,x2,x3,x4,vol,xc,radius,geometrical_hmin,geometrical_hmax);
307 
308  //calculate the prescribed h
309  double prescribed_h = (nodes_begin + out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(NODAL_H);
310  prescribed_h += (nodes_begin + out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(NODAL_H);
311  prescribed_h += (nodes_begin + out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(NODAL_H);
312  prescribed_h += (nodes_begin + out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(NODAL_H);
313  prescribed_h *= 0.25;
314 
315  //check the number of nodes on the wall
316  int nb = int( (nodes_begin + out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
317  nb += int( (nodes_begin + out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
318  nb += int( (nodes_begin + out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
319  nb += int((nodes_begin + out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
320 
321  //check the number of nodes of bo node_erase.Execute();undary
322  int nfs = int( (nodes_begin + out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
323  nfs += int( (nodes_begin + out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
324  nfs += int( (nodes_begin + out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
325  nfs += int((nodes_begin + out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
326 
327  //check the number of nodes of boundary
328  int nfluid = int( (nodes_begin + out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_FLUID) );
329  nfluid += int( (nodes_begin + out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_FLUID) );
330  nfluid += int( (nodes_begin + out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_FLUID) );
331  nfluid += int((nodes_begin + out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_FLUID) );
332 
333  int n_lag = 0;
334  if (ThisModelPart.GetNodalSolutionStepVariablesList().Has(IS_LAGRANGIAN_INLET))
335  {
336  n_lag = int( (nodes_begin + out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
337  n_lag += int( (nodes_begin + out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
338  n_lag += int( (nodes_begin + out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
339  n_lag += int((nodes_begin + out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
340  }
341 
342 
343  //cases:
344  //4 nodes on the wall - elminate
345  // at least one node of boundary OR at least one node NOT of fluid --> pass alpha shape
346  /*
347  if(nb == 4) // 4 nodes on the wall
348  preserved_list[el] = false;
349  else if (nboundary != 0 || nfluid != 4) //close to the free surface or external
350  {
351  if( radius < prescribed_h * alpha_param && //alpha shape says to preserve
352  nb!=4) //the nodes are not all on the boundary
353  {
354  preserved_list[el] = true; //preserve!!
355  number_of_preserved_elems += 1;
356  }
357  }
358  else
359  {
360  preserved_list[el] = true;
361  number_of_preserved_elems += 1;
362  }
363  */
364  if(nb == 4) // 4 nodes on the wall
365  {
366  preserved_list[el] = false;
367 
368  }
369  else
370  {
371  //if (nfs != 0 || nfluid != 4) //close to the free surface or external
372  if (n_lag >= 2 && nb==3 ) //close to the free surface or external
373  {
374  if( radius < prescribed_h * 1.5*alpha_param ) //alpha shape says to preserve
375  //the nodes are not all on the boundary
376  {
377  preserved_list[el] = true; //preserve!!
378  number_of_preserved_elems += 1;
379  }
380  else
381  {
382  preserved_list[el] = false;
383 
384  }
385  }
386  else if (nfs != 0 ) //close to the free surface or external
387  {
388  if( radius < prescribed_h * alpha_param ) //alpha shape says to preserve
389  //the nodes are not all on the boundary
390  {
391  preserved_list[el] = true; //preserve!!
392  number_of_preserved_elems += 1;
393  }
394  else
395  {
396  preserved_list[el] = false;
397 
398  }
399  }
400 
401 
402  else if (nfluid < 3.9) //close to the free surface or external
403  {
404  if( radius < prescribed_h * alpha_param ) //alpha shape says to preserve
405  //the nodes are not all on the boundary
406  {
407  preserved_list[el] = true; //preserve!!
408  number_of_preserved_elems += 1;
409  }
410  else
411  {
412  preserved_list[el] = false;
413 
414  }
415  }
416 
417  else //internal elements should be preserved as much as possible not to create holes
418  {
419  if( radius < prescribed_h * alpha_param * 5.0 )
420  {
421 //std::cout << "element not deleted" <<std::endl;
422  preserved_list[el] = true; //preserve!!
423  number_of_preserved_elems += 1;
424  }
425  else
426  {
427 //std::cout << "sliver removed" << std::endl;
428  preserved_list[el] = false;
429 
430  }
431 // preserved_list[el] = true;
432 // number_of_preserved_elems += 1;
433  }
434  }
435 
436  }
437  std::cout << "time for passing alpha shape" << alpha_shape_time.elapsed() << std::endl;
438 
439  //freeing unnecessary memory
440  in.deinitialize();
441  in.initialize();
442 
443  //setting the new new ids in a aux vector
444 // tetgenio in2;
445 
446  in2.firstnumber = 1;
447  in2.numberofpoints = ThisModelPart.Nodes().size();
448  in2.pointlist = new REAL[in2.numberofpoints * 3];
449 
450  for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
451  {
452  int base = i*3;
453  in2.pointlist[base] = (nodes_begin + i)->X();
454  in2.pointlist[base+1] = (nodes_begin + i)->Y();
455  in2.pointlist[base+2] = (nodes_begin + i)->Z();
456  }
457  std::cout << "qui" << std::endl;
458  in2.numberoftetrahedra = number_of_preserved_elems;
459  in2.tetrahedronlist = new int[in2.numberoftetrahedra * 4];
460  in2.tetrahedronvolumelist = new double[in2.numberoftetrahedra];
461 
462  int counter = 0;
463  for(int el = 0; el< el_number; el++)
464  {
465  if( preserved_list[el] )
466  {
467  //saving the compact element list
468  int new_base = counter*4;
469  int old_base = el*4;
470  in2.tetrahedronlist[new_base] = out.tetrahedronlist[old_base];
471  in2.tetrahedronlist[new_base+1] = out.tetrahedronlist[old_base+1];
472  in2.tetrahedronlist[new_base+2] = out.tetrahedronlist[old_base+2];
473  in2.tetrahedronlist[new_base+3] = out.tetrahedronlist[old_base+3];
474 
475  //calculate the prescribed h
476  double prescribed_h = (nodes_begin + out.tetrahedronlist[old_base]-1)->FastGetSolutionStepValue(NODAL_H);
477  prescribed_h += (nodes_begin + out.tetrahedronlist[old_base+1]-1)->FastGetSolutionStepValue(NODAL_H);
478  prescribed_h += (nodes_begin + out.tetrahedronlist[old_base+2]-1)->FastGetSolutionStepValue(NODAL_H);
479  prescribed_h += (nodes_begin + out.tetrahedronlist[old_base+3]-1)->FastGetSolutionStepValue(NODAL_H);
480  prescribed_h *= 0.25;
481  //if h is the height of a perfect tetrahedra, the edge size is edge = sqrt(3/2) h
482  //filling in the list of "IDEAL" tetrahedron volumes=1/12 * (edge)^3 * sqrt(2)~0.11785* h^3=
483  //0.2165063509*h^3
484 
485 
486  //chapuza to be checked - I just didnt want so much of refinement,,,,
487  in2.tetrahedronvolumelist[counter] = 2000.0*0.217*prescribed_h*prescribed_h*prescribed_h;
488 
489 
490  //in2.tetrahedronvolumelist[counter] = 0.217*prescribed_h*prescribed_h*prescribed_h;
491  //in2.tetrahedronvolumelist[counter] = 0.0004;
492  //KRATOS_WATCH(in2.tetrahedronvolumelist[counter])
493  counter += 1;
494  }
495 
496  }
497 
498  //NOW WE SHALL IDENTIFY FLYING NODES
499  /*
500  std::vector<double> aux_before_ref(ThisModelPart.Nodes().size());
501 
502  for (unsigned int i=0; i<aux_before_ref.size();i++)
503  {
504  aux_before_ref[i]=0.0;
505  }
506 
507  for(int el = 0; el< el_number; el++)
508  {
509  int old_base=el*4;
510  if (preserved_list[el]==true)
511  {
512 
513  //we add a non-zero value for the node of the element that has a non-zero value
514  aux_before_ref[out.tetrahedronlist[old_base]-1]+=1;
515  aux_before_ref[out.tetrahedronlist[old_base+1]-1]+=1;
516  aux_before_ref[out.tetrahedronlist[old_base+2]-1]+=1;
517  aux_before_ref[out.tetrahedronlist[old_base+3]-1]+=1;
518  }
519  }
520  */
521  //creating a new mesh
522  boost::timer mesh_recreation_time;
523 
524  //freeing unnecessary memory
525  out.deinitialize();
526  out.initialize();
527 
528 
529  //HERE WE ADD THE VOLUME CONSTRAINT -the desired volume of the equidistant tertrahedra
530  //based upon average nodal_h (that is the "a" switch
531  //char regeneration_options[] = "rQJYq1.8anS", "rQJYq1.4/20nS";
532  if (add_nodes==true)
533  {
534  char mesh_regen_opts[] = "rQJYq1.4/20nS";
535  tetrahedralize(mesh_regen_opts, &in2, &outnew);
536  KRATOS_WATCH("Adaptive remeshing executed")
537  }
538  else
539  {
540  char mesh_regen_opts[] = "rQJYnS";
541  tetrahedralize(mesh_regen_opts, &in2, &outnew);
542  KRATOS_WATCH("Non-Adaptive remeshing executed")
543  }
544 
545  //q - creates quality mesh, with the default radius-edge ratio set to 2.0
546 
547  std::cout << "mesh recreation time" << mesh_recreation_time.elapsed() << std::endl;
548 
549 
550  //PAVEL
551 
552  //putting the new nodes in a spatial container
553  //PointerVector< Element >& neighb
554 
555  /*
556  typedef Node PointType;
557  typedef Node::Pointer PointPointerType;
558  typedef std::vector<PointType::Pointer> PointVector;
559  typedef PointVector::iterator PointIterator;
560  typedef std::vector<double> DistanceVector;
561  typedef std::vector<double>::iterator DistanceIterator;
562 
563  int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
564 
565  // bucket types
566  typedef Bucket<3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType;
567  typedef Bins< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > StaticBins;
568 
569  // DynamicBins;
570  typedef Tree< KDTreePartition<BucketType> > kd_tree; //Kdtree;
571  //typedef Tree< StaticBins > Bin; //Binstree;
572  */
573  PointVector list_of_new_nodes;
574 
575  Node::DofsContainerType& reference_dofs = (ThisModelPart.NodesBegin())->GetDofs();
576 
577  int n_points_before_refinement = in2.numberofpoints;
578  //if the refinement was performed, we need to add it to the model part.
579  if (outnew.numberofpoints>n_points_before_refinement)
580  {
581  for(int i = n_points_before_refinement; i<outnew.numberofpoints;i++)
582  {
583  int id=i+1;
584  int base = i*3;
585  double& x= outnew.pointlist[base];
586  double& y= outnew.pointlist[base+1];
587  double& z= outnew.pointlist[base+2];
588 
589  Node::Pointer pnode = ThisModelPart.CreateNewNode(id,x,y,z);
590 
591  //putting the new node also in an auxiliary list
592  //KRATOS_WATCH("adding nodes to list")
593  list_of_new_nodes.push_back( pnode );
594 
595  //std::cout << "new node id = " << pnode->Id() << std::endl;
596  //generating the dofs
597  for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
598  {
599  Node::DofType &rDof = **iii;
600  Node::DofType::Pointer p_new_dof = pnode->pAddDof( rDof );
601 
602  (p_new_dof)->FreeDof();
603  }
604 
605  }
606  }
607 
608  std::cout << "During refinement we added " << outnew.numberofpoints-n_points_before_refinement<< "nodes " <<std::endl;
609 
610 
611  bucket_size = 50;
612 
613  //performing the interpolation - all of the nodes in this list will be preserved
614  max_results = 800;
615  PointVector results(max_results);
616  DistanceVector results_distances(max_results);
618  //int data_size = ThisModelPart.GetNodalSolutionStepTotalDataSize();
619 
620 
621  //double* work_array;
622  //Node work_point(0,0.0,0.0,0.0);
623 
624 
625 
626  //NOW WE SHALL IDENTIFY LONELY BUT NOT FLYING NODES (compare with the flying nodes identification above (after first remeshing step))
627  /*
628  std::vector<double> aux_after_ref(outnew.numberofpoints);
629 
630  for (unsigned int i=0; i<aux_after_ref.size();i++)
631  {
632  aux_after_ref[i]=0.0;
633  }
634  int el_number_ref= outnew.numberoftetrahedra;
635 
636  for(int el = 0; el< el_number_ref; el++)
637  {
638  int base=el*4;
639 
640 
641  //we add a non-zero value for the node of the element that has a non-zero value
642  aux_after_ref[outnew.tetrahedronlist[base]-1]+=1;
643  aux_after_ref[outnew.tetrahedronlist[base+1]-1]+=1;
644  aux_after_ref[outnew.tetrahedronlist[base+2]-1]+=1;
645  aux_after_ref[outnew.tetrahedronlist[base+3]-1]+=1;
646 
647  }
648  */
649 
650  if(outnew.numberofpoints-n_points_before_refinement > 0) //if we added points
651  {
652  kd_tree nodes_tree2(list_of_new_nodes.begin(),list_of_new_nodes.end(),bucket_size);
653  //std::cout<<nodes_tree2<<std::endl;
654  nodes_begin = ThisModelPart.NodesBegin();
655 
656  for(int el = 0; el< in2.numberoftetrahedra; el++)
657  //for(unsigned int el = 0; el< outnew.numberoftetrahedra; el++)
658  {
659 
660  int base = el * 4;
661  //coordinates
662 
663  point_base = (in2.tetrahedronlist[base] - 1)*3;
664  x1[0] = in2.pointlist[point_base];
665  x1[1] = in2.pointlist[point_base+1];
666  x1[2] = in2.pointlist[point_base+2];
667 
668  point_base = (in2.tetrahedronlist[base+1] - 1)*3;
669  x2[0] = in2.pointlist[point_base];
670  x2[1] = in2.pointlist[point_base+1];
671  x2[2] = in2.pointlist[point_base+2];
672 
673  point_base = (in2.tetrahedronlist[base+2] - 1)*3;
674  x3[0] = in2.pointlist[point_base];
675  x3[1] = in2.pointlist[point_base+1];
676  x3[2] = in2.pointlist[point_base+2];
677 
678  point_base = (in2.tetrahedronlist[base+3] - 1)*3;
679  x4[0] = in2.pointlist[point_base];
680  x4[1] = in2.pointlist[point_base+1];
681  x4[2] = in2.pointlist[point_base+2];
682 
683  //calculate the geometrical data
684  //degenerate elements are given a very high radius
685  double geometrical_hmin, geometrical_hmax;
686  double radius;
687  double vol;
688 
690 
691  //it calculates the data necessary to perfom the search, like the element center etc., search radius
692  //and writes it to radius, xc etc etc
693 
694  CalculateElementData(x1,x2,x3,x4,vol,xc,radius,geometrical_hmin,geometrical_hmax);
695 
696  //find all of the nodes in a radius centered in xc
697 
698  std::size_t number_of_points_in_radius;
699  work_point.X() = xc[0]; work_point.Y() = xc[1]; work_point.Z() = xc[2];
700 
701  number_of_points_in_radius = nodes_tree2.SearchInRadius(work_point, radius, results.begin(),results_distances.begin(), max_results);
702 
703 
704  //for each of the nodes in radius find if it is inside the element or not
705  //if inside interpolate
706 
707 
709  *( (nodes_begin + in2.tetrahedronlist[base]-1).base() ),
710  *( (nodes_begin + in2.tetrahedronlist[base+1]-1).base() ),
711  *( (nodes_begin + in2.tetrahedronlist[base+2]-1).base() ),
712  *( (nodes_begin + in2.tetrahedronlist[base+3]-1).base() )
713  );
714 
715 
716  //KRATOS_WATCH(results.size())
717  for(auto it=results.begin(); it!=results.begin() + number_of_points_in_radius; it++)
718  {
719  bool is_inside=false;
720 
721  is_inside = CalculatePosition(x1[0],x1[1],x1[2],
722  x2[0],x2[1],x2[2],
723  x3[0],x3[1],x3[2],
724  x4[0],x4[1],x4[2],
725  (*it)->X(),(*it)->Y(), (*it)->Z(),N);
726 
727  if(is_inside == true)
728  {
729 
730  Interpolate( geom, N, step_data_size, *(it) );
731 
732  }
733 
734  }
735  }
736  }
737 
738  ThisModelPart.Elements().clear();
739  ThisModelPart.Conditions().clear();
740 
741  //set the coordinates to the original value
742  for(auto & list_of_new_node : list_of_new_nodes)
743  {
744  const array_1d<double,3>& disp = list_of_new_node->FastGetSolutionStepValue(DISPLACEMENT);
745  list_of_new_node->X0() = list_of_new_node->X() - disp[0];
746  list_of_new_node->Y0() = list_of_new_node->Y() - disp[1];
747  list_of_new_node->Z0() = list_of_new_node->Z() - disp[2];
748  }
749  //cleaning unnecessary data
750  in2.deinitialize();
751  in2.initialize();
752 
753  //***********************************************************************************
754  //***********************************************************************************
755  boost::timer adding_elems;
756  //add preserved elements to the kratos
757  Properties::Pointer properties = ThisModelPart.GetMesh().pGetProperties(1);
758  nodes_begin = ThisModelPart.NodesBegin();
759  (ThisModelPart.Elements()).reserve(outnew.numberoftetrahedra);
760 
761  for(int iii = 0; iii< outnew.numberoftetrahedra; iii++)
762  {
763  int id = iii + 1;
764  int base = iii * 4;
766  *( (nodes_begin + outnew.tetrahedronlist[base]-1).base() ),
767  *( (nodes_begin + outnew.tetrahedronlist[base+1]-1).base() ),
768  *( (nodes_begin + outnew.tetrahedronlist[base+2]-1).base() ),
769  *( (nodes_begin + outnew.tetrahedronlist[base+3]-1).base() )
770  );
771 
772 
773 #ifdef _DEBUG
774 ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
775  if( *(ModelNodes).find( outnew.tetrahedronlist[base]).base() == *(ThisModelPart.Nodes().end()).base() )
776  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
777  if( *(ModelNodes).find( outnew.tetrahedronlist[base+1]).base() == *(ThisModelPart.Nodes().end()).base() )
778  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
779  if( *(ModelNodes).find( outnew.tetrahedronlist[base+2]).base() == *(ThisModelPart.Nodes().end()).base() )
780  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
781  if( *(ModelNodes).find( outnew.tetrahedronlist[base+3]).base() == *(ThisModelPart.Nodes().end()).base() )
782  KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
783 #endif
784 
785  Element::Pointer p_element = rReferenceElement.Create(id, geom, properties);
786  (ThisModelPart.Elements()).push_back(p_element);
787 
788  }
789  std::cout << "time for adding elems" << adding_elems.elapsed() << std::endl;;
790  ThisModelPart.Elements().Sort();
791 
792  boost::timer adding_neighb;
793 // //filling the neighbour list
794  ModelPart::ElementsContainerType::const_iterator el_begin = ThisModelPart.ElementsBegin();
795  for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.ElementsBegin();
796  iii != ThisModelPart.ElementsEnd(); iii++)
797  {
798  //Geometry< Node >& geom = iii->GetGeometry();
799  int base = ( iii->Id() - 1 )*4;
800 
801  (iii->GetValue(NEIGHBOUR_ELEMENTS)).resize(4);
802  GlobalPointersVector< Element >& neighb = iii->GetValue(NEIGHBOUR_ELEMENTS);
803 
804  for(int i = 0; i<4; i++)
805  {
806  int index = outnew.neighborlist[base+i];
807  if(index > 0)
808  neighb(i) = GlobalPointer<Element>(&*(el_begin + index-1));
809  else
810  neighb(i) = Element::WeakPointer();
811  }
812  }
813  std::cout << "time for adding neigbours" << adding_neighb.elapsed() << std::endl;;
814 
815 
816 
817 
818 
819  //***********************************************************************************
820  //***********************************************************************************
821  //mark boundary nodes
822  //reset the boundary flag
823  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in!=ThisModelPart.NodesEnd(); in++)
824  {
825  in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
826  }
827 
828 
829  //***********************************************************************************
830  //***********************************************************************************
831  boost::timer adding_faces;
832 
833  (ThisModelPart.Conditions()).reserve(outnew.numberoftrifaces );
834 
835  //creating the faces
836  for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.ElementsBegin();
837  iii != ThisModelPart.ElementsEnd(); iii++)
838  {
839 
840  int base = ( iii->Id() - 1 )*4;
841 
842  //create boundary faces and mark the boundary nodes
843  //each face is opposite to the corresponding node number so
844  // 0 ----- 1 2 3
845  // 1 ----- 0 3 2
846  // 2 ----- 0 1 3
847  // 3 ----- 0 2 1
848 
849  //node 1
850 
851  //if( neighb(0).expired() );
852  if( outnew.neighborlist[base] == -1)
853  {
854  CreateBoundaryFace(1, 2, 3, ThisModelPart, 0, *(iii.base()), rReferenceBoundaryCondition, properties );
855  }
856  //if(neighb(1).expired() );
857  if( outnew.neighborlist[base+1] == -1)
858  {
859  CreateBoundaryFace(0,3,2, ThisModelPart, 1, *(iii.base()), rReferenceBoundaryCondition, properties );
860  }
861  if( outnew.neighborlist[base+2] == -1)
862  //if(neighb(2).expired() );
863  {
864  CreateBoundaryFace(0,1,3, ThisModelPart, 2, *(iii.base()), rReferenceBoundaryCondition, properties );
865  }
866  if( outnew.neighborlist[base+3] == -1)
867  //if(neighb(3).expired() );
868  {
869  CreateBoundaryFace(0,2,1, ThisModelPart, 3, *(iii.base()), rReferenceBoundaryCondition, properties );
870  }
871 
872  }
873  outnew.deinitialize();
874  outnew.initialize();
875  std::cout << "time for adding faces" << adding_faces.elapsed() << std::endl;;
876 
877 
878 
879  //here we remove lonely nodes that are not teh flying nodes, but are the lonely nodes inside water vol
880  /*
881  for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in!=ThisModelPart.NodesEnd(); in++)
882  {
883  //if this is a node added during refinement, but isnt contained in any element
884  if (aux_after_ref[in->GetId()]==0 && in->GetId()>aux_before_ref.size() && in->GetId()<aux_after_ref.size())
885  {
886  in->Set(TO_ERASE,true);
887  KRATOS_WATCH("This is that ugly ugly lonely interior node a666a. IT SHALL BE TERRRRMINATED!!!!!!")
888  KRATOS_WATCH(in->GetId())
889  }
890  //if this was an interior node after first step and became single due to derefinement - erase it
891  if (aux_after_ref[in->GetId()]==0 && in->GetId()<aux_before_ref.size() && aux_before_ref[in->GetId()]!=0)
892  {
893  in->Set(TO_ERASE,true);
894  KRATOS_WATCH("This is that ugly ugly lonely interior node. IT SHALL BE TERRRRMINATED!!!!!!")
895  }
896  }
897  */
898 
899 
900  double second_part_time = auxiliary.elapsed();
901  std::cout << "second part time = " << second_part_time - first_part_time << std::endl;
902 
903 
904  KRATOS_CATCH("")
905  }
906 
910 
911 
915 
916 
920 
922  std::string Info() const override{return "";}
923 
925  void PrintInfo(std::ostream& rOStream) const override{}
926 
928  void PrintData(std::ostream& rOStream) const override{}
929 
930 
934 
935 
937 
938  protected:
941 
942 
946 
947 
951 
952 
956 
957 
961 
962 
966 
967 
971 
972 
974 
975  private:
978 
979 
983  boost::numeric::ublas::bounded_matrix<double,3,3> mJ; //local jacobian
984  boost::numeric::ublas::bounded_matrix<double,3,3> mJinv; //inverse jacobian
985  array_1d<double,3> mc; //center pos
986  array_1d<double,3> mRhs; //center pos
987 
988 
989  void CreateBoundaryFace(const int& i1, const int& i2, const int& i3, ModelPart& ThisModelPart, const int& outer_node_id, Element::Pointer origin_element, Condition const& rReferenceBoundaryCondition, Properties::Pointer properties)
990  {
991  KRATOS_TRY
992 
993  Geometry<Node >& geom = origin_element->GetGeometry();
994  //mark the nodes as free surface
995  geom[i1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
996  geom[i2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
997  geom[i3].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
998 
999  //generate a face condition
1001  temp.reserve(3);
1002  temp.push_back(geom(i1));
1003  temp.push_back(geom(i2));
1004  temp.push_back(geom(i3));
1006  //Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Triangle3D< Node >(temp) );
1007  int id = (origin_element->Id()-1)*4;
1008  Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp, properties);
1009 
1010  //assigning the neighbour node
1011  (p_cond->GetValue(NEIGHBOUR_NODES)).clear();
1012  (p_cond->GetValue(NEIGHBOUR_NODES)).push_back( Node::WeakPointer( geom(outer_node_id) ) );
1013  (p_cond->GetValue(NEIGHBOUR_ELEMENTS)).clear();
1014  (p_cond->GetValue(NEIGHBOUR_ELEMENTS)).push_back( Element::WeakPointer( origin_element ) );
1015  ThisModelPart.Conditions().push_back(p_cond);
1016  KRATOS_CATCH("")
1017 
1018  }
1019 
1022  void Interpolate( Tetrahedra3D4<Node >& geom, const array_1d<double,4>& N,
1023  unsigned int step_data_size,
1024  Node::Pointer pnode)
1025  {
1026  unsigned int buffer_size = pnode->GetBufferSize();
1027 
1028 
1029  for(unsigned int step = 0; step<buffer_size; step++)
1030  {
1031 
1032  //getting the data of the solution step
1033  double* step_data = (pnode)->SolutionStepData().Data(step);
1034 
1035 
1036  double* node0_data = geom[0].SolutionStepData().Data(step);
1037  double* node1_data = geom[1].SolutionStepData().Data(step);
1038  double* node2_data = geom[2].SolutionStepData().Data(step);
1039  double* node3_data = geom[3].SolutionStepData().Data(step);
1040 
1041  //copying this data in the position of the vector we are interested in
1042  for(unsigned int j= 0; j<step_data_size; j++)
1043  {
1044 
1045  step_data[j] = N[0]*node0_data[j] + N[1]*node1_data[j] + N[2]*node2_data[j] + N[3]*node3_data[j];
1046 
1047 
1048  }
1049  }
1050  //now we assure that the flag variables are set coorect!! since we add nodes inside of the fluid volume only
1051  //we manually reset the IS_BOUNDARY, IS_FLUID, IS_STRUCTURE, IS_FREE_SURFACE values in a right way
1052  //not to have values, like 0.33 0.66 resulting if we would have been interpolating them in the same way
1053  //as the normal variables, like Velocity etc
1054 
1055  pnode->FastGetSolutionStepValue(IS_BOUNDARY)=0.0;
1056  pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1057  pnode->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET)=0.0;
1058  pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=0.0;
1059  pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1060  }
1061 
1062  inline double CalculateVol( const double x0, const double y0, const double z0,
1063  const double x1, const double y1, const double z1,
1064  const double x2, const double y2, const double z2,
1065  const double x3, const double y3, const double z3
1066  )
1067  {
1068  double x10 = x1 - x0;
1069  double y10 = y1 - y0;
1070  double z10 = z1 - z0;
1071 
1072  double x20 = x2 - x0;
1073  double y20 = y2 - y0;
1074  double z20 = z2 - z0;
1075 
1076  double x30 = x3 - x0;
1077  double y30 = y3 - y0;
1078  double z30 = z3 - z0;
1079 
1080  double detJ = x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30;
1081  return detJ*0.1666666666666666666667;
1082 
1083  //return 0.5*( (x1-x0)*(y2-y0)- (y1-y0)*(x2-x0) );
1084  }
1085 
1086  inline bool CalculatePosition( const double x0, const double y0, const double z0,
1087  const double x1, const double y1, const double z1,
1088  const double x2, const double y2, const double z2,
1089  const double x3, const double y3, const double z3,
1090  const double xc, const double yc, const double zc,
1091  array_1d<double,4>& N
1092  )
1093  {
1094 
1095 
1096  double vol = CalculateVol(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3);
1097 
1098  double inv_vol = 0.0;
1099  if(vol < 1e-26)
1100  {
1101  //KRATOS_THROW_ERROR(std::logic_error,"element with zero vol found","");
1102  KRATOS_WATCH("WARNING!!!!!!!!!! YOU HAVE A SLIVER HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
1103  return false;
1104  }
1105  else
1106  {
1107  inv_vol = 1.0 / vol;
1108  N[0] = CalculateVol(x1,y1,z1,x3,y3,z3,x2,y2,z2,xc,yc,zc) * inv_vol;
1109  N[1] = CalculateVol(x3,y3,z3,x0,y0,z0,x2,y2,z2,xc,yc,zc) * inv_vol;
1110  N[2] = CalculateVol(x3,y3,z3,x1,y1,z1,x0,y0,z0,xc,yc,zc) * inv_vol;
1111  N[3] = CalculateVol(x0,y0,z0,x1,y1,z1,x2,y2,z2,xc,yc,zc) * inv_vol;
1112 
1113 
1114  if(N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[3] >=0.0 && N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0 && N[3] <=1.0)
1115  {
1116  //if the xc yc zc is inside the tetrahedron return true
1117  return true;
1118  }
1119  else
1120  return false;
1121  }
1122 
1123 
1124  }
1125  //***************************************
1126  //***************************************
1127  inline void CalculateCenterAndSearchRadius(const double x0, const double y0, const double z0,
1128  const double x1, const double y1, const double z1,
1129  const double x2, const double y2, const double z2,
1130  const double x3, const double y3, const double z3,
1131  double& xc, double& yc, double& zc, double& R
1132  )
1133  {
1134  xc = 0.25*(x0+x1+x2+x3);
1135  yc = 0.25*(y0+y1+y2+y3);
1136  zc = 0.25*(z0+z1+z2+z3);
1137 
1138  double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0) + (zc-z0)*(zc-z0);
1139  double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1) + (zc-z1)*(zc-z1);
1140  double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2) + (zc-z2)*(zc-z2);
1141  double R4 = (xc-x3)*(xc-x3) + (yc-y3)*(yc-y3) + (zc-z3)*(zc-z3);
1142 
1143  R = R1;
1144  if(R2 > R) R = R2;
1145  if(R3 > R) R = R3;
1146  if(R4 > R) R = R4;
1147 
1148  R = sqrt(R);
1149  }
1150  //*****************************************************************************************
1151  //*****************************************************************************************
1152  void CalculateElementData(const array_1d<double,3>& c1,
1153  const array_1d<double,3>& c2,
1154  const array_1d<double,3>& c3,
1155  const array_1d<double,3>& c4,
1156  double& volume,
1157  array_1d<double,3>& xc,
1158  double& radius,
1159  double& hmin,
1160  double& hmax
1161  )
1162  {
1163  KRATOS_TRY
1164  const double x0 = c1[0]; const double y0 = c1[1]; const double z0 = c1[2];
1165  const double x1 = c2[0]; const double y1 = c2[1]; const double z1 = c2[2];
1166  const double x2 = c3[0]; const double y2 = c3[1]; const double z2 = c3[2];
1167  const double x3 = c4[0]; const double y3 = c4[1]; const double z3 = c4[2];
1168 
1169 // KRATOS_WATCH("111111111111111111111");
1170  //calculate min side lenght
1171  //(use xc as a auxiliary vector) !!!!!!!!!!!!
1172  double aux;
1173  noalias(xc) = c4; noalias(xc)-=c1; aux = inner_prod(xc,xc); hmin = aux; hmax = aux;
1174  noalias(xc) = c3; noalias(xc)-=c1; aux = inner_prod(xc,xc); if(aux<hmin) hmin = aux; else if(aux>hmax) hmax = aux;
1175  noalias(xc) = c2; noalias(xc)-=c1; aux = inner_prod(xc,xc); if(aux<hmin) hmin = aux; else if(aux>hmax) hmax = aux;
1176  noalias(xc) = c4; noalias(xc)-=c2; aux = inner_prod(xc,xc); if(aux<hmin) hmin = aux; else if(aux>hmax) hmax = aux;
1177  noalias(xc) = c3; noalias(xc)-=c2; aux = inner_prod(xc,xc); if(aux<hmin) hmin = aux; else if(aux>hmax) hmax = aux;
1178  noalias(xc) = c4; noalias(xc)-=c3; aux = inner_prod(xc,xc); if(aux<hmin) hmin = aux; else if(aux>hmax) hmax = aux;
1179  hmin = sqrt(hmin);
1180  hmax = sqrt(hmax);
1181 
1182  //calculation of the jacobian
1183  mJ(0,0) = x1-x0; mJ(0,1) = y1-y0; mJ(0,2) = z1-z0;
1184  mJ(1,0) = x2-x0; mJ(1,1) = y2-y0; mJ(1,2) = z2-z0;
1185  mJ(2,0) = x3-x0; mJ(2,1) = y3-y0; mJ(2,2) = z3-z0;
1186 // KRATOS_WATCH("33333333333333333333333");
1187 
1188  //inverse of the jacobian
1189  //first column
1190  mJinv(0,0) = mJ(1,1)*mJ(2,2) - mJ(1,2)*mJ(2,1);
1191  mJinv(1,0) = -mJ(1,0)*mJ(2,2) + mJ(1,2)*mJ(2,0);
1192  mJinv(2,0) = mJ(1,0)*mJ(2,1) - mJ(1,1)*mJ(2,0);
1193  //second column
1194  mJinv(0,1) = -mJ(0,1)*mJ(2,2) + mJ(0,2)*mJ(2,1);
1195  mJinv(1,1) = mJ(0,0)*mJ(2,2) - mJ(0,2)*mJ(2,0);
1196  mJinv(2,1) = -mJ(0,0)*mJ(2,1) + mJ(0,1)*mJ(2,0);
1197  //third column
1198  mJinv(0,2) = mJ(0,1)*mJ(1,2) - mJ(0,2)*mJ(1,1);
1199  mJinv(1,2) = -mJ(0,0)*mJ(1,2) + mJ(0,2)*mJ(1,0);
1200  mJinv(2,2) = mJ(0,0)*mJ(1,1) - mJ(0,1)*mJ(1,0);
1201  //calculation of determinant (of the input matrix)
1202 
1203 // KRATOS_WATCH("44444444444444444444444444");
1204  double detJ = mJ(0,0)*mJinv(0,0)
1205  + mJ(0,1)*mJinv(1,0)
1206  + mJ(0,2)*mJinv(2,0);
1207 
1208  volume = detJ * 0.16666666667;
1209 // KRATOS_WATCH("55555555555555555555555");
1210 
1211  if(volume < 1e-3 * hmax*hmax*hmax) //this is a sliver and we will remove it
1212  {
1213  //KRATOS_WATCH("666666666666666666666666666");
1214  //very bad element // ser a very high center for it to be eliminated
1215  xc[0]=0.0; xc[1]=0.0; xc[2]=0.0;
1216  radius = 10000000;
1217 // KRATOS_WATCH("777777777777777777777777");
1218  }
1219  else
1220  {
1221 
1222 // KRATOS_WATCH("888888888888888888888888");
1223  double x0_2 = x0*x0 + y0*y0 + z0*z0;
1224  double x1_2 = x1*x1 + y1*y1 + z1*z1;
1225  double x2_2 = x2*x2 + y2*y2 + z2*z2;
1226  double x3_2 = x3*x3 + y3*y3 + z3*z3;
1227 
1228  //finalizing the calculation of the inverted matrix
1229  mJinv /= detJ;
1230 
1231  //calculating the RHS
1232  //calculating the RHS
1233  mRhs[0] = 0.5*(x1_2 - x0_2);
1234  mRhs[1] = 0.5*(x2_2 - x0_2);
1235  mRhs[2] = 0.5*(x3_2 - x0_2);
1236 
1237  //calculate position of the center
1238  noalias(xc) = prod(mJinv,mRhs);
1239 
1240  //calculate radius
1241  radius = pow(xc[0] - x0,2);
1242  radius += pow(xc[1] - y0,2);
1243  radius += pow(xc[2] - z0,2);
1244  radius = sqrt(radius);
1245 // KRATOS_WATCH("999999999999999999999999999");
1246  }
1247  KRATOS_CATCH("")
1248  }
1249 
1250  template< class T, std::size_t dim >
1251  class PointDistance{
1252  public:
1253  double operator()( T const& p1, T const& p2 ){
1254  double dist = 0.0;
1255  for( std::size_t i = 0 ; i < dim ; i++){
1256  double tmp = p1[i] - p2[i];
1257  dist += tmp*tmp;
1258  }
1259  return sqrt(dist);
1260  }
1261  };
1262 
1263  template< class T, std::size_t dim >
1264  class DistanceCalculator{
1265  public:
1266  double operator()( T const& p1, T const& p2 ){
1267  double dist = 0.0;
1268  for( std::size_t i = 0 ; i < dim ; i++){
1269  double tmp = p1[i] - p2[i];
1270  dist += tmp*tmp;
1271  }
1272  return dist;
1273  }
1274  };
1275 
1279 
1283 
1284 
1288 
1289 
1293 
1294 
1298 
1300  TetGenPfemModelerVms& operator=(TetGenPfemModelerVms const& rOther);
1301 
1302 
1304 
1305  }; // Class TetGenPfemModelerVms
1306 
1308 
1311 
1312 
1316 
1317 
1319  inline std::istream& operator >> (std::istream& rIStream,
1320  TetGenPfemModelerVms& rThis);
1321 
1323  inline std::ostream& operator << (std::ostream& rOStream,
1324  const TetGenPfemModelerVms& rThis)
1325  {
1326  rThis.PrintInfo(rOStream);
1327  rOStream << std::endl;
1328  rThis.PrintData(rOStream);
1329 
1330  return rOStream;
1331  }
1333 
1334 
1335 } // namespace Kratos.
1336 
1337 #endif // KRATOS_TETGEN_PFEM_MODELER_VMS_H_INCLUDED defined
1338 
1339 
Short class definition.
Definition: bucket.h:57
Base class for all Conditions.
Definition: condition.h:59
virtual Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const
It creates a new condition pointer.
Definition: condition.h:205
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
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
VariablesList & GetNodalSolutionStepVariablesList()
Definition: model_part.h:549
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
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
Short class definition.
Definition: tetgen_pfem_refine.h:67
Short class definition.
Definition: tetgen_pfem_refine_vms.h:68
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: tetgen_pfem_refine_vms.h:925
~TetGenPfemModelerVms() override=default
Destructor.
TetGenPfemModelerVms()
Default constructor.
Definition: tetgen_pfem_refine_vms.h:81
void ReGenerateMesh(ModelPart &ThisModelPart, Element const &rReferenceElement, Condition const &rReferenceBoundaryCondition, EntitiesEraseProcess< Node > &node_erase, bool rem_nodes=true, bool add_nodes=true, double alpha_param=1.4, double h_factor=0.5)
Definition: tetgen_pfem_refine_vms.h:103
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: tetgen_pfem_refine_vms.h:928
KRATOS_CLASS_POINTER_DEFINITION(TetGenPfemModelerVms)
Pointer definition of TetGenPfemModeler.
std::string Info() const override
Turn back information as a string.
Definition: tetgen_pfem_refine_vms.h:922
A four node tetrahedra geometry with linear shape functions.
Definition: tetrahedra_3d_4.h:59
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
A three node 3D triangle geometry with linear shape functions.
Definition: triangle_3d_3.h:77
bool Has(const VariableData &rThisVariable) const
Definition: variables_list.h:372
#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
z
Definition: GenerateWind.py:163
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
DistanceVector::iterator DistanceIterator
Definition: internal_variables_interpolation_process.h:56
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
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
Point PointType
Geometric definitions.
Definition: mortar_classes.h:36
float dist
Definition: edgebased_PureConvection.py:89
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
x2
Definition: generate_frictional_mortar_condition.py:122
x1
Definition: generate_frictional_mortar_condition.py:121
tuple tmp
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:98
out
Definition: isotropic_damage_automatic_differentiation.py:200
R
Definition: isotropic_damage_automatic_differentiation.py:172
float radius
Definition: mesh_to_mdpa_converter.py:18
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
int dim
Definition: sensitivityMatrix.py:25
volume
Definition: sp_statistics.py:15
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31