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.
particle_utilities.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Author Julio Marti.
11 //
12 
13 #if !defined(KRATOS_PARTICLES_UTILITIES_INCLUDED )
14 #define KRATOS_PARTICLES_UTILITIES_INCLUDED
15 
16 #define PRESSURE_ON_EULERIAN_MESH
17 #define USE_FEW_PARTICLES
18 
19 // System includes
20 #include <string>
21 #include <iostream>
22 #include <algorithm>
23 
24 // External includes
25 
26 
27 // Project includes
28 #include "includes/define.h"
29 #include "includes/model_part.h"
30 #include "includes/node.h"
31 #include "includes/kratos_flags.h"
32 #include "utilities/geometry_utilities.h"
33 #include "geometries/tetrahedra_3d_4.h"
34 #include "pfem_2_application_variables.h"
35 #include "spatial_containers/spatial_containers.h"
36 #include "utilities/timer.h"
37 #include "utilities/binbased_fast_point_locator.h"
38 //#include "utilities/enrichment_utilities.h"
39 
40 #include <boost/timer.hpp>
41 #include "utilities/timer.h"
42 
43 #ifdef _OPENMP
44 #include "omp.h"
45 #endif
46 
47 namespace Kratos
48 {
49 
50  template< class T, std::size_t dim >
51  class DistanceCalculator1
52  {
53  public:
54 
55  double operator()(T const& p1, T const& p2)
56  {
57  double dist = 0.0;
58  for (std::size_t i = 0; i < dim; i++)
59  {
60  double tmp = p1[i] - p2[i];
61  dist += tmp*tmp;
62  }
63  return dist; //square distance because it is easier to work without the square root//
64  }
65 
66  };
67 
68  template<std::size_t TDim> class ParticleUtils
69  {
70  public:
71  KRATOS_CLASS_POINTER_DEFINITION(ParticleUtils<TDim>);
72 
73 
74  void EstimateTime(ModelPart& rEulerianModelPart,const double max_dt)
75  {
76 
77  KRATOS_TRY
78  // KRATOS_ERROR(std::logic_error, "NEGATIVE VALUE OF Time step estimated" , "");
79  //initializee dt with max dt
80  //initialize dt with incredible value
81  double /*dt, glob_min_dt,*/ dummy;
82  // double h, nu;
83  array_1d<double,3> N = ZeroVector(3);
84  array_1d<double,3> aux = ZeroVector(3); //dimension = number of nodes
85  array_1d<double,3> vel = ZeroVector(3); //dimension = number of nodes
86  BoundedMatrix<double,3,2> DN_DX = ZeroMatrix(3,2);
87  array_1d<double,2> ms_vel_gauss = ZeroVector(2); //dimesion coincides with space dimension
88 
89  //initialize it with given value
90  // glob_min_dt=max_dt;
91 
92 
93  // dt=0.0;
94  for(ModelPart::ElementsContainerType::iterator im = rEulerianModelPart.ElementsBegin() ; im !=rEulerianModelPart.ElementsEnd() ; ++im)
95  {
96  GeometryUtils::CalculateGeometryData(im->GetGeometry(),DN_DX,N,dummy);
97 
98  double h = sqrt(2.00*dummy);
99 
100  array_1d<double,3> const& v = im->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY);
101  ms_vel_gauss[0] = v[0];
102  ms_vel_gauss[1] = v[1];
103 
104  //direction of the height is stored in the auxilliary vector
105  for (unsigned int i=1; i<3; i++)
106  {
107  array_1d<double,3> const& vi = im->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
108  ms_vel_gauss[0] += vi[0];
109  ms_vel_gauss[1] += vi[1];
110  }
111  ms_vel_gauss *=0.3333;
112 
113  double norm_u = ms_vel_gauss[0]*ms_vel_gauss[0] + ms_vel_gauss[1]*ms_vel_gauss[1];
114  norm_u = sqrt(norm_u);
115 
116  double courant= norm_u * max_dt / h;
117 
118  double& counter = im->GetValue(POISSON_RATIO);
119  counter = courant;
120 
121  }
122  KRATOS_CATCH("");
123  }
124 
125  void VisualizationModelPart(ModelPart& rCompleteModelPart, ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart)
126  {
127  KRATOS_TRY;
128 
129  rCompleteModelPart.Elements() = rEulerianModelPart.Elements();
130  rCompleteModelPart.Nodes() = rEulerianModelPart.Nodes();
131 
132  unsigned int id;
133  if(rEulerianModelPart.Nodes().size()!= 0)
134  id = (rEulerianModelPart.Nodes().end() - 1)->Id() + 1;
135  else
136  id = 1;
137 
138  //preallocate the memory needed
139  int tot_nodes = rEulerianModelPart.Nodes().size() + rLagrangianModelPart.Nodes().size();
140  rCompleteModelPart.Nodes().reserve( tot_nodes );
141 
142  //note that here we renumber the nodes
143  for (ModelPart::NodesContainerType::iterator node_it = rLagrangianModelPart.NodesBegin();
144  node_it != rLagrangianModelPart.NodesEnd(); node_it++)
145  {
146  node_it->SetId(id++);
147  rCompleteModelPart.AddNode(*(node_it.base()));
148  }
149 
150  KRATOS_CATCH("");
151  }
152 
153 
154  void TransferToEulerianMesh_Face_Heat_Flux(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart)
155  {
156  KRATOS_TRY
157  //defintions for spatial search
158  typedef Node PointType;
159  typedef Node ::Pointer PointTypePointer;
160  typedef std::vector<PointType::Pointer> PointVector;
161  typedef std::vector<PointType::Pointer>::iterator PointIterator;
162  typedef std::vector<double> DistanceVector;
163  typedef std::vector<double>::iterator DistanceIterator;
164 
165  //creating an auxiliary list for the new nodes
166  PointVector list_of_nodes;
167 
168  //*************
169  // Bucket types
170  typedef Bucket< TDim, PointType, PointVector, PointTypePointer, PointIterator, DistanceIterator > BucketType;
171 
172  typedef Tree< KDTreePartition<BucketType> > tree; //Kdtree;
173 
174  //starting calculating time of construction of the kdtree
175  boost::timer kdtree_construction;
176 
177  for (ModelPart::NodesContainerType::iterator node_it = rLagrangianModelPart.NodesBegin();
178  node_it != rLagrangianModelPart.NodesEnd(); ++node_it)
179  {
180  PointTypePointer pnode = *(node_it.base());
181 
182  //putting the nodes of the destination_model part in an auxiliary list
183  list_of_nodes.push_back(pnode);
184  }
185 
186  std::cout << "kdt constructin time " << kdtree_construction.elapsed() << std::endl;
187 
188  //create a spatial database with the list of new nodes
189  unsigned int bucket_size = 20;
190  tree nodes_tree(list_of_nodes.begin(), list_of_nodes.end(), bucket_size);
191 
192  //work arrays
193  Node work_point(0, 0.0, 0.0, 0.0);
194  unsigned int MaximumNumberOfResults = 10000;
195  PointVector Results(MaximumNumberOfResults);
196  DistanceVector SquaredResultsDistances(MaximumNumberOfResults);
197 
198 
199  if (rEulerianModelPart.NodesBegin()->SolutionStepsDataHas(NODAL_H) == false)
200  KRATOS_ERROR<<"Add ----NODAL_H---- variable!!!!!! ERROR";
201 
202  double sigma = 0.0;
203  if constexpr (TDim == 2)
204  sigma = 10.0 / (7.0 * 3.1415926);
205  else
206  sigma = 1.0 / 3.1415926;
207 
208  for (ModelPart::NodesContainerType::iterator node_it = rEulerianModelPart.NodesBegin(); node_it != rEulerianModelPart.NodesEnd(); node_it++)
209  {
210  if( (node_it)->FastGetSolutionStepValue(IS_FREE_SURFACE)==true || (node_it)->FastGetSolutionStepValue(IS_WATER)==1 )
211  { //IS_FREE_SURFACE
212  work_point.X() = node_it->X();
213  work_point.Y() = node_it->Y();
214  work_point.Z() = node_it->Z();
215 
216  double radius = 1.5 * node_it->FastGetSolutionStepValue(NODAL_H);
217 
218  //find all of the new nodes within the radius
219  int number_of_points_in_radius;
220 
221  //look between the new nodes which of them is inside the radius of the circumscribed cyrcle
222  number_of_points_in_radius = nodes_tree.SearchInRadius(work_point, radius, Results.begin(), SquaredResultsDistances.begin(), MaximumNumberOfResults);
223 
224  if (number_of_points_in_radius > 0)
225  {
226 
227  double& temperature = (node_it)->FastGetSolutionStepValue(FACE_HEAT_FLUX);
228  //double temperature=0.0;
229 
230  double temperature_aux = 0.0;
231 
232  double tot_weight = 0.0;
233 
234  for (int k = 0; k < number_of_points_in_radius; k++)
235  {
236  double distance = sqrt(*(SquaredResultsDistances.begin() + k));
237  double weight = SPHCubicKernel(sigma, distance, radius);
238 
239  PointIterator it_found = Results.begin() + k;
240 
241  if((*it_found)->FastGetSolutionStepValue(IS_INTERFACE)==1) //MATERIAL_VARIABLE
242  {
243  temperature_aux += weight * (*it_found)->FastGetSolutionStepValue(INCIDENT_RADIATION_FUNCTION);//);//FACE_HEAT_FLUX
244  tot_weight += weight;
245  }
246  }
247  if(tot_weight>0.0)
248  {
249  temperature_aux /= tot_weight;
250 
251  temperature +=(0.5 * temperature_aux * 1.00); //1.5 //1.25
252  }
253  }
254  }
255 
256  }
257  KRATOS_CATCH("")
258  }
259 
260  void TransferToEulerianMesh(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart)
261  {
262  KRATOS_TRY
263 
264  //defintions for spatial search
265  typedef Node PointType;
266  typedef Node ::Pointer PointTypePointer;
267  typedef std::vector<PointType::Pointer> PointVector;
268  typedef std::vector<PointType::Pointer>::iterator PointIterator;
269  typedef std::vector<double> DistanceVector;
270  typedef std::vector<double>::iterator DistanceIterator;
271 
272  //creating an auxiliary list for the new nodes
273  PointVector list_of_nodes;
274  //*************
275  // Bucket types
276  typedef Bucket< TDim, PointType, PointVector, PointTypePointer, PointIterator, DistanceIterator > BucketType;
277 
278  typedef Tree< KDTreePartition<BucketType> > tree; //Kdtree;
279 
280  //starting calculating time of construction of the kdtree
281  boost::timer kdtree_construction;
282 
283  for (ModelPart::NodesContainerType::iterator node_it = rLagrangianModelPart.NodesBegin();
284  node_it != rLagrangianModelPart.NodesEnd(); ++node_it)
285  {
286  PointTypePointer pnode = *(node_it.base());
287 
288  //putting the nodes of the destination_model part in an auxiliary list
289  list_of_nodes.push_back(pnode);
290  }
291 
292  std::cout << "kdt constructin time " << kdtree_construction.elapsed() << std::endl;
293 
294  //create a spatial database with the list of new nodes
295  unsigned int bucket_size = 20;
296  tree nodes_tree(list_of_nodes.begin(), list_of_nodes.end(), bucket_size);
297 
298  //work arrays
299  Node work_point(0, 0.0, 0.0, 0.0);
300  unsigned int MaximumNumberOfResults = 10000;
301  PointVector Results(MaximumNumberOfResults);
302  DistanceVector SquaredResultsDistances(MaximumNumberOfResults);
303 
304  if (rEulerianModelPart.NodesBegin()->SolutionStepsDataHas(NODAL_H) == false)
305  KRATOS_ERROR<<"Add ----NODAL_H---- variable!!!!!! ERROR";
306 
307  double sigma = 0.0;
308  if constexpr (TDim == 2)
309  sigma = 10.0 / (7.0 * 3.1415926);
310  else
311  sigma = 1.0 / 3.1415926;
312 
313  for (ModelPart::NodesContainerType::iterator node_it = rEulerianModelPart.NodesBegin(); node_it != rEulerianModelPart.NodesEnd(); node_it++)
314  {
315  if((node_it)->FastGetSolutionStepValue(IS_INTERFACE)==1)
316  { //IS_FREE_SURFACE
317  work_point.X() = node_it->X();
318  work_point.Y() = node_it->Y();
319  work_point.Z() = node_it->Z();
320  //KRATOS_ERROR(std::logic_error, "Add ----NODAL_H---- variable!!!!!! ERROR", "");
321  double radius = 2.0 * node_it->FastGetSolutionStepValue(NODAL_H);
322 
323  //find all of the new nodes within the radius
324  int number_of_points_in_radius;
325 
326  //look between the new nodes which of them is inside the radius of the circumscribed cyrcle
327  number_of_points_in_radius = nodes_tree.SearchInRadius(work_point, radius, Results.begin(), SquaredResultsDistances.begin(), MaximumNumberOfResults);
328 
329  if (number_of_points_in_radius > 0)
330  {
331  //double& temperature = (node_it)->FastGetSolutionStepValue(TEMPERATURE);
332  double temperature_aux = 0.0;
333  double tot_weight = 0.0;
334  for (int k = 0; k < number_of_points_in_radius; k++)
335  {
336  double distance = sqrt(*(SquaredResultsDistances.begin() + k));
337  double weight = SPHCubicKernel(sigma, distance, radius);
338  PointIterator it_found = Results.begin() + k;
339  //if((*it_found)->FastGetSolutionStepValue(IS_BOUNDARY)>0.5) //MATERIAL_VARIABLE
340  if((*it_found)->FastGetSolutionStepValue(IS_FREE_SURFACE) ==1 || (*it_found)->FastGetSolutionStepValue(IS_WATER) ==1 ) //MATERIAL_VARIABLE
341  {
342  double tempp=0.0;
343  tempp=(*it_found)->FastGetSolutionStepValue(YCH4);
344  //KRATOS_ERROR(std::logic_error, "nodo without temperature", "");
345  if(tempp<298.0) tempp=298.0;
346  //else tempp=(*it_found)->FastGetSolutionStepValue(YCH4);
347  temperature_aux += weight * tempp;//temperature
348  tot_weight += weight;
349  //KRATOS_ERROR(std::logic_error, "Add ----NODAL_H---- variable!!!!!! ERROR", "");
350  }
351  }
352  if(tot_weight>0.0)
353  {
354  temperature_aux /= tot_weight;
355  (node_it)->FastGetSolutionStepValue(FUEL)=temperature_aux;
356  }
357  else
358  {
359  KRATOS_WATCH(tot_weight);
360  KRATOS_WATCH((node_it)->X());
361  KRATOS_WATCH((node_it)->Y());
362  if((node_it)->FastGetSolutionStepValue(TEMPERATURE)<298.0) (node_it)->FastGetSolutionStepValue(FUEL)=298.0;
363  else (node_it)->FastGetSolutionStepValue(FUEL)=(node_it)->FastGetSolutionStepValue(TEMPERATURE);
364  }
365  }
366  }
367  else
368  {
369  //(node_it)->FastGetSolutionStepValue(FUEL)=(node_it)->FastGetSolutionStepValue(TEMPERATURE);
370  }
371  }
372  KRATOS_CATCH("")
373  }
374 
375  void TransferToEulerianMeshShapeBased_aux(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart, BinBasedFastPointLocator<TDim>& node_locator)
376  {
377  KRATOS_TRY
378 
379  Vector N;
380  const int max_results = 1000;
381  typename BinBasedFastPointLocator<TDim>::ResultContainerType results(max_results);
382  const int nparticles = rLagrangianModelPart.Nodes().size();
383 
384 #pragma omp parallel for firstprivate(results,N)
385  for (int i = 0; i < nparticles; i++)
386  {
387  ModelPart::NodesContainerType::iterator iparticle = rLagrangianModelPart.NodesBegin() + i;
388 
389  Node ::Pointer pparticle = *(iparticle.base());
390  typename BinBasedFastPointLocator<TDim>::ResultIteratorType result_begin = results.begin();
391 
392  Element::Pointer pelement;
393 
394  bool is_found = node_locator.FindPointOnMesh(pparticle->Coordinates(), N, pelement, result_begin, max_results);
395 
396  if (is_found == true)
397  {
398  Geometry<Node >& geom = pelement->GetGeometry();
399  BoundedMatrix<double, 3, 2 > msDN_DX;
400  array_1d<double, 3 > N;
401  //array_1d<double, 3 > N;
402  double Area=0.0;
403  GeometryUtils::CalculateGeometryData(geom, msDN_DX, N, Area);
404 
405  int s0=0;
406  int s1=0;
407  int s2=0;
408  int sum=0;
409  if(geom[0].FastGetSolutionStepValue(IS_INTERFACE)>0.5) s0=1; //IS_INTERFACE
410  if(geom[1].FastGetSolutionStepValue(IS_INTERFACE)>0.5) s1=1;
411  if(geom[2].FastGetSolutionStepValue(IS_INTERFACE)>0.5) s2=1;
412  sum=s0 + s1 + s2;
413  array_1d<double, 2 > qrad=ZeroVector(2);
414  array_1d<double, 2 > qrad_P1=ZeroVector(2);
415  array_1d<double,2> interface_segment=ZeroVector(2);
416  array_1d<double,2> normaledge1=ZeroVector(2);
417  for (unsigned int jj = 0; jj < 2; jj++)
418  {
419  for (unsigned int kk = 0; kk < 3; kk++)
420  {
421  qrad[jj] += msDN_DX(kk, jj) * geom[kk].FastGetSolutionStepValue(TEMPERATURE);
422  qrad_P1[jj] += msDN_DX(kk, jj) * geom[kk].FastGetSolutionStepValue(INCIDENT_RADIATION_FUNCTION);
423  }
424  }
425 
426  double faceheatflux=0.0;
427  //double faceheatflux_P1=0.0;
428  if(sum==2)
429  {
430  if((geom[1].FastGetSolutionStepValue(IS_INTERFACE)>0.5 && geom[0].FastGetSolutionStepValue(IS_INTERFACE)>0.5)) //IS_INTERFACE
431  {
432  double norm=0.0;
433  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
434  interface_segment[0] = (geom[0].X()-geom[1].X());
435  interface_segment[1] = (geom[0].Y()-geom[1].Y());
436  norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
437  //double area1=norm;
438  normaledge1(0)= -interface_segment[1]/norm;
439  normaledge1(1)= interface_segment[0]/norm;
440  faceheatflux += abs(1.0*(qrad[0]*normaledge1(0)+qrad[1]*normaledge1(1))*0.0131);
441  }
442 
443  if((geom[1].FastGetSolutionStepValue(IS_INTERFACE)>0.5 && geom[2].FastGetSolutionStepValue(IS_INTERFACE)>0.5))
444  {
445  double norm=0.0;
446  interface_segment[0] = (geom[1].X()-geom[2].X());
447  interface_segment[1] = (geom[1].Y()-geom[2].Y());
448  norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
449  //double area1=norm;
450  normaledge1(0)= -interface_segment[1]/norm;
451  normaledge1(1)= interface_segment[0]/norm;
452  faceheatflux += abs(1.0*(qrad[0]*normaledge1(0)+qrad[1]*normaledge1(1))*0.0131);
453  }
454  if((geom[2].FastGetSolutionStepValue(IS_INTERFACE)>0.5 && geom[0].FastGetSolutionStepValue(IS_INTERFACE)>0.5))
455  {
456  double norm=0.0;
457  interface_segment[0] = (geom[2].X()-geom[0].X());
458  interface_segment[1] = (geom[2].Y()-geom[0].Y());
459  norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
460  normaledge1(0)= -interface_segment[1]/norm;
461  normaledge1(1)= interface_segment[0]/norm;
462  faceheatflux += abs(1.0*(qrad[0]*normaledge1(0)+qrad[1]*normaledge1(1))*0.0131);
463  }
464  }
465  if(sum==1)
466  {
467  if((geom[1].FastGetSolutionStepValue(IS_INTERFACE)<0.5 && geom[0].FastGetSolutionStepValue(IS_INTERFACE)<0.5)) //IS_INTERFACE
468  {
469  double norm=0.0;
470  interface_segment[0] = (geom[0].X()-geom[1].X());
471  interface_segment[1] = (geom[0].Y()-geom[1].Y());
472  norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
473  normaledge1(0)= -interface_segment[1]/norm;
474  normaledge1(1)= interface_segment[0]/norm;
475  faceheatflux += abs(1.0*(qrad[0]*normaledge1(0)+qrad[1]*normaledge1(1))*0.0131);
476  }
477  if((geom[1].FastGetSolutionStepValue(IS_INTERFACE)<0.5 && geom[2].FastGetSolutionStepValue(IS_INTERFACE)<0.5))
478  {
479  double norm=0.0;
480  interface_segment[0] = (geom[1].X()-geom[2].X());
481  interface_segment[1] = (geom[1].Y()-geom[2].Y());
482  norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
483  normaledge1(0)= -interface_segment[1]/norm;
484  normaledge1(1)= interface_segment[0]/norm;
485  faceheatflux += abs(1.0*(qrad[0]*normaledge1(0)+qrad[1]*normaledge1(1))*0.0131);
486  }
487  if((geom[2].FastGetSolutionStepValue(IS_INTERFACE)<0.5 && geom[0].FastGetSolutionStepValue(IS_INTERFACE)<0.5))
488  {
489  double norm=0.0;
490  interface_segment[0] = (geom[2].X()-geom[0].X());
491  interface_segment[1] = (geom[2].Y()-geom[0].Y());
492  norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
493  normaledge1(0)= -interface_segment[1]/norm;
494  normaledge1(1)= interface_segment[0]/norm;
495  faceheatflux += abs(1.0*(qrad[0]*normaledge1(0)+qrad[1]*normaledge1(1))*0.0131);
496  }
497  }
498 
499  (iparticle)->FastGetSolutionStepValue(FACE_HEAT_FLUX)+=(faceheatflux /*+ fhf+ faceheatflux_P1*/);
500  }
501  }
502  KRATOS_CATCH("")
503  }
505  void CalculateNormal(ModelPart& full_model_part)
506  {
507  KRATOS_TRY
508  //resetting the normals
509  array_1d<double,3> zero;
510  noalias(zero) = ZeroVector(3);
511 
512  for(ModelPart::NodesContainerType::const_iterator in = full_model_part.NodesBegin(); in!=full_model_part.NodesEnd(); in++)
513  {
514  in->FastGetSolutionStepValue(NORMAL) = zero;
515  }
516 
517  array_1d<double,3> v1;
518  array_1d<double,3> v2;
519  //array_1d<double,3>& An =zero;
520  //double area_normal=0.0;
521  array_1d<double,3> area_normal;
522 
523  for(ModelPart::ElementsContainerType::iterator iii = full_model_part.ElementsBegin(); iii != full_model_part.ElementsEnd(); iii++)
524  {
525  if( iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[3].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0)
526  {
527  v1[0] = iii->GetGeometry()[1].X() -iii->GetGeometry()[3].X();
528  v1[1] = iii->GetGeometry()[1].Y() - iii->GetGeometry()[3].Y();
529  v1[2] = iii->GetGeometry()[1].Z() - iii->GetGeometry()[3].Z();
530 
531  v2[0] = iii->GetGeometry()[2].X() - iii->GetGeometry()[3].X();
532  v2[1] = iii->GetGeometry()[2].Y() - iii->GetGeometry()[3].Y();
533  v2[2] = iii->GetGeometry()[2].Z() - iii->GetGeometry()[3].Z();
534 
535  MathUtils<double>::CrossProduct(area_normal,v1,v2);
536  //area_normal *= -0.5;
537 
538  array_1d<double,3> msAuxVec = ZeroVector(3);
539  double c0 = abs(area_normal[0]);
540  double c1 = abs(area_normal[1]);
541  double c2 = abs(area_normal[2]);
542  msAuxVec[0]=c0;
543  msAuxVec[1]=c1;
544  msAuxVec[2]=c2;
545  // double norm_c =norm_2(msAuxVec);
546 
547  double norm_u = msAuxVec[0]*msAuxVec[0] + msAuxVec[1]*msAuxVec[1] + msAuxVec[2]*msAuxVec[2];
548  double norm_c =sqrt(norm_u);
549 
550  iii->GetGeometry()[1].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
551  iii->GetGeometry()[2].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
552  iii->GetGeometry()[3].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
553 
554  }
555  if( iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[3].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0)
556  {
557  v1[0] = iii->GetGeometry()[0].X() -iii->GetGeometry()[2].X();
558  v1[1] = iii->GetGeometry()[0].Y() - iii->GetGeometry()[2].Y();
559  v1[2] = iii->GetGeometry()[0].Z() - iii->GetGeometry()[2].Z();
560 
561  v2[0] = iii->GetGeometry()[3].X() - iii->GetGeometry()[2].X();
562  v2[1] = iii->GetGeometry()[3].Y() - iii->GetGeometry()[2].Y();
563  v2[2] = iii->GetGeometry()[3].Z() - iii->GetGeometry()[2].Z();
564  MathUtils<double>::CrossProduct(area_normal,v1,v2);
565  //area_normal *= -0.5;
566  array_1d<double,3> msAuxVec = ZeroVector(3);
567  double c0 = abs(area_normal[0]);
568  double c1 = abs(area_normal[1]);
569  double c2 = abs(area_normal[2]);
570  msAuxVec[0]=c0;
571  msAuxVec[1]=c1;
572  msAuxVec[2]=c2;
573  //double norm_c =norm_2(msAuxVec);
574 
575  double norm_u = msAuxVec[0]*msAuxVec[0] + msAuxVec[1]*msAuxVec[1] + msAuxVec[2]*msAuxVec[2];
576  double norm_c =sqrt(norm_u);
577 
578  iii->GetGeometry()[0].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
579  iii->GetGeometry()[3].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
580  iii->GetGeometry()[2].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
581  }
582 
583  if( iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[3].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0)
584  {
585  v1[0] = iii->GetGeometry()[0].X() -iii->GetGeometry()[3].X();
586  v1[1] = iii->GetGeometry()[0].Y() - iii->GetGeometry()[3].Y();
587  v1[2] = iii->GetGeometry()[0].Z() - iii->GetGeometry()[3].Z();
588 
589  v2[0] = iii->GetGeometry()[1].X() - iii->GetGeometry()[3].X();
590  v2[1] = iii->GetGeometry()[1].Y() - iii->GetGeometry()[3].Y();
591  v2[2] = iii->GetGeometry()[1].Z() - iii->GetGeometry()[3].Z();
592 
593  MathUtils<double>::CrossProduct(area_normal,v1,v2);
594  //area_normal *= -0.5;
595  array_1d<double,3> msAuxVec = ZeroVector(3);
596  double c0 = abs(area_normal[0]);
597  double c1 = abs(area_normal[1]);
598  double c2 = abs(area_normal[2]);
599  msAuxVec[0]=c0;
600  msAuxVec[1]=c1;
601  msAuxVec[2]=c2;
602 
603  double norm_u = msAuxVec[0]*msAuxVec[0] + msAuxVec[1]*msAuxVec[1] + msAuxVec[2]*msAuxVec[2];
604  double norm_c =sqrt(norm_u);
605 
606  iii->GetGeometry()[0].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
607  iii->GetGeometry()[1].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
608  iii->GetGeometry()[3].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
609  }
610  if( iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0 && iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) == 1.0)
611  {
612  v1[0] = iii->GetGeometry()[0].X() -iii->GetGeometry()[1].X();
613  v1[1] = iii->GetGeometry()[0].Y() - iii->GetGeometry()[1].Y();
614  v1[2] = iii->GetGeometry()[0].Z() - iii->GetGeometry()[1].Z();
615 
616  v2[0] = iii->GetGeometry()[2].X() - iii->GetGeometry()[1].X();
617  v2[1] = iii->GetGeometry()[2].Y() - iii->GetGeometry()[1].Y();
618  v2[2] = iii->GetGeometry()[2].Z() - iii->GetGeometry()[1].Z();
619 
620  MathUtils<double>::CrossProduct(area_normal,v1,v2);
621  //area_normal *= -0.5;
622  array_1d<double,3> msAuxVec = ZeroVector(3);
623  double c0 = abs(area_normal[0]);
624  double c1 = abs(area_normal[1]);
625  double c2 = abs(area_normal[2]);
626  msAuxVec[0]=c0;
627  msAuxVec[1]=c1;
628  msAuxVec[2]=c2;
629  // double norm_c =norm_2(msAuxVec);
630 
631  double norm_u = msAuxVec[0]*msAuxVec[0] + msAuxVec[1]*msAuxVec[1] + msAuxVec[2]*msAuxVec[2];
632  double norm_c =sqrt(norm_u);
633 
634  iii->GetGeometry()[0].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
635  iii->GetGeometry()[2].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
636  iii->GetGeometry()[1].FastGetSolutionStepValue(NORMAL) += area_normal/ norm_c;
637  }
638 
639  }
640  for(ModelPart::NodesContainerType::iterator iii = full_model_part.NodesBegin(); iii != full_model_part.NodesEnd(); iii++)
641  {
642 
643  if(iii->FastGetSolutionStepValue(IS_BOUNDARY)==1.0){
644  array_1d<double,3>& value_y1 = iii->FastGetSolutionStepValue(NORMAL);
645  double norm_y1 =norm_2(value_y1);
646  value_y1 /=(norm_y1 + 1e-9);
647  }
648  }
649 
650  KRATOS_CATCH("")
651  }
652 
653  void TransferToEulerianMeshShapeBased_aux_3D(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart, BinBasedFastPointLocator<TDim>& node_locator)
654  {
655  KRATOS_TRY
656  //typedef Node PointType;
657  //typedef Node ::Pointer PointTypePointer;
658  Vector N;
659  const int max_results = 1000;
660  typename BinBasedFastPointLocator<TDim>::ResultContainerType results(max_results);
661  const int nparticles = rLagrangianModelPart.Nodes().size();
662 #pragma omp parallel for firstprivate(results,N)
663  for (int i = 0; i < nparticles; i++)
664  {
665  ModelPart::NodesContainerType::iterator iparticle = rLagrangianModelPart.NodesBegin() + i;
666  Node ::Pointer pparticle = *(iparticle.base());
667  typename BinBasedFastPointLocator<TDim>::ResultIteratorType result_begin = results.begin();
668  Element::Pointer pelement;
669  bool is_found = node_locator.FindPointOnMesh(pparticle->Coordinates(), N, pelement, result_begin, max_results);
670  if (is_found == true)
671  {
672  Geometry<Node >& geom = pelement->GetGeometry();
673  BoundedMatrix<double, 4, 3 > msDN_DX;
674  array_1d<double, 4 > N;
675  double Area=0.0;
676  GeometryUtils::CalculateGeometryData(geom, msDN_DX, N, Area);
677  array_1d<double, 3 > qrad=ZeroVector(3);
678  double temmp=0.0;
679  for (unsigned int jj = 0; jj < 3; jj++)
680  {
681  for (unsigned int kk = 0; kk < 4; kk++)
682  {
683  temmp=geom[kk].FastGetSolutionStepValue(TEMPERATURE);
684  if(temmp<298.0) temmp=298.0;
685  qrad[jj] += msDN_DX(kk, jj) * temmp;//geom[kk].FastGetSolutionStepValue(TEMPERATURE);
686  }
687  }
688  //double faceheatflux=0.0;
689  (iparticle)->FastGetSolutionStepValue(NORMAL) *=(-1.0);
690  (iparticle)->FastGetSolutionStepValue(FACE_HEAT_FLUX) += abs( (iparticle)->FastGetSolutionStepValue(NORMAL_X) * qrad[0] + (iparticle)->FastGetSolutionStepValue(NORMAL_Y) * qrad[1] + (iparticle)->FastGetSolutionStepValue(NORMAL_Z) * qrad[2]) *0.0131;
691 
692  }
693  }
694  KRATOS_CATCH("")
695  }
696  //restarting the step from the beginning
697  void RestartStep(ModelPart & rModelPart)
698  {
699  KRATOS_TRY;
700 
701  //setting the variables to their value at the beginning of the time step
702  rModelPart.OverwriteSolutionStepData(1, 0);
703 
704  //setting the coordinates to their value at the beginning of the step
705  for (ModelPart::NodesContainerType::iterator node_it = rModelPart.NodesBegin();node_it != rModelPart.NodesEnd(); node_it++)
706  {
707  array_1d<double, 3 > & coords = node_it->Coordinates();
708  const array_1d<double, 3 > & old_disp = node_it->FastGetSolutionStepValue(DISPLACEMENT, 1);
709 
710  coords[0] = node_it->X0() + old_disp[0];
711  coords[1] = node_it->Y0() + old_disp[1];
712  coords[2] = node_it->Z0() + old_disp[2];
713  }
714 
715 
716  KRATOS_CATCH("");
717  }
718 
719 
720  void MoveMesh_Streamlines_freesurfaceflows(ModelPart& rModelPart, unsigned int substeps)
721  {
722  const double dt = rModelPart.GetProcessInfo()[DELTA_TIME];
723  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
724  BinBasedFastPointLocator<TDim> SearchStructure(rModelPart);
725  SearchStructure.UpdateSearchDatabase();
726 
727  //do movement
728  array_1d<double, 3 > veulerian;
729  //double temperature=0.0;
730  array_1d<double, 3 > acc_particle;
731  Vector N;
732  const int max_results = 10000;
733  typename BinBasedFastPointLocator<TDim>::ResultContainerType results(max_results);
734 
735  const int nparticles = rModelPart.Nodes().size();
736 
737 #pragma omp parallel for firstprivate(results,N,veulerian,acc_particle)
738  for (int i = 0; i < nparticles; i++)
739  {
740  //int substep = 0;
741  int subdivisions = 5;
742  //double temperature=0.0;
743  ModelPart::NodesContainerType::iterator iparticle = rModelPart.NodesBegin() + i;
744  Node ::Pointer pparticle = *(iparticle.base());
745  //small_dt = dt / subdivisions;
746 
747  bool do_move = true;
748  bool first_time=false;
749  iparticle->FastGetSolutionStepValue(DISTANCE)=0.0;
750 
751  iparticle->FastGetSolutionStepValue(EMBEDDED_VELOCITY) = iparticle->FastGetSolutionStepValue(VELOCITY,1); //AUX_VEL
752 
753  if(iparticle->Is(SLIP)) do_move = false;
754 
755  //iparticle->FastGetSolutionStepValue(TEMPERATURE) = 298.0;
756  if( do_move == true ) //note that we suppose the velocity components to be all fixed
757  {
758  array_1d<double,3> old_position = pparticle->Coordinates();
759  array_1d<double,3> current_position = pparticle->Coordinates();
760  noalias(iparticle->GetInitialPosition()) = old_position;
761  iparticle->FastGetSolutionStepValue(DISPLACEMENT,1) = ZeroVector(3);
762  //array_1d<double, 3 > & vel_particle = iparticle->FastGetSolutionStepValue(VELOCITY);
763  //subdivisions=10;
764  const double small_dt = dt / subdivisions;
765  //
766  for (int substep = 0; substep < subdivisions; substep++)
767  {
768  typename BinBasedFastPointLocator<TDim>::ResultIteratorType result_begin = results.begin();
769  Element::Pointer pelement;
770  bool is_found = SearchStructure.FindPointOnMesh(current_position, N, pelement, result_begin, max_results);
771  iparticle->Set(TO_ERASE, true);
772  //(iparticle)->GetValue(ERASE_FLAG) = true;
773  //KRATOS_WATCH(is_found);
774  if (is_found == true)
775  {
776  Geometry< Node >& geom = pelement->GetGeometry();
777  //int nn=0;
778  noalias(veulerian) = ZeroVector(3); //0.0;//N[0] * geom[0].FastGetSolutionStepValue(VELOCITY,1);
779  //temperature=0.0;//N[0] * geom[0].FastGetSolutionStepValue(TEMPERATURE);
780  for (unsigned int k = 0; k < geom.size(); k++)
781  {
782  noalias(veulerian) += N[k] * geom[k].FastGetSolutionStepValue(VELOCITY,1);
783  }
784  /*if(iparticle->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET)==1)
785  {
786  veulerian(0)*=0.0;
787  veulerian(2)*=0.0;
788  }*/
789  first_time=true;
790  noalias(current_position) += small_dt*veulerian;
791  pparticle->Set(TO_ERASE, false);
792  iparticle->FastGetSolutionStepValue(DISTANCE) += small_dt;
793  iparticle->FastGetSolutionStepValue(EMBEDDED_VELOCITY)=veulerian;
794  }
795  else
796  {
797  double time1=iparticle->FastGetSolutionStepValue(DISTANCE);
798  array_1d<double,3> acc;
799  acc[0] = 0.0;
800  acc[1] = -10.0;
801  acc[2] = 0.0;
802  if( first_time == false /*&& iparticle->Is(SLIP) == false*/ )
803  {
804  noalias(current_position) += small_dt *iparticle->FastGetSolutionStepValue(EMBEDDED_VELOCITY);
805  //noalias(current_position) += small_dt * small_dt * acc;
806  pparticle->Set(TO_ERASE, false);
807  }
808  else
809  {
810  time1 -=small_dt;
811  //double tiempo_restante=dt-time1;
812  noalias(current_position) += small_dt *iparticle->FastGetSolutionStepValue(EMBEDDED_VELOCITY);
813  //noalias(current_position) += small_dt * small_dt * acc;
814  pparticle->Set(TO_ERASE, false);
815  }
816 
817  }
818  }//for
819 
820  //update the displacement BUT DO NOT OVERWRITE THE POSITION!!
821  iparticle->FastGetSolutionStepValue(DISPLACEMENT) = current_position - iparticle->GetInitialPosition();
822  //KRATOS_WATCH(iparticle->FastGetSolutionStepValue(DISPLACEMENT));
823  }//move
824  }
825  //compute mesh velocity
826  for(ModelPart::NodesContainerType::iterator it = rModelPart.NodesBegin(); it!=rModelPart.NodesEnd(); it++)
827  {
828  //array_1d<double,3>& dn = it->FastGetSolutionStepValue(DISPLACEMENT,1);
829  array_1d<double,3>& dn1 = it->FastGetSolutionStepValue(DISPLACEMENT);
830  noalias(it->Coordinates()) = it->GetInitialPosition();
831  noalias(it->Coordinates()) += dn1;
832  }
833  }
834 
835  void MoveLonelyNodes(ModelPart& ThisModelPart)
836  {
837  KRATOS_TRY;
838  double Dt = ThisModelPart.GetProcessInfo()[DELTA_TIME];
839  array_1d<double,3> DeltaDisp, acc;
840 
841  for(ModelPart::NodeIterator i = ThisModelPart.NodesBegin() ;
842  i != ThisModelPart.NodesEnd() ; ++i)
843  {
844  if(
845  (i)->Is(SLIP) == false &&
846  (i)->GetValue(NEIGHBOUR_ELEMENTS).size() == 0 &&
847  ((i)->GetDof(VELOCITY_X).IsFixed() == false || (i)->GetDof(VELOCITY_Y).IsFixed() == false || (i)->GetDof(VELOCITY_Z).IsFixed() == false)
848  )
849  {
850  //i->Set(TO_ERASE,true);
851  //set to zero the pressure
852  (i)->FastGetSolutionStepValue(PRESSURE) = 0;
853  const array_1d<double,3>& old_vel = (i)->FastGetSolutionStepValue(VELOCITY,1);
854  array_1d<double,3>& vel = (i)->FastGetSolutionStepValue(VELOCITY);
855  //array_1d<double,3>& acc = (i)->FastGetSolutionStepValue(ACCELERATION);
856  noalias(acc) = (i)->FastGetSolutionStepValue(BODY_FORCE);
857  acc[0]= 0.0;
858  acc[1]= -10.0;
859  acc[2]= 0.0;
860  noalias(vel) = old_vel;
861  noalias(vel) += Dt * acc ;
862 
863  //calculate displacements
864  //noalias(DeltaDisp) = Dt * vel;
865 
866  //array_1d<double,3>& disp = i->FastGetSolutionStepValue(DISPLACEMENT);
867  //noalias(disp) = i->FastGetSolutionStepValue(DISPLACEMENT,1);
868  //noalias(disp) += DeltaDisp;
869  noalias(i->Coordinates()) += Dt * Dt * acc;
870 
871  }
872 
873  }
874 
875  KRATOS_CATCH("")
876  }
877 
878 
879  void MarkExcessivelyCloseNodes(ModelPart::NodesContainerType& rNodes)
880  {
881  KRATOS_TRY;
882  KRATOS_WATCH("ENTERD Mark close nodes")
883  //double fact2 = admissible_distance_factor*admissible_distance_factor;
884 
885  for(ModelPart::NodesContainerType::iterator in = rNodes.begin(); in!=rNodes.end(); in++)
886  {
887  if(in->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) ==1) //if it is not a wall node i can erase
888  {
889 
890  int nf=0;
891  //loop on neighbours and erase if they are too close
892  for( GlobalPointersVector< Node >::iterator i = in->GetValue(NEIGHBOUR_NODES).begin(); i != in->GetValue(NEIGHBOUR_NODES).end(); i++)
893  {
894 
895  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
896  if( /*i->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) ==1 and*/ i->FastGetSolutionStepValue(IS_FREE_SURFACE) ==1) //we can erase the current node only if the neighb is not to be erased
897  {
898  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
899  nf++;
900  //KRATOS_WATCH(nf)
901  }
902  if(nf>=2) {in->FastGetSolutionStepValue(IS_WATER)= 1;
903 
904  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
905  }
906  }
907  }
908  }
909 
910  KRATOS_CATCH("")
911  }
912 
913  void TransferToParticlesAirVelocity(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart, BinBasedFastPointLocator<TDim>& node_locator)
914  {
915  KRATOS_TRY
916 
917  //defintions for spatial search
918  //typedef Node PointType;
919  //typedef Node ::Pointer PointTypePointer;
920 
921 
922  Vector N;
923  const int max_results = 1000;
924  typename BinBasedFastPointLocator<TDim>::ResultContainerType results(max_results);
925  const int nparticles = rLagrangianModelPart.Nodes().size();
926 
927 #pragma omp parallel for firstprivate(results,N)
928  for (int i = 0; i < nparticles; i++)
929  {
930  ModelPart::NodesContainerType::iterator iparticle = rLagrangianModelPart.NodesBegin() + i;
931 
932  Node ::Pointer pparticle = *(iparticle.base());
933  typename BinBasedFastPointLocator<TDim>::ResultIteratorType result_begin = results.begin();
934 
935  Element::Pointer pelement;
936 
937  bool is_found = node_locator.FindPointOnMesh(pparticle->Coordinates(), N, pelement, result_begin, max_results);
938 
939  if (is_found == true)
940  {
941  Geometry<Node >& geom = pelement->GetGeometry();
942 
943  BoundedMatrix<double, 4, 3 > msDN_DX;
944 
945  array_1d<double, 4 > N;
946 
947  double Area=0.0;
948  GeometryUtils::CalculateGeometryData(geom, msDN_DX, N, Area);
949 
950  array_1d<double, 3 > velocity=ZeroVector(3);
951 
952  array_1d<double, 3 > temmp=ZeroVector(3);
953  //double temmp=0.0;
954  for (unsigned int jj = 0; jj < 3; jj++)
955  {
956  temmp=geom[jj].FastGetSolutionStepValue(VELOCITY);
957  velocity =N(jj) * temmp;
958  }
959  //KRATOS_WATCH(qrad);
960  //double faceheatflux=0.0;
961  (iparticle)->FastGetSolutionStepValue(ANGULAR_VELOCITY) = velocity;
962  }
963  }
964  KRATOS_CATCH("")
965  }
966 
967  double Calculate_Vol(ModelPart & rLagrangianModelPart)
968  {
969  KRATOS_TRY
970 
971  //defintions for spatial search
972  //typedef Node PointType;
973  //typedef Node ::Pointer PointTypePointer;
974 
975  //particles
976  for (ModelPart::NodesContainerType::iterator node_it = rLagrangianModelPart.NodesBegin(); node_it != rLagrangianModelPart.NodesEnd(); node_it++)
977  {
978  if( node_it->GetValue(NEIGHBOUR_ELEMENTS).size() != 0) (node_it)->FastGetSolutionStepValue(K0) = 0.0;
979  //if( node_it->FastGetSolutionStepValue(NODAL_MASS) == 0.0) KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
980  }
981 
982  for (ModelPart::ElementsContainerType::iterator el_it = rLagrangianModelPart.ElementsBegin();el_it != rLagrangianModelPart.ElementsEnd(); el_it++)
983  {
984 
985  Geometry<Node >& geom = el_it->GetGeometry();
986  double x0 = geom[0].X();
987  double y0 = geom[0].Y();
988  double z0 = geom[0].Z();
989  double x1 = geom[1].X();
990  double y1 = geom[1].Y();
991  double z1 = geom[1].Z();
992  double x2 = geom[2].X();
993  double y2 = geom[2].Y();
994  double z2 = geom[2].Z();
995  double x3 = geom[3].X();
996  double y3 = geom[3].Y();
997  double z3 = geom[3].Z();
998  double area=0.0;
999 
1000  area=CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3);
1001 
1002  geom[0].FastGetSolutionStepValue(K0) += area * 0.25;
1003  geom[1].FastGetSolutionStepValue(K0) += area * 0.25;
1004  geom[2].FastGetSolutionStepValue(K0) += area * 0.25;
1005  geom[3].FastGetSolutionStepValue(K0) += area * 0.25;
1006 
1007  }
1008 
1009  double sum=0.0;
1010 
1011  for (ModelPart::NodesContainerType::iterator node_it = rLagrangianModelPart.NodesBegin(); node_it != rLagrangianModelPart.NodesEnd(); node_it++)
1012  {
1013 
1014  sum +=(node_it)->FastGetSolutionStepValue(K0) ;
1015  }
1016 
1017  return sum;
1018 
1019  KRATOS_CATCH("")
1020  }
1021 
1022  void DetectAllOilClusters(ModelPart & mp_local_model_part)
1023  {
1024  int mnumber_of_oil_clusters=0;
1025  for (ModelPart::NodesContainerType::iterator inode = mp_local_model_part.NodesBegin(); inode != mp_local_model_part.NodesEnd(); inode++)
1026  {
1027  inode->FastGetSolutionStepValue(DIAMETER) = -1; //OIL_CLUSTER
1028  }
1029 
1030  for (ModelPart::ElementsContainerType::iterator ielem = mp_local_model_part.ElementsBegin();ielem != mp_local_model_part.ElementsEnd(); ielem++)
1031  {
1032  Geometry< Node >& geom = ielem->GetGeometry();
1033  if(geom.size()>1)
1034  {
1035  ielem->GetValue(DIAMETER) = -1;
1036  }
1037  }
1038 
1039 
1040  //fist we paint all the nodes connected to the outlet:
1041  int color = 0;
1042  for (ModelPart::NodesContainerType::iterator inode = mp_local_model_part.NodesBegin(); inode != mp_local_model_part.NodesEnd(); inode++)
1043  {
1044  if(inode->IsFixed(POROSITY) && inode->FastGetSolutionStepValue(DIAMETER)!=0) //nodes connected to the outlet are flagged as cluster zero. // if(inode->IsFixed(CONNECTED_TO_OUTLET) && inode->FastGetSolutionStepValue(OIL_CLUSTER)!=0)
1045  {
1046  ColorOilClusters(inode, 0);
1047  }
1048  }
1049 
1050  //having painted those nodes, we proceed with the rest of the colours
1051  for (ModelPart::NodesContainerType::iterator inode = mp_local_model_part.NodesBegin(); inode != mp_local_model_part.NodesEnd(); inode++)
1052  {
1053  if(inode->FastGetSolutionStepValue(DIAMETER) < 0 )
1054  {
1055  color++;
1056  ColorOilClusters(inode, color);
1057  }
1058  }
1059 
1060  for (ModelPart::ElementsContainerType::iterator ielem = mp_local_model_part.ElementsBegin(); ielem != mp_local_model_part.ElementsEnd(); ielem++)
1061  {
1062  Geometry< Node >& geom = ielem->GetGeometry();
1063  if(geom.size()>1 && ielem->GetValue(DIAMETER) < 0 )
1064  {
1065  color++;
1066  ielem->GetValue(DIAMETER) = color;
1067  }
1068  }
1069 
1070  //finally we flag the nodes with cluster=0 as connected to outlet
1071  for (ModelPart::NodesContainerType::iterator inode = mp_local_model_part.NodesBegin(); inode != mp_local_model_part.NodesEnd(); inode++)
1072  {
1073  if(inode->FastGetSolutionStepValue(DIAMETER) == 0)
1074  inode->FastGetSolutionStepValue(POROSITY)=1.0;
1075  }
1076 
1077  mnumber_of_oil_clusters = color;
1078 
1079  double area=0.0;
1080  array_1d<double, 3 > velocity_a=ZeroVector(3);
1081  array_1d<double, 3 > velocity_p=ZeroVector(3);
1082  array_1d<double, 3 > temmp=ZeroVector(3);
1083  array_1d<double, 3 > drag_coefficient=ZeroVector(3);
1084 
1085  //KRATOS_WATCH(mnumber_of_oil_clusters);
1086  int zz= mnumber_of_oil_clusters + 1;
1087  for(int jj=0; jj< zz; jj++ )
1088  {
1089  if(jj!=0)
1090  {
1091  if(jj==0) KRATOS_ERROR<<"element with zero vol found";
1092  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
1093  area=0.0;
1094  velocity_a=ZeroVector(3);
1095  velocity_p=ZeroVector(3);
1096  drag_coefficient=ZeroVector(3);
1097  int nn=0;
1098  for (ModelPart::NodesContainerType::iterator inode = mp_local_model_part.NodesBegin(); inode != mp_local_model_part.NodesEnd(); inode++)
1099  {
1100  int colour_p = (inode)->FastGetSolutionStepValue(DIAMETER);
1101  if(colour_p==jj)
1102  {
1103  area += (inode)->FastGetSolutionStepValue(K0);
1104  velocity_a += (inode)->FastGetSolutionStepValue(ANGULAR_VELOCITY);
1105  velocity_p += (inode)->FastGetSolutionStepValue(VELOCITY);
1106  nn++;
1107  }
1108  }
1109  velocity_a *=(1.0/nn);
1110  velocity_p *=(1.0/nn);
1111  //KRATOS_WATCH("AREA_ANTES");
1112  //KRATOS_WATCH("area");
1113  ComputedDragCoefficient(area, velocity_a, velocity_p, drag_coefficient );
1114 
1115  for (ModelPart::NodesContainerType::iterator inode = mp_local_model_part.NodesBegin(); inode != mp_local_model_part.NodesEnd(); inode++)
1116  {
1117  int colour_p = (inode)->FastGetSolutionStepValue(DIAMETER);
1118  if(colour_p==jj)
1119  {
1120  inode->FastGetSolutionStepValue(DRAG_FORCE_X)=drag_coefficient(0);
1121  inode->FastGetSolutionStepValue(DRAG_FORCE_Y)=drag_coefficient(1);
1122  inode->FastGetSolutionStepValue(DRAG_FORCE_Z)=drag_coefficient(2);
1123  }
1124  }
1125 
1126  }
1127 
1128  }
1129 
1130  }
1131 
1132  void ComputedDragCoefficient(double nodal_mass, array_1d<double, 3> velocity_air, array_1d<double, 3> velocity_polymer, array_1d<double, 3> & drag_coefficient )
1133  {
1134  KRATOS_TRY
1135 
1136  double drag_coeff=0.0;
1137  //array_1d<double, 3> drag_coefficient=ZeroVector(3);
1138  //nodal_mass=0.0;
1139  array_1d<double, 3> vrelative;
1140 
1141  //nodal_mass=(node_it)->FastGetSolutionStepValue(NODAL_MASS);
1142  double aux=nodal_mass * 3.0/(3.0*3.1416);
1143  double Radius= pow(aux, 0.3333333);
1144  double area=4.0 * 3.1416 * Radius * Radius;
1145  noalias(vrelative)=velocity_air-velocity_polymer;
1146  double norm_u = norm_2(vrelative);
1147  double reynolds = 2 * Radius * norm_u / 0.00001; // 2 * mRadius * mNormOfSlipVel / mKinematicViscosity
1148  if (reynolds < 0.01)
1149  {
1150  reynolds = 0.01;
1151  }
1152  CalculateNewtonianDragCoefficient(reynolds, drag_coeff);
1153  noalias(drag_coefficient) = 0.5 * 1.0 * area * drag_coeff * norm_u* vrelative * (1.0 / nodal_mass); //drag_coeff = 0.5 * mFluidDensity * area * drag_coeff * mNormOfSlipVel;
1154  KRATOS_CATCH("")
1155  }
1156 
1157 
1158 
1159  void CalculateNewtonianDragCoefficient(const double reynolds, double& drag_coeff)
1160  {
1161  KRATOS_TRY
1162 
1163  if (reynolds < 1){
1164  drag_coeff = 24.0; // Reynolds;
1165  }
1166  else {
1167  if (reynolds > 1000){
1168  drag_coeff = 0.44;
1169  }
1170  else{
1171  drag_coeff = 24.0 / reynolds * (1.0 + 0.15 * pow(reynolds, 0.687));
1172  }
1173  }
1174 
1175  KRATOS_CATCH("")
1176  }
1177 
1178 
1179  void ColorOilClusters(ModelPart::NodesContainerType::iterator iNode, const int color)
1180  {
1181  if(iNode->GetSolutionStepValue(DIAMETER) < 0 ) // if(iNode->GetSolutionStepValue(OIL_CLUSTER) < 0 && ( water_fraction<0.99999999999999 || theta>1.57079632679) )
1182  iNode->GetSolutionStepValue(DIAMETER)=color;
1183 
1184  ModelPart::NodesContainerType front_nodes;
1185  GlobalPointersVector<Element >& r_neighbour_elements = iNode->GetValue(NEIGHBOUR_ELEMENTS);
1186  for(GlobalPointersVector<Element >::iterator i_neighbour_element = r_neighbour_elements.begin() ; i_neighbour_element != r_neighbour_elements.end() ; i_neighbour_element++)
1187  {
1188  if(i_neighbour_element->GetValue(DIAMETER) < 0 )
1189  {
1190  i_neighbour_element->SetValue(DIAMETER, color);
1191 
1192  Element::GeometryType& p_geometry = i_neighbour_element->GetGeometry();
1193 
1194  for(unsigned int i = 0; i < p_geometry.size(); i++)
1195  {
1196  if(p_geometry[i].GetSolutionStepValue(DIAMETER) < 0 )
1197  {
1198  p_geometry[i].GetSolutionStepValue(DIAMETER) = color;
1199  front_nodes.push_back(p_geometry(i));
1200  }
1201  }
1202  }
1203  }
1204  while(!front_nodes.empty())
1205  {
1206  ModelPart::NodesContainerType new_front_nodes;
1207  for(ModelPart::NodesContainerType::iterator i_node = front_nodes.begin() ; i_node != front_nodes.end() ; i_node++)
1208  {
1209  GlobalPointersVector<Element >& r_neighbour_elements = i_node->GetValue(NEIGHBOUR_ELEMENTS);
1210  for(GlobalPointersVector<Element >::iterator i_neighbour_element = r_neighbour_elements.begin() ; i_neighbour_element != r_neighbour_elements.end() ; i_neighbour_element++)
1211  {
1212  if(i_neighbour_element->GetValue(DIAMETER) < 0 )
1213  {
1214  i_neighbour_element->SetValue(DIAMETER, color);
1215 
1216  Element::GeometryType& p_geometry = i_neighbour_element->GetGeometry();
1217 
1218  for(unsigned int i = 0; i < p_geometry.size(); i++)
1219  {
1220  if(p_geometry[i].GetSolutionStepValue(DIAMETER) < 0 )
1221  {
1222  p_geometry[i].GetSolutionStepValue(DIAMETER) = color;
1223  new_front_nodes.push_back(p_geometry(i));
1224  }
1225  }
1226  }
1227  }
1228  }
1229  front_nodes.clear();// (o resize ( 0 ), si clear no existe)
1230  for( ModelPart::NodesContainerType::iterator i_node = new_front_nodes.begin() ; i_node != new_front_nodes.end() ; i_node++)
1231  front_nodes.push_back(*(i_node.base()));
1232 
1233  }
1234  }
1235 
1236 
1237  void movethermocouples(ModelPart& rEulerianModelPart, ModelPart& rLagrangianModelPart, BinBasedFastPointLocator<TDim>& node_locator)
1238  {
1239  KRATOS_TRY
1240  array_1d<double, 3 > veulerian;
1241  double temperature;
1242  Vector N;
1243  //double G;
1244  const int max_results = 1000;
1245  typename BinBasedFastPointLocator<TDim>::ResultContainerType results(max_results);
1246  double dt =0.01;
1247  const int nparticles = rLagrangianModelPart.Nodes().size();
1248 
1249 #pragma omp parallel for firstprivate(results,N,veulerian,temperature)
1250  for (int i = 0; i < nparticles; i++)
1251  {
1252  ModelPart::NodesContainerType::iterator iparticle = rLagrangianModelPart.NodesBegin() + i;
1253 
1254  int subdivisions=5.0;
1255  const double small_dt = dt / subdivisions;
1256  for (unsigned int substep = 0; substep < subdivisions; substep++)
1257  {
1258  Node ::Pointer pparticle = *(iparticle.base());
1259  typename BinBasedFastPointLocator<TDim>::ResultIteratorType result_begin = results.begin();
1260  Element::Pointer pelement;
1261  bool is_found = node_locator.FindPointOnMesh(pparticle->Coordinates(), N, pelement, result_begin, max_results);
1262  if (is_found == true)
1263  {
1264  Geometry< Node >& geom = pelement->GetGeometry();
1265 
1266  //move according to the streamline
1267  noalias(veulerian) = N[0] * geom[0].FastGetSolutionStepValue(VELOCITY, 1);
1268  temperature = N[0] * geom[0].FastGetSolutionStepValue(YCH4);
1269  for (unsigned int k = 1; k < geom.size(); k++)
1270  {
1271  noalias(veulerian) += N[k] * geom[k].FastGetSolutionStepValue(VELOCITY, 1);
1272  temperature += N[k] * geom[k].FastGetSolutionStepValue(YCH4);
1273  //KRATOS_WATCH(geom[k].FastGetSolutionStepValue(YCH4));
1274  }
1275  double & temp = (iparticle)->FastGetSolutionStepValue(YCH4);
1276  temp =temperature;
1277  veulerian(0) *=0.0;
1278  veulerian(2) *=0.0;
1279  array_1d<double, 3 > & disp = (iparticle)->FastGetSolutionStepValue(DISPLACEMENT);
1280  noalias(disp) += small_dt*veulerian;
1281  noalias(iparticle->Coordinates()) = iparticle->GetInitialPosition();
1282  noalias(iparticle->Coordinates()) += iparticle->FastGetSolutionStepValue(DISPLACEMENT);
1283  }
1284  }
1285  }
1286 
1287  KRATOS_CATCH("")
1288  }
1289 
1290  void TransferToEulerianMesh_2(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart)
1291  {
1292  KRATOS_TRY
1293 
1294  //defintions for spatial search
1295  typedef Node PointType;
1296  typedef Node ::Pointer PointTypePointer;
1297  typedef std::vector<PointType::Pointer> PointVector;
1298  typedef std::vector<PointType::Pointer>::iterator PointIterator;
1299  typedef std::vector<double> DistanceVector;
1300  typedef std::vector<double>::iterator DistanceIterator;
1301 
1302  //creating an auxiliary list for the new nodes
1303  PointVector list_of_nodes;
1304 
1305  //*************
1306  // Bucket types
1307  typedef Bucket< TDim, PointType, PointVector, PointTypePointer, PointIterator, DistanceIterator > BucketType;
1308 
1309  typedef Tree< KDTreePartition<BucketType> > tree; //Kdtree;
1310 
1311 
1312  //starting calculating time of construction of the kdtree
1313  boost::timer kdtree_construction;
1314 
1315  for (ModelPart::NodesContainerType::iterator node_it = rLagrangianModelPart.NodesBegin();
1316  node_it != rLagrangianModelPart.NodesEnd(); ++node_it)
1317  {
1318  PointTypePointer pnode = *(node_it.base());
1319 
1320  //putting the nodes of the destination_model part in an auxiliary list
1321  list_of_nodes.push_back(pnode);
1322  }
1323 
1324  std::cout << "kdt constructin time " << kdtree_construction.elapsed() << std::endl;
1325 
1326  //create a spatial database with the list of new nodes
1327  unsigned int bucket_size = 20;
1328  tree nodes_tree(list_of_nodes.begin(), list_of_nodes.end(), bucket_size);
1329 
1330  //work arrays
1331  Node work_point(0, 0.0, 0.0, 0.0);
1332  unsigned int MaximumNumberOfResults = 10000;
1333  PointVector Results(MaximumNumberOfResults);
1334  DistanceVector SquaredResultsDistances(MaximumNumberOfResults);
1335 
1336 
1337  if (rEulerianModelPart.NodesBegin()->SolutionStepsDataHas(NODAL_H) == false)
1338  KRATOS_ERROR<<"Add ----NODAL_H---- variable!!!!!! ERROR";
1339 
1340  double sigma = 0.0;
1341 
1342  if constexpr (TDim == 2)
1343  sigma = 10.0 / (7.0 * 3.1415926);
1344  else
1345  sigma = 1.0 / 3.1415926;
1346 
1347  for (ModelPart::NodesContainerType::iterator node_it = rEulerianModelPart.NodesBegin(); node_it != rEulerianModelPart.NodesEnd(); node_it++)
1348  {
1349  if((node_it)->Y()< -0.15102 )
1350  {
1351  work_point.X() = node_it->X();
1352  work_point.Y() = node_it->Y();
1353  work_point.Z() = node_it->Z();
1354 
1355  double radius = 2.0 * node_it->FastGetSolutionStepValue(NODAL_H);
1356 
1357  //find all of the new nodes within the radius
1358  int number_of_points_in_radius;
1359 
1360  //look between the new nodes which of them is inside the radius of the circumscribed cyrcle
1361  number_of_points_in_radius = nodes_tree.SearchInRadius(work_point, radius, Results.begin(), SquaredResultsDistances.begin(), MaximumNumberOfResults);
1362 
1363  if (number_of_points_in_radius > 0)
1364  {
1365  //double& temperature = (node_it)->FastGetSolutionStepValue(TEMPERATURE);
1366  double temperature_aux = 0.0;
1367  double tot_weight = 0.0;
1368  double C = 1.19e15;
1369  double E_over_R = 24067.0;
1370 
1371  for (int k = 0; k < number_of_points_in_radius; k++)
1372  {
1373  double distance = sqrt(*(SquaredResultsDistances.begin() + k));
1374  double weight = SPHCubicKernel(sigma, distance, radius);
1375  PointIterator it_found = Results.begin() + k;
1376 
1377  if( (*it_found)->FastGetSolutionStepValue(DIAMETER) >0 && (*it_found)->FastGetSolutionStepValue(YN2) ==0.0) //MATERIAL_VARIABLE
1378  {
1379  double tempp=0.0;
1380  tempp=(*it_found)->FastGetSolutionStepValue(YCH4);
1381  //if(tempp<298.0) tempp=298.0;
1382  //else tempp=(*it_found)->FastGetSolutionStepValue(YCH4);
1383 
1384  temperature_aux += weight * 27400.0 * 1.0 * C * exp(-E_over_R / tempp ) * 905.0 ;//temperature
1385  tot_weight += weight;
1386 
1387 
1388  }
1389  }
1390  if(tot_weight>0.0)
1391  {
1392  //(node_it)->FastGetSolutionStepValue(FUEL)=1500.0;
1393  //double nodal_mass = node_it->FastGetSolutionStepValue(NODAL_MASS);
1394  //KRATOS_WATCH(node_it->FastGetSolutionStepValue(NODAL_MASS))
1395  //double& heat = node_it->FastGetSolutionStepValue(HEAT_FLUX);
1396  temperature_aux /= (tot_weight );
1397  //temperature= 1500.0;//temperature_aux;
1398  if(temperature_aux>1e9) (node_it)->FastGetSolutionStepValue(HEAT_FLUX) += 1e9;
1399  else (node_it)->FastGetSolutionStepValue(HEAT_FLUX) += temperature_aux;
1400  }
1401  }
1402  }
1403  }
1404  KRATOS_CATCH("")
1405  }
1406 
1407  void TransferToEulerianMeshShapeBased(ModelPart& rEulerianModelPart, ModelPart & rLagrangianModelPart, BinBasedFastPointLocator<TDim>& node_locator)
1408  {
1409  KRATOS_TRY
1410 
1411  for (ModelPart::NodesContainerType::iterator node_it = rEulerianModelPart.NodesBegin(); node_it != rEulerianModelPart.NodesEnd(); node_it++)
1412  {
1413  (node_it)->GetValue(POISSON_RATIO) = 0.0;
1414  (node_it)->GetValue(YOUNG_MODULUS) = 0.0;
1415  (node_it)->GetValue(NODAL_MASS) = 0.0;
1416  (node_it)->GetValue(HEAT_FLUX) = 0.0;
1417  }
1418 
1419  for (ModelPart::ElementsContainerType::iterator el_it = rEulerianModelPart.ElementsBegin();el_it != rEulerianModelPart.ElementsEnd(); el_it++)
1420  {
1421  el_it->SetValue(YOUNG_MODULUS,0.0);
1422  }
1423  Vector N;
1424  const int max_results = 1000;
1425  typename BinBasedFastPointLocator<TDim>::ResultContainerType results(max_results);
1426  const int nparticles = rLagrangianModelPart.Nodes().size();
1427 
1428  double C = 1.19e15;
1429  double E_over_R = 24067.0;
1430  double A=0.0;
1431 
1432 #pragma omp parallel for firstprivate(results,N)
1433  for (int i = 0; i < nparticles; i++)
1434  {
1435  ModelPart::NodesContainerType::iterator iparticle = rLagrangianModelPart.NodesBegin() + i;
1436  //KRATOS_ERROR(std::logic_error, "Add ----FORCE---- variable!!!!!! ERROR", "");
1437  Node ::Pointer pparticle = *(iparticle.base());
1438  typename BinBasedFastPointLocator<TDim>::ResultIteratorType result_begin = results.begin();
1439 
1440  Element::Pointer pelement;
1441 
1442  bool is_found = node_locator.FindPointOnMesh(pparticle->Coordinates(), N, pelement, result_begin, max_results);
1443 
1444  if (is_found == true)
1445  {
1446  Geometry<Node >& geom = pelement->GetGeometry();
1447  const array_1d<double, 3 > & vel_particle = (iparticle)->FastGetSolutionStepValue(VELOCITY);
1448  double density_particle = (iparticle)->FastGetSolutionStepValue(DENSITY);
1449  double Tp=(iparticle)->FastGetSolutionStepValue(YCH4); //HEAT_FLUX
1450  if(Tp>1000.0) Tp=1000.0;
1451  double temperature = 0.3e+8;//1000000.0;//(iparticle)->FastGetSolutionStepValue(TEMPERATURE);
1452 
1453  temperature = 2e+7;
1454 
1455  if( (iparticle)->FastGetSolutionStepValue(DIAMETER)>0)// ((iparticle)->GetValue(NEIGHBOUR_ELEMENTS)).size() == 0)
1456  {
1457  (iparticle)->FastGetSolutionStepValue(YN2)=1.0;
1458  // KRATOS_ERROR(std::logic_error, "Add ----FORCE---- variable!!!!!! ERROR", "");
1459  KRATOS_WATCH("aloneeeeeeeeeeeeeeeeeeeee");
1460  KRATOS_WATCH((iparticle)->FastGetSolutionStepValue(DIAMETER));
1461  for (unsigned int k = 0; k < geom.size(); k++)
1462  {
1463 
1464  //KRATOS_ERROR(std::logic_error, "Add ----FORCE---- variable!!!!!! ERROR", "");
1465  geom[k].SetLock();
1466  geom[k].GetValue(YOUNG_MODULUS) += N[k] * 27400.0 * 1.0 * C * exp(-E_over_R/(Tp)) *905.0 * (iparticle)->FastGetSolutionStepValue(K0);//0.0001;
1467 
1468  KRATOS_WATCH("K0");
1469  KRATOS_WATCH((iparticle)->FastGetSolutionStepValue(K0));
1470 
1471  geom[k].GetValue(POISSON_RATIO) += N[k];
1472 
1473  geom[k].UnSetLock();
1474  }
1475  }
1476  }
1477  }
1478 
1479  for (ModelPart::ElementsContainerType::iterator el_it = rEulerianModelPart.ElementsBegin();el_it != rEulerianModelPart.ElementsEnd(); el_it++)
1480  {
1481 
1482  Geometry<Node >& geom = el_it->GetGeometry();
1483  double x0 = geom[0].X();
1484  double y0 = geom[0].Y();
1485  double z0 = geom[0].Z();
1486  double x1 = geom[1].X();
1487  double y1 = geom[1].Y();
1488  double z1 = geom[1].Z();
1489  double x2 = geom[2].X();
1490  double y2 = geom[2].Y();
1491  double z2 = geom[2].Z();
1492  double x3 = geom[3].X();
1493  double y3 = geom[3].Y();
1494  double z3 = geom[3].Z();
1495  double area=0.0;
1496  //if constexpr (TDim==2) area=CalculateVol(x0, y0, x1, y1, x2, y2);
1497  //else
1498  area=CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3);
1499 
1500  geom[0].FastGetSolutionStepValue(NODAL_MASS) += area * 0.25;
1501  geom[1].FastGetSolutionStepValue(NODAL_MASS) += area * 0.25;
1502  geom[2].FastGetSolutionStepValue(NODAL_MASS) += area * 0.25;
1503  geom[3].FastGetSolutionStepValue(NODAL_MASS) += area * 0.25;
1504 
1505  }
1506 
1507 
1508  for (ModelPart::NodesContainerType::iterator node_it = rEulerianModelPart.NodesBegin(); node_it != rEulerianModelPart.NodesEnd(); node_it++)
1509  {
1510  const double NN = (node_it)->GetValue(POISSON_RATIO);
1511  const double tt = (node_it)->GetValue(YOUNG_MODULUS);
1512  if (NN != 0.0)
1513  {
1514  //KRATOS_ERROR(std::logic_error, "element with zero vol found", "");
1515  //KRATOS_WATCH(tt);
1516  //KRATOS_WATCH(NN);
1517  double nodal_mass = node_it->FastGetSolutionStepValue(NODAL_MASS);
1518  double& heat = node_it->FastGetSolutionStepValue(HEAT_FLUX);
1519  double heat_value=tt/NN * (1.0/nodal_mass);
1520  if (heat_value>1e9 /*0.5e+7*/) heat_value = 1e9 /*0.5e+7*/; //0.5e+7;
1521  heat= heat_value;
1522  KRATOS_WATCH(heat);
1523  }
1524  //}
1525  }
1526 
1527 
1528  // Timer::Stop("Interpolacion");
1529  //KRATOS_WATCH(time)
1530 
1531  KRATOS_CATCH("")
1532  }
1533 
1534  private:
1535 
1536  inline double SPHCubicKernel(const double sigma, const double r, const double hmax)
1537  {
1538  double h_half = 0.5 * hmax;
1539  const double s = r / h_half;
1540  const double coeff = sigma / pow(h_half, static_cast<int>(TDim));
1541 
1542  if (s <= 1.0)
1543  return coeff * (1.0 - 1.5 * s * s + 0.75 * s * s * s);
1544  else if (s <= 2.0)
1545  return 0.25 * coeff * pow(2.0 - s, 3);
1546  else
1547  return 0.0;
1548  }
1549 
1550  inline void CalculateCenterAndSearchRadius(Geometry<Node >&geom, double& xc, double& yc, double& zc, double& R, array_1d<double, 3 > & N )
1551  {
1552  double x0 = geom[0].X();
1553  double y0 = geom[0].Y();
1554  double x1 = geom[1].X();
1555  double y1 = geom[1].Y();
1556  double x2 = geom[2].X();
1557  double y2 = geom[2].Y();
1558 
1559  xc = 0.3333333333333333333 * (x0 + x1 + x2);
1560  yc = 0.3333333333333333333 * (y0 + y1 + y2);
1561  zc = 0.0;
1562 
1563  double R1 = (xc - x0)*(xc - x0) + (yc - y0)*(yc - y0);
1564  double R2 = (xc - x1)*(xc - x1) + (yc - y1)*(yc - y1);
1565  double R3 = (xc - x2)*(xc - x2) + (yc - y2)*(yc - y2);
1566 
1567  R = R1;
1568  if (R2 > R) R = R2;
1569  if (R3 > R) R = R3;
1570 
1571  R = 1.01 * sqrt(R);
1572  }
1573 
1574  inline void CalculateCenterAndSearchRadius(Geometry<Node >&geom, double& xc, double& yc, double& zc, double& R, array_1d<double, 4 > & N )
1575  {
1576  double x0 = geom[0].X();
1577  double y0 = geom[0].Y();
1578  double z0 = geom[0].Z();
1579  double x1 = geom[1].X();
1580  double y1 = geom[1].Y();
1581  double z1 = geom[1].Z();
1582  double x2 = geom[2].X();
1583  double y2 = geom[2].Y();
1584  double z2 = geom[2].Z();
1585  double x3 = geom[3].X();
1586  double y3 = geom[3].Y();
1587  double z3 = geom[3].Z();
1588 
1589 
1590  xc = 0.25 * (x0 + x1 + x2 + x3);
1591  yc = 0.25 * (y0 + y1 + y2 + y3);
1592  zc = 0.25 * (z0 + z1 + z2 + z3);
1593 
1594  double R1 = (xc - x0)*(xc - x0) + (yc - y0)*(yc - y0) + (zc - z0)*(zc - z0);
1595  double R2 = (xc - x1)*(xc - x1) + (yc - y1)*(yc - y1) + (zc - z1)*(zc - z1);
1596  double R3 = (xc - x2)*(xc - x2) + (yc - y2)*(yc - y2) + (zc - z2)*(zc - z2);
1597  double R4 = (xc - x3)*(xc - x3) + (yc - y3)*(yc - y3) + (zc - z3)*(zc - z3);
1598 
1599  R = R1;
1600  if (R2 > R) R = R2;
1601  if (R3 > R) R = R3;
1602  if (R4 > R) R = R4;
1603 
1604  R = sqrt(R);
1605  }
1606 
1607 
1608  inline bool CalculatePosition(Geometry<Node >&geom,const double xc, const double yc, const double zc, array_1d<double, 4 > & N )
1609  {
1610 
1611  double x0 = geom[0].X();
1612  double y0 = geom[0].Y();
1613  double z0 = geom[0].Z();
1614  double x1 = geom[1].X();
1615  double y1 = geom[1].Y();
1616  double z1 = geom[1].Z();
1617  double x2 = geom[2].X();
1618  double y2 = geom[2].Y();
1619  double z2 = geom[2].Z();
1620  double x3 = geom[3].X();
1621  double y3 = geom[3].Y();
1622  double z3 = geom[3].Z();
1623 
1624  double vol = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3);
1625 
1626  double inv_vol = 0.0;
1627  if (vol < 0.0000000000001)
1628  {
1629  KRATOS_ERROR<<"element with zero vol found";
1630  }
1631  else
1632  {
1633  inv_vol = 1.0 / vol;
1634  }
1635 
1636  N[0] = CalculateVol(x1, y1, z1, x3, y3, z3, x2, y2, z2, xc, yc, zc) * inv_vol;
1637  N[1] = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, xc, yc, zc) * inv_vol;
1638  N[2] = CalculateVol(x3, y3, z3, x1, y1, z1, x0, y0, z0, xc, yc, zc) * inv_vol;
1639  N[3] = CalculateVol(x3, y3, z3, x0, y0, z0, x2, y2, z2, xc, yc, zc) * inv_vol;
1640 
1641 
1642  if (N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[3] >= 0.0 &&
1643  N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0 && N[3] <= 1.0)
1644  //if the xc yc zc is inside the tetrahedron return true
1645  return true;
1646 
1647  return false;
1648  }
1649 
1650  inline double CalculateVol(const double x0, const double y0,
1651  const double x1, const double y1,
1652  const double x2, const double y2
1653  )
1654  {
1655  return 0.5 * ((x1 - x0)*(y2 - y0)- (y1 - y0)*(x2 - x0));
1656  }
1657 
1658 
1659  inline double CalculateVol(const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double x2, const double y2, const double z2,const double x3, const double y3, const double z3 )
1660  {
1661  double x10 = x1 - x0;
1662  double y10 = y1 - y0;
1663  double z10 = z1 - z0;
1664 
1665  double x20 = x2 - x0;
1666  double y20 = y2 - y0;
1667  double z20 = z2 - z0;
1668 
1669  double x30 = x3 - x0;
1670  double y30 = y3 - y0;
1671  double z30 = z3 - z0;
1672 
1673  double detJ = x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30;
1674  return detJ * 0.1666666666666666666667;
1675  }
1676 
1677  void ComputeGaussPointPositions(Geometry< Node >& geom, BoundedMatrix<double, 4, 3 > & pos, BoundedMatrix<double, 4, 3 > & N)
1678  {
1679  double one_third = 1.0 / 3.0;
1680  double one_sixt = 1.0 / 6.0;
1681  double two_third = 2.0 * one_third;
1682 
1683  N(0, 0) = one_sixt;
1684  N(0, 1) = one_sixt;
1685  N(0, 2) = two_third;
1686  N(1, 0) = two_third;
1687  N(1, 1) = one_sixt;
1688  N(1, 2) = one_sixt;
1689  N(2, 0) = one_sixt;
1690  N(2, 1) = two_third;
1691  N(2, 2) = one_sixt;
1692  N(3, 0) = one_third;
1693  N(3, 1) = one_third;
1694  N(3, 2) = one_third;
1695 
1696 
1697  //first
1698  pos(0, 0) = one_sixt * geom[0].X() + one_sixt * geom[1].X() + two_third * geom[2].X();
1699  pos(0, 1) = one_sixt * geom[0].Y() + one_sixt * geom[1].Y() + two_third * geom[2].Y();
1700  pos(0, 2) = one_sixt * geom[0].Z() + one_sixt * geom[1].Z() + two_third * geom[2].Z();
1701 
1702  //second
1703  pos(1, 0) = two_third * geom[0].X() + one_sixt * geom[1].X() + one_sixt * geom[2].X();
1704  pos(1, 1) = two_third * geom[0].Y() + one_sixt * geom[1].Y() + one_sixt * geom[2].Y();
1705  pos(1, 2) = two_third * geom[0].Z() + one_sixt * geom[1].Z() + one_sixt * geom[2].Z();
1706 
1707  //third
1708  pos(2, 0) = one_sixt * geom[0].X() + two_third * geom[1].X() + one_sixt * geom[2].X();
1709  pos(2, 1) = one_sixt * geom[0].Y() + two_third * geom[1].Y() + one_sixt * geom[2].Y();
1710  pos(2, 2) = one_sixt * geom[0].Z() + two_third * geom[1].Z() + one_sixt * geom[2].Z();
1711 
1712  //fourth
1713  pos(3, 0) = one_third * geom[0].X() + one_third * geom[1].X() + one_third * geom[2].X();
1714  pos(3, 1) = one_third * geom[0].Y() + one_third * geom[1].Y() + one_third * geom[2].Y();
1715  pos(3, 2) = one_third * geom[0].Z() + one_third * geom[1].Z() + one_third * geom[2].Z();
1716 
1717  }
1718 
1719  void ComputeGaussPointPositions(Geometry< Node >& geom, BoundedMatrix<double, 16, 3 > & pos, BoundedMatrix<double, 16, 3 > & N)
1720  {
1721  //lower diagonal terms
1722  double ypos = 1.0 / 12.0;
1723  int pos_counter = 0;
1724  for (unsigned int i = 0; i < 4; i++)
1725  {
1726  double xpos = 1.0 / 12.0;
1727  for (unsigned int j = 0; j < 4 - i; j++)
1728  {
1729  double N1 = xpos;
1730  double N2 = ypos;
1731  double N3 = 1.0 - xpos - ypos;
1732 
1733  pos(pos_counter, 0) = N1 * geom[0].X() + N2 * geom[1].X() + N3 * geom[2].X();
1734  pos(pos_counter, 1) = N1 * geom[0].Y() + N2 * geom[1].Y() + N3 * geom[2].Y();
1735  pos(pos_counter, 2) = N1 * geom[0].Z() + N2 * geom[1].Z() + N3 * geom[2].Z();
1736 
1737  N(pos_counter, 0) = N1;
1738  N(pos_counter, 1) = N2;
1739  N(pos_counter, 2) = N3;
1740 
1741  xpos += 1.0 / 4.0;
1742  pos_counter += 1;
1743 
1744  }
1745  ypos += 1.0 / 4.0;
1746  }
1747 
1748  //lower diagonal terms
1749  ypos = 2.0 / 12.0;
1750  // pos_counter = 8;
1751  for (unsigned int i = 0; i < 3; i++)
1752  {
1753  double xpos = 2.0 / 12.0;
1754  for (unsigned int j = 0; j < 4 - i; j++)
1755  {
1756  double N1 = xpos;
1757  double N2 = ypos;
1758  double N3 = 1.0 - xpos - ypos;
1759 
1760  pos(pos_counter, 0) = N1 * geom[0].X() + N2 * geom[1].X() + N3 * geom[2].X();
1761  pos(pos_counter, 1) = N1 * geom[0].Y() + N2 * geom[1].Y() + N3 * geom[2].Y();
1762  pos(pos_counter, 2) = N1 * geom[0].Z() + N2 * geom[1].Z() + N3 * geom[2].Z();
1763 
1764  N(pos_counter, 0) = N1;
1765  N(pos_counter, 1) = N2;
1766  N(pos_counter, 2) = N3;
1767 
1768  xpos += 1.0 / 4.0;
1769  pos_counter += 1;
1770 
1771  }
1772  ypos += 1.0 / 4.0;
1773  }
1774  }
1775 
1776  void ConsistentMassMatrix(const double A, BoundedMatrix<double, 3, 3 > & M)
1777  {
1778  double c1 = A / 12.0;
1779  double c2 = 2.0 * c1;
1780  M(0, 0) = c2;
1781  M(0, 1) = c1;
1782  M(0, 2) = c1;
1783  M(1, 0) = c1;
1784  M(1, 1) = c2;
1785  M(1, 2) = c1;
1786  M(2, 0) = c1;
1787  M(2, 1) = c1;
1788  M(2, 2) = c2;
1789  }
1790 
1791  void CalculateInterfaceNormal(BoundedMatrix<double, 3, 2 >& rPoints, array_1d<double,3>& rDistances, array_1d<double,2>& normal, double & interface_area, array_1d<double,3>& Ninterface, BoundedMatrix<double, 2, 2 >& rInterfacePoints)
1792  {
1793  double sign_correction=1.0;
1794 
1795 
1796 
1797  BoundedMatrix<double, 2, 2 > InterfacePoints;
1798 
1799  array_1d<bool,3> cut_edges;
1800 
1801  array_1d<double,2> interface_segment=ZeroVector(2);
1802 
1803  if ((rDistances(0)*rDistances(1))<0.0) cut_edges[0]=true;//edge 12 is cut
1804 
1805  else cut_edges[0]=false;
1806 
1807 
1808 
1809  if ((rDistances(1)*rDistances(2))<0.0) cut_edges[1]=true;//edge 23 is cut.
1810 
1811  else cut_edges[1]=false;
1812 
1813 
1814 
1815  if ((rDistances(2)*rDistances(0))<0.0) cut_edges[2]=true;//edge 13 is cut.
1816 
1817  else cut_edges[2]=false;
1818 
1819 
1820 
1821 
1822 
1823  if (cut_edges[0])
1824 
1825  {
1826 
1827  if (rDistances(0)>0.0) sign_correction=1.0;
1828 
1829  else sign_correction=-1.0;
1830 
1831 
1832 
1833  const double relative_position = abs(rDistances(1)/(rDistances(1)-rDistances(0) ) );
1834 
1835  InterfacePoints(0,0) = relative_position*rPoints(0,0) + (1.0-relative_position)*rPoints(1,0);
1836 
1837  InterfacePoints(0,1) = relative_position*rPoints(0,1) + (1.0-relative_position)*rPoints(1,1);
1838 
1839 
1840 
1841  if (cut_edges[1])
1842 
1843  {
1844 
1845  const double relative_position2 = abs(rDistances(2)/(rDistances(1)-rDistances(2) ) );
1846 
1847  InterfacePoints(1,0) = relative_position2*rPoints(1,0) + (1.0-relative_position2)*rPoints(2,0);
1848 
1849  InterfacePoints(1,1) = relative_position2*rPoints(1,1) + (1.0-relative_position2)*rPoints(2,1);
1850 
1851  }
1852 
1853  else
1854 
1855  {
1856 
1857  const double relative_position2 = abs(rDistances(0)/(rDistances(2)-rDistances(0) ) );
1858 
1859  InterfacePoints(1,0) = relative_position2*rPoints(2,0) + (1.0-relative_position2)*rPoints(0,0);
1860 
1861  InterfacePoints(1,1) = relative_position2*rPoints(2,1) + (1.0-relative_position2)*rPoints(0,1);
1862 
1863  }
1864 
1865  }
1866 
1867  else
1868 
1869  {
1870 
1871  if (rDistances(1)>0.0) sign_correction=1.0;
1872 
1873  else sign_correction=-1.0;
1874 
1875 
1876 
1877  const double relative_position = abs(rDistances(2)/(rDistances(2)-rDistances(1) ) );
1878 
1879  InterfacePoints(0,0) = relative_position*rPoints(1,0) + (1.0-relative_position)*rPoints(2,0);
1880 
1881  InterfacePoints(0,1) = relative_position*rPoints(1,1) + (1.0-relative_position)*rPoints(2,1);
1882 
1883 
1884 
1885  const double relative_position2 = abs(rDistances(0)/(rDistances(2)-rDistances(0) ) );
1886 
1887  InterfacePoints(1,0) = relative_position2*rPoints(2,0) + (1.0-relative_position2)*rPoints(0,0);
1888 
1889  InterfacePoints(1,1) = relative_position2*rPoints(2,1) + (1.0-relative_position2)*rPoints(0,1);
1890 
1891  }
1892 
1893  interface_segment[0] = (InterfacePoints(1,0)-InterfacePoints(0,0));
1894 
1895  interface_segment[1] = (InterfacePoints(1,1)-InterfacePoints(0,1));
1896 
1897 
1898  const double norm = sqrt( pow((interface_segment[0]),2) + pow((interface_segment[1]),2));
1899 
1900 
1901 
1902  normal(0)= -interface_segment[1]*sign_correction/norm;
1903 
1904  normal(1)= interface_segment[0]*sign_correction/norm;
1905 
1906  //KRATOS_WATCH(interface_segment)
1907 
1908  //KRATOS_WATCH(InterfacePoints)
1909 
1910  interface_area=norm;
1911 
1912  rInterfacePoints(0,0)=InterfacePoints(0,0);
1913  rInterfacePoints(0,1)=InterfacePoints(0,1);
1914  rInterfacePoints(1,0)=InterfacePoints(1,0);
1915  rInterfacePoints(1,1)=InterfacePoints(1,1);
1916 
1917  const double x_interface = 0.5*(InterfacePoints(0,0)+InterfacePoints(1,0));
1918 
1919  const double y_interface = 0.5*(InterfacePoints(0,1)+InterfacePoints(1,1));
1920 
1921  // CalculatePosition(rPoints, x_interface, y_interface, 0.0, Ninterface );
1922 
1924  double x0 = rPoints(0,0);
1925  double y0 = rPoints(0,1);
1926  double x1 = rPoints(1,0);
1927  double y1 = rPoints(1,1);
1928  double x2 = rPoints(2,0);
1929  double y2 = rPoints(2,1);
1930  double area = CalculateVol(x0, y0, x1, y1, x2, y2);
1931  double inv_area = 0.0;
1932  if (area == 0.0)
1933  {
1934  KRATOS_ERROR<<"element with zero area found";
1935  }
1936  else
1937  {
1938  inv_area = 1.0 / area;
1939  }
1940 
1941  Ninterface[0]= CalculateVol(x1, y1, x2, y2, x_interface, y_interface) * inv_area;
1942  Ninterface[1] = CalculateVol(x2, y2, x0, y0, x_interface, y_interface) * inv_area;
1943  Ninterface[2] = CalculateVol(x0, y0, x1, y1, x_interface, y_interface) * inv_area;
1944 
1945  }
1946 
1947  bool CalculatePosition(Geometry<Node >&geom,const double xc, const double yc, const double zc,array_1d<double, 3 > & N )
1948  {
1949  double x0 = geom[0].X();
1950  double y0 = geom[0].Y();
1951  double x1 = geom[1].X();
1952  double y1 = geom[1].Y();
1953  double x2 = geom[2].X();
1954  double y2 = geom[2].Y();
1955 
1956  double area = CalculateVol(x0, y0, x1, y1, x2, y2);
1957  double inv_area = 0.0;
1958  if (area == 0.0)
1959  {
1960  KRATOS_ERROR<<"element with zero area found";
1961  }
1962  else
1963  {
1964  inv_area = 1.0 / area;
1965  }
1966 
1967 
1968  N[0] = CalculateVol(x1, y1, x2, y2, xc, yc) * inv_area;
1969  N[1] = CalculateVol(x2, y2, x0, y0, xc, yc) * inv_area;
1970  N[2] = CalculateVol(x0, y0, x1, y1, xc, yc) * inv_area;
1971 
1972 
1973  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
1974  return true;
1975 
1976  return false;
1977  }
1978 
1979  template<class T>
1980  bool InvertMatrix(const T& input, T& inverse)
1981 
1982  {
1983  typedef permutation_matrix<std::size_t> pmatrix;
1984 
1985  // create a working copy of the input
1986 
1987  T A(input);
1988 
1989  // create a permutation matrix for the LU-factorization
1990 
1991  pmatrix pm(A.size1());
1992 
1993  // perform LU-factorization
1994 
1995  int res = lu_factorize(A, pm);
1996 
1997  if (res != 0)
1998 
1999  return false;
2000 
2001  // create identity matrix of "inverse"
2002 
2003  inverse.assign(identity_matrix<double> (A.size1()));
2004 
2005  // backsubstitute to get the inverse
2006 
2007  lu_substitute(A, pm, inverse);
2008 
2009  return true;
2010 
2011 
2012  }
2013 
2014  };
2015 
2016 } // namespace Kratos.
2017 
2018 #endif // KRATOS_LAGRANGIAN_PARTICLES_UTILITIES_INCLUDED defined
binbased_fast_point_locator.h
Kratos::BinBasedFastPointLocator
This class is designed to allow the fast location of MANY points on the top of a 3D mesh.
Definition: binbased_fast_point_locator.h:68
Kratos::BinBasedFastPointLocator::ResultIteratorType
ConfigureType::ResultIteratorType ResultIteratorType
Definition: binbased_fast_point_locator.h:82
Kratos::BinBasedFastPointLocator::ResultContainerType
ConfigureType::ResultContainerType ResultContainerType
Definition: binbased_fast_point_locator.h:81
Kratos::Bucket
Short class definition.
Definition: bucket.h:57
Kratos::DistanceCalculator1
Definition: particle_utilities.h:52
Kratos::DistanceCalculator1::operator()
double operator()(T const &p1, T const &p2)
Definition: particle_utilities.h:55
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::size
SizeType size() const
Definition: geometry.h:518
Kratos::Geometry::GetValue
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometry.h:627
Kratos::GeometryUtils::CalculateGeometryData
static void CalculateGeometryData(const GeometryType &rGeometry, BoundedMatrix< double, 4, 3 > &rDN_DX, array_1d< double, 4 > &rN, double &rVolume)
This function is designed to compute the shape function derivatives, shape functions and volume in 3D...
Definition: geometry_utilities.h:176
Kratos::GlobalPointersVector::iterator
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: global_pointers_vector.h:79
Kratos::GlobalPointersVector::begin
iterator begin()
Definition: global_pointers_vector.h:221
Kratos::GlobalPointersVector::end
iterator end()
Definition: global_pointers_vector.h:229
Kratos::Internals::Matrix
Definition: amatrix_interface.h:41
Kratos::MathUtils::CrossProduct
static T CrossProduct(const T &a, const T &b)
Performs the vector product of the two input vectors a,b.
Definition: math_utils.h:762
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::ElementsBegin
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
Kratos::ModelPart::AddNode
void AddNode(NodeType::Pointer pNewNode, IndexType ThisIndex=0)
Definition: model_part.cpp:211
Kratos::ModelPart::OverwriteSolutionStepData
void OverwriteSolutionStepData(IndexType SourceSolutionStepIndex, IndexType DestinationSourceSolutionStepIndex)
Definition: model_part.cpp:180
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::NodeIterator
MeshType::NodeIterator NodeIterator
Definition: model_part.h:134
Kratos::ModelPart::Elements
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
Kratos::ModelPart::NodesContainerType
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
Kratos::ModelPart::ElementsEnd
ElementIterator ElementsEnd(IndexType ThisIndex=0)
Definition: model_part.h:1179
Kratos::ModelPart::NodesEnd
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
Kratos::ModelPart::Nodes
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::Parameters::SetValue
void SetValue(const std::string &rEntry, const Parameters &rOtherValue)
This method sets an existing parameter with a given parameter.
Definition: kratos_parameters.cpp:443
Kratos::ParticleUtils
Definition: particle_utilities.h:69
Kratos::ParticleUtils::DetectAllOilClusters
void DetectAllOilClusters(ModelPart &mp_local_model_part)
Definition: particle_utilities.h:1022
Kratos::ParticleUtils::TransferToEulerianMesh
void TransferToEulerianMesh(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart)
Definition: particle_utilities.h:260
Kratos::ParticleUtils::EstimateTime
void EstimateTime(ModelPart &rEulerianModelPart, const double max_dt)
Definition: particle_utilities.h:74
Kratos::ParticleUtils::TransferToEulerianMeshShapeBased_aux
void TransferToEulerianMeshShapeBased_aux(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart, BinBasedFastPointLocator< TDim > &node_locator)
Definition: particle_utilities.h:375
Kratos::ParticleUtils::TransferToEulerianMeshShapeBased
void TransferToEulerianMeshShapeBased(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart, BinBasedFastPointLocator< TDim > &node_locator)
Definition: particle_utilities.h:1407
Kratos::ParticleUtils::TransferToParticlesAirVelocity
void TransferToParticlesAirVelocity(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart, BinBasedFastPointLocator< TDim > &node_locator)
Definition: particle_utilities.h:913
Kratos::ParticleUtils::MarkExcessivelyCloseNodes
void MarkExcessivelyCloseNodes(ModelPart::NodesContainerType &rNodes)
Definition: particle_utilities.h:879
Kratos::ParticleUtils::ColorOilClusters
void ColorOilClusters(ModelPart::NodesContainerType::iterator iNode, const int color)
Definition: particle_utilities.h:1179
Kratos::ParticleUtils::VisualizationModelPart
void VisualizationModelPart(ModelPart &rCompleteModelPart, ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart)
Definition: particle_utilities.h:125
Kratos::ParticleUtils::movethermocouples
void movethermocouples(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart, BinBasedFastPointLocator< TDim > &node_locator)
Definition: particle_utilities.h:1237
Kratos::ParticleUtils::MoveLonelyNodes
void MoveLonelyNodes(ModelPart &ThisModelPart)
Definition: particle_utilities.h:835
Kratos::ParticleUtils::MoveMesh_Streamlines_freesurfaceflows
void MoveMesh_Streamlines_freesurfaceflows(ModelPart &rModelPart, unsigned int substeps)
Definition: particle_utilities.h:720
Kratos::ParticleUtils::TransferToEulerianMesh_Face_Heat_Flux
void TransferToEulerianMesh_Face_Heat_Flux(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart)
Definition: particle_utilities.h:154
Kratos::ParticleUtils::TransferToEulerianMesh_2
void TransferToEulerianMesh_2(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart)
Definition: particle_utilities.h:1290
Kratos::ParticleUtils::TransferToEulerianMeshShapeBased_aux_3D
void TransferToEulerianMeshShapeBased_aux_3D(ModelPart &rEulerianModelPart, ModelPart &rLagrangianModelPart, BinBasedFastPointLocator< TDim > &node_locator)
Definition: particle_utilities.h:653
Kratos::ParticleUtils::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(ParticleUtils< TDim >)
Kratos::ParticleUtils::Calculate_Vol
double Calculate_Vol(ModelPart &rLagrangianModelPart)
Definition: particle_utilities.h:967
Kratos::ParticleUtils::CalculateNewtonianDragCoefficient
void CalculateNewtonianDragCoefficient(const double reynolds, double &drag_coeff)
Definition: particle_utilities.h:1159
Kratos::ParticleUtils::CalculateNormal
void CalculateNormal(ModelPart &full_model_part)
3D
Definition: particle_utilities.h:505
Kratos::ParticleUtils::ComputedDragCoefficient
void ComputedDragCoefficient(double nodal_mass, array_1d< double, 3 > velocity_air, array_1d< double, 3 > velocity_polymer, array_1d< double, 3 > &drag_coefficient)
Definition: particle_utilities.h:1132
Kratos::ParticleUtils::RestartStep
void RestartStep(ModelPart &rModelPart)
Definition: particle_utilities.h:697
Kratos::Point::Y
double Y() const
Definition: point.h:187
Kratos::Point::Z
double Z() const
Definition: point.h:193
Kratos::Point::X
double X() const
Definition: point.h:181
Kratos::SearchStructure
Definition: search_structure.h:309
Kratos::SearchStructure::Set
void Set(IndexVector const &IndexCell, SizeVector const &_MaxSize, IteratorIteratorType const &IteratorBegin)
Definition: search_structure.h:363
Kratos::Tree
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
Kratos::array_1d< double, 3 >
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_WATCH
#define KRATOS_WATCH(variable)
Definition: define.h:806
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
KRATOS_ERROR
#define KRATOS_ERROR
Definition: exception.h:161
kratos_flags.h
model_part.h
DEM_benchmarks.dt
dt
Definition: DEM_benchmarks.py:173
GenerateCN.im
im
Definition: GenerateCN.py:100
Kratos::Python::GetValue
Parameters GetValue(Parameters &rParameters, const std::string &rEntry)
Definition: add_kratos_parameters_to_python.cpp:53
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::norm_2
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::sum
TExpressionType::data_type sum(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:675
Kratos::DistanceIterator
DistanceVector::iterator DistanceIterator
Definition: internal_variables_interpolation_process.h:56
Kratos::PointVector
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::PointTypePointer
PointType::Pointer PointTypePointer
Definition: internal_variables_interpolation_process.h:52
Kratos::PointIterator
PointVector::iterator PointIterator
Definition: internal_variables_interpolation_process.h:54
Kratos::DistanceVector
std::vector< double > DistanceVector
Definition: internal_variables_interpolation_process.h:55
Kratos::PointType
Point PointType
Geometric definitions.
Definition: mortar_classes.h:36
PecletTest.velocity
float velocity
Definition: PecletTest.py:54
all_t_win_vs_m_fixed_error.color
color
Definition: all_t_win_vs_m_fixed_error.py:230
bombardelli_test.coeff
list coeff
Definition: bombardelli_test.py:41
coupling_interface_data.GetSolutionStepValue
def GetSolutionStepValue(entity, variable, solution_step_index)
Definition: coupling_interface_data.py:253
edgebased_PureConvection.dist
float dist
Definition: edgebased_PureConvection.py:89
face_heat.temperature
float temperature
Definition: face_heat.py:58
face_heat.Dt
Dt
Definition: face_heat.py:78
generate_convection_diffusion_explicit_element.res
res
Definition: generate_convection_diffusion_explicit_element.py:211
generate_convection_diffusion_explicit_element.v
v
Definition: generate_convection_diffusion_explicit_element.py:114
generate_droplet_dynamics.h
h
Definition: generate_droplet_dynamics.py:91
generate_frictional_mortar_condition.input
input
Definition: generate_frictional_mortar_condition.py:435
generate_frictional_mortar_condition.x2
x2
Definition: generate_frictional_mortar_condition.py:122
generate_frictional_mortar_condition.x1
x1
Definition: generate_frictional_mortar_condition.py:121
generate_hyper_elastic_simo_taylor_neo_hookean.C
int C
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:27
generate_total_lagrangian_mixed_volumetric_strain_element.tmp
tuple tmp
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:98
isotropic_damage_automatic_differentiation.R
R
Definition: isotropic_damage_automatic_differentiation.py:172
mesh_to_mdpa_converter.radius
float radius
Definition: mesh_to_mdpa_converter.py:18
pure_conduction.vel
vel
Definition: pure_conduction.py:76
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
rotating_cone.sigma
float sigma
Definition: rotating_cone.py:79
rotating_cone.temp
float temp
Definition: rotating_cone.py:85
rotating_cone.xc
float xc
Definition: rotating_cone.py:77
rotating_cone.yc
float yc
Definition: rotating_cone.py:78
script_THERMAL_CORRECT.counter
int counter
Definition: script_THERMAL_CORRECT.py:218
script.dummy
dummy
Definition: script.py:194
sensitivityMatrix.A
A
Definition: sensitivityMatrix.py:70
sensitivityMatrix.N
N
Definition: sensitivityMatrix.py:29
sensitivityMatrix.dim
int dim
Definition: sensitivityMatrix.py:25
tensors::i
integer i
Definition: TensorModule.f:17
test_pureconvectionsolver_benchmarking.zero
zero
Definition: test_pureconvectionsolver_benchmarking.py:94
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
pfem_2_application_variables.h
spatial_containers.h
tetrahedra_3d_4.h
Generated on Sat Feb 3 2024 02:06:39 for KratosMultiphysics by doxygen 1.9.1