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.
move_particle_utility_pfem2.h
Go to the documentation of this file.
1 /*
2 ==============================================================================
3 KratosIncompressibleFluidApplication
4 A library based on:
5 Kratos
6 A General Purpose Software for Multi-Physics Finite Element Analysis
7 Version 1.0 (Released on march 05, 2007).
8 
9 Copyright 2007
10 Pooyan Dadvand, Riccardo Rossi
11 pooyan@cimne.upc.edu
12 rrossi@cimne.upc.edu
13 - CIMNE (International Center for Numerical Methods in Engineering),
14 Gran Capita' s/n, 08034 Barcelona, Spain
15 
16 
17 Permission is hereby granted, free of charge, to any person obtaining
18 a copy of this software and associated documentation files (the
19 "Software"), to deal in the Software without restriction, including
20 without limitation the rights to use, copy, modify, merge, publish,
21 distribute, sublicense and/or sell copies of the Software, and to
22 permit persons to whom the Software is furnished to do so, subject to
23 the following condition:
24 
25 Distribution of this code for any commercial purpose is permissible
26 ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNERS.
27 
28 The above copyright notice and this permission notice shall be
29 included in all copies or substantial portions of the Software.
30 
31 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
32 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
33 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
34 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
35 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
36 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
37 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
38 
39 ==============================================================================
40 */
41 
42 //
43 // Project Name: Kratos
44 // Last Modified by: $Author: pbecker $
45 // Date: $Date: 2011-09-21 12:30:32 $
46 // Revision: $Revision: 1.0 $
47 //
48 //
49 
50 
51 #if !defined(KRATOS_MOVE_PARTICLE_UTILITY_PFEM2_INCLUDED)
52 #define KRATOS_MOVE_PARTICLE_UTILITY_FLUID_PFEM2_INCLUDED
53 
54 
55 
56 // System includes
57 #include <string>
58 #include <iostream>
59 #include <algorithm>
60 
61 // External includes
62 
63 
64 // Project includes
65 #include "includes/define.h"
66 #include "includes/node.h"
67 
69 #include "includes/dof.h"
70 #include "includes/variables.h"
71 #include "includes/cfd_variables.h"
72 #include "includes/deprecated_variables.h"
73 #include "includes/global_pointer_variables.h"
74 #include "containers/array_1d.h"
75 #include "containers/data_value_container.h"
76 #include "includes/mesh.h"
77 #include "utilities/math_utils.h"
79 
80 #include "utilities/geometry_utilities.h"
81 
82 #include "includes/model_part.h"
83 
84 
85 #include "spatial_containers/spatial_containers.h"
86 #include "spatial_containers/cell.h"
87 #include "spatial_containers/bins_dynamic_objects.h"
88 
89 #include "utilities/spatial_containers_configure.h"
90 
91 #include "geometries/line_2d_2.h"
92 #include "geometries/triangle_2d_3.h"
93 #include "geometries/triangle_3d_3.h"
94 #include "geometries/point.h"
95 
96 #include "pfem_2_application_variables.h"
97 #include "pfem_particle_fluidonly.h"
98 
99 //#include "utilities/enrich_2d_2dofs.h"
100 #include "utilities/enrichment_utilities.h"
101 #include "utilities/openmp_utils.h"
102 
103 #include "time.h"
104 
105 //#include "processes/process.h"
106 
107 namespace Kratos
108 {
109  //this class is to be modified by the user to customize the interpolation process
110  template< unsigned int TDim>
111  class MoveParticleUtilityPFEM2
112  {
113  public:
114 
115  typedef SpatialContainersConfigure<TDim> Configure;
116  typedef typename Configure::PointType PointType;
117  //typedef PointType::CoordinatesArrayType CoordinatesArrayType;
118  typedef typename Configure::ContainerType ContainerType;
119  //typedef Configure::PointerType PointerType;
120  typedef typename Configure::IteratorType IteratorType;
121  typedef typename Configure::ResultContainerType ResultContainerType;
122  //typedef Configure::ResultPointerType ResultPointerType;
123  typedef typename Configure::ResultIteratorType ResultIteratorType;
124  typedef PointerVector< PFEM_Particle_Fluid, PFEM_Particle_Fluid*, std::vector<PFEM_Particle_Fluid*> > ParticlePointerVector;
125  //typedef Configure::ContactPairType ContactPairType;
126  //typedef Configure::ContainerContactType ContainerContactType;
127  //typedef Configure::IteratorContactType IteratorContactType;
128  //typedef Configure::PointerContactType PointerContactType;
129  //typedef Configure::PointerTypeIterator PointerTypeIterator;
130 
131  KRATOS_CLASS_POINTER_DEFINITION(MoveParticleUtilityPFEM2);
132 
133  //template<unsigned int TDim>
134  MoveParticleUtilityPFEM2(ModelPart& model_part, int maximum_number_of_particles)
135  : mr_model_part(model_part) , mmaximum_number_of_particles(maximum_number_of_particles)
136  {
137  KRATOS_INFO("MoveParticleUtilityPfem2") << "Initializing utility" << std::endl;
138 
139  Check();
140 
141 
142  //tools to move the domain, in case we are using a moving domain approach.
143  mintialized_transfer_tool=false;
144  mcalculation_domain_complete_displacement=ZeroVector(3);
145  mcalculation_domain_added_displacement=ZeroVector(3);
146 
147  //storing water and air density and their inverses, just in case it is needed for the streamline integration
148  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
149  mDENSITY_AIR = CurrentProcessInfo[DENSITY_AIR];
150  mDENSITY_WATER = CurrentProcessInfo[DENSITY_WATER];
151 
152  //mmaximum_number_of_particles = maximum_number_of_particles;
153 
154  //loop in elements to change their ID to their position in the array. Easier to get information later.
155  //DO NOT PARALELIZE THIS! IT MUST BE SERIAL!!!!!!!!!!!!!!!!!!!!!!
156  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
157  for(unsigned int ii=0; ii<mr_model_part.Elements().size(); ii++)
158  {
159  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
160  ielem->SetId(ii+1);
161  }
162  mlast_elem_id= (mr_model_part.ElementsEnd()-1)->Id();
163  int node_id=0;
164  // we look for the smallest edge. could be used as a weighting function when going lagrangian->eulerian instead of traditional shape functions(method currently used)
165  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
166  vector<unsigned int> node_partition;
167  #ifdef _OPENMP
168  int number_of_threads = omp_get_max_threads();
169  #else
170  int number_of_threads = 1;
171  #endif
172  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
173 
174  #pragma omp parallel for
175  for(int kkk=0; kkk<number_of_threads; kkk++)
176  {
177  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
178  {
179  ModelPart::NodesContainerType::iterator pnode = inodebegin+ii;
180  array_1d<double,3> position_node;
181  double distance=0.0;
182  position_node = pnode->Coordinates();
183  GlobalPointersVector< Node >& rneigh = pnode->GetValue(NEIGHBOUR_NODES);
184  //we loop all the nodes to check all the edges
185  const double number_of_neighbours = double(rneigh.size());
186  for( GlobalPointersVector<Node >::iterator inode = rneigh.begin(); inode!=rneigh.end(); inode++)
187  {
188  array_1d<double,3> position_difference;
189  position_difference = inode->Coordinates() - position_node;
190  double current_distance= sqrt(pow(position_difference[0],2)+pow(position_difference[1],2)+pow(position_difference[2],2));
191  //if (current_distance>distance)
192  // distance=current_distance;
193  distance += current_distance / number_of_neighbours;
194  }
195  //and we save the largest edge.
196  pnode->FastGetSolutionStepValue(MEAN_SIZE)=distance;
197 
198  node_id=pnode->GetId();
199  }
200  }
201  mlast_node_id=node_id;
202 
203  //we also calculate the element mean size in the same way, for the courant number
204  //also we set the right size to the LHS column for the pressure enrichments, in order to recover correctly the enrichment pressure
205  vector<unsigned int> element_partition;
206  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
207 
208  //before doing anything we must reset the vector of nodes contained by each element (particles that are inside each element.
209  #pragma omp parallel for
210  for(int kkk=0; kkk<number_of_threads; kkk++)
211  {
212  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
213  {
214  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
215 
216  double elem_size;
217  array_1d<double,3> Edge(3,0.0);
218  Edge = ielem->GetGeometry()[1].Coordinates() - ielem->GetGeometry()[0].Coordinates();
219  elem_size = Edge[0]*Edge[0];
220  for (unsigned int d = 1; d < TDim; d++)
221  elem_size += Edge[d]*Edge[d];
222 
223  for (unsigned int i = 2; i < (TDim+1); i++)
224  for(unsigned int j = 0; j < i; j++)
225  {
226  Edge = ielem->GetGeometry()[i].Coordinates() - ielem->GetGeometry()[j].Coordinates();
227  double Length = Edge[0]*Edge[0];
228  for (unsigned int d = 1; d < TDim; d++)
229  Length += Edge[d]*Edge[d];
230  if (Length < elem_size) elem_size = Length;
231  }
232  elem_size = sqrt(elem_size);
233  ielem->SetValue(MEAN_SIZE, elem_size);
234 
235  //and the matrix column for the enrichments in the pressure.
236  if constexpr (TDim==3)
237  ielem->SetValue(ENRICH_LHS_ROW_3D, ZeroVector(4));
238  // {
239  // Vector & lhs_enrich = ielem->GetValue(ENRICH_LHS_ROW_3D);
240  // lhs_enrich.resize(4);
241  // lhs_enrich=ZeroVector(4);
242  // }
243  else
244  ielem->SetValue(ENRICH_LHS_ROW, ZeroVector(3));
245  //KRATOS_WATCH(mElemSize)
246  }
247  }
248 
249 
250  //matrix containing the position of the 4/15/45 particles that we will seed at the beggining
251  BoundedMatrix<double, 5*(1+TDim), 3 > pos;
252  BoundedMatrix<double, 5*(1+TDim), (1+TDim) > N;
253 
254  int particle_id=0;
255  mnelems = mr_model_part.Elements().size();
256 
257  KRATOS_INFO("MoveParticleUtilityPfem2") << "About to resize vectors" << std::endl;
258 
259 
260  //setting the right size to the vector containing the particles assigned to each element
261  //particles vector. this vector contains ALL the particles in the simulation.
262  mparticles_vector.resize(mnelems*mmaximum_number_of_particles);
263 
264  //and this vector contains the current number of particles that are in each element (currently zero)
265  mnumber_of_particles_in_elems.resize(mnelems);
266  mnumber_of_particles_in_elems=ZeroVector(mnelems);
267 
268  //when moving the particles, an auxiliary vector is necessary (to store the previous number)
269  mnumber_of_particles_in_elems_aux.resize(mnelems);
270 
271  //each element will have a list of pointers to all the particles that are inside.
272  //this vector contains the pointers to the vector of (particle) pointers of each element.
273  mpointers_to_particle_pointers_vectors.resize(mnelems);
274  KRATOS_INFO("MoveParticleUtilityPfem2") << "About to create particles" << std::endl;
275  //now we seed: LOOP IN ELEMENTS
276  //using loop index, DO NOT paralelize this! change lines : mparticles_in_elems_pointers((ii*mmaximum_number_of_particles)+mparticles_in_elems_integers(ii)) = pparticle; and the next one
277 
278  for(unsigned int ii=0; ii<mr_model_part.Elements().size(); ii++)
279  {
280  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
281  ielem->SetValue(FLUID_PARTICLE_POINTERS, ParticlePointerVector( mmaximum_number_of_particles*2) );//, &firstparticle );
282  ParticlePointerVector& particle_pointers = ielem->GetValue(FLUID_PARTICLE_POINTERS);
283  //now we link the mpointers_to_particle_pointers_vectors to the corresponding element
284  mpointers_to_particle_pointers_vectors(ii) = &particle_pointers;
285  //now we resize the vector of particle pointers. it is double sized because we move the particles from an initial position (first half) to a final position (second half).
286  //for(int j=0; j<(mmaximum_number_of_particles*2); j++)
287  // particle_pointers.push_back(&firstparticle);
288  int & number_of_particles = ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
289  number_of_particles=0;
290  //int & number_of_water_particles = ielem->GetValue(NUMBER_OF_WATER_PARTICLES);
291 
292  Geometry< Node >& geom = ielem->GetGeometry();
293  //unsigned int elem_id = ielem->Id();
294  //mareas_vector[i_int]=CalculateArea(geom); UNUSED SO COMMENTED
295  ComputeGaussPointPositions_initial(geom, pos, N); //we also have the standard (4), and 45
296  //now we seed the particles in the current element
297  for (unsigned int j = 0; j < pos.size1(); j++)
298  {
299  ++particle_id;
300 
301  PFEM_Particle_Fluid& pparticle = mparticles_vector[particle_id-1];
302  pparticle.X()=pos(j,0);
303  pparticle.Y()=pos(j,1);
304  pparticle.Z()=pos(j,2);
305 
306  pparticle.GetEraseFlag()=false;
307 
308  array_1d<float, 3 > & vel = pparticle.GetVelocity();
309  float & distance = pparticle.GetDistance();
310  noalias(vel) = ZeroVector(3);
311  distance = 0.0;
312 
313  for (unsigned int k = 0; k < (TDim+1); k++)
314  {
315  noalias(vel) += (N(j, k) * geom[k].FastGetSolutionStepValue(VELOCITY));
316  distance += N(j, k) * geom[k].FastGetSolutionStepValue(DISTANCE);
317  }
318 
319  if (distance <= 0.0)
320  distance = -1.0;
321  else
322  distance = 1.0;
323 
324  particle_pointers(j) = &pparticle;
325  number_of_particles++;
326  }
327  }
328 
329 
330  bool nonzero_mesh_velocity = false;
331  //seeing if we have to use the mesh_velocity or not
332  for(ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
333  inode!=mr_model_part.NodesEnd(); inode++)
334  {
335  const array_1d<double, 3 > velocity = inode->FastGetSolutionStepValue(MESH_VELOCITY);
336  for(unsigned int i = 0; i!=3; i++)
337  {
338  if (fabs(velocity[i])>1.0e-9)
339  nonzero_mesh_velocity=true;
340  }
341  if( nonzero_mesh_velocity==true)
342  break;
343  }
344 
345  if ( nonzero_mesh_velocity==true)
346  muse_mesh_velocity_to_convect = true; // if there is mesh velocity, then we have to take it into account when moving the particles
347  else
348  muse_mesh_velocity_to_convect = false; //otherwise, we can avoid reading the values since we know it is zero everywhere (to save time!)
349 
350 
351 
352  m_nparticles=particle_id; //we save the last particle created as the total number of particles we have. For the moment this is true.
353  KRATOS_INFO("MoveParticleUtilityPfem2") << "Number of particles created : " << m_nparticles << std::endl;
354  mparticle_printing_tool_initialized=false;
355  }
356 
357 
358  ~MoveParticleUtilityPFEM2()
359  {}
360 
361  void MountBin()
362  {
363  KRATOS_TRY
364 
365  //copy the elements to a new container, as the list will
366  //be shuffled duringthe construction of the tree
367  ContainerType& rElements = mr_model_part.ElementsArray();
368  IteratorType it_begin = rElements.begin();
369  IteratorType it_end = rElements.end();
370  //const int number_of_elem = rElements.size();
371 
372  typename BinsObjectDynamic<Configure>::Pointer paux = typename BinsObjectDynamic<Configure>::Pointer(new BinsObjectDynamic<Configure>(it_begin, it_end ) );
373  paux.swap(mpBinsObjectDynamic);
374  //BinsObjectDynamic<Configure> mpBinsObjectDynamic(it_begin, it_end );
375 
376  KRATOS_INFO("MoveParticleUtilityPfem2") << "Finished mounting Bins" << std::endl;
377 
378  KRATOS_CATCH("")
379  }
380 
381 
382  //TOOL TO TRANSFER INFORMATION INITIALLY FROM ONE DOMAIN TO OTHER.
383  void IntializeTransferTool(ModelPart* topographic_model_part, array_1d<double, 3 > initial_domains_offset, bool ovewrite_particle_data)
384  //mtopographic_model_part(topographic_model_part)
385  {
386  KRATOS_TRY
387 
388  mintialized_transfer_tool=true;
389  const unsigned int max_results = 1000;
390  std::cout << "initializing transfer utility" << std::endl;
391  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
392  mcalculation_domain_complete_displacement=initial_domains_offset;
393 
394  mtopographic_model_part_pointer = topographic_model_part; //copying the pointer.
395 
396  //CONSTRUCTING BIN STRUCTURE
397  ContainerType& rElements_topo = mtopographic_model_part_pointer->ElementsArray();
398  IteratorType it_begin_topo = rElements_topo.begin();
399  IteratorType it_end_topo = rElements_topo.end();
400  typename BinsObjectDynamic<Configure>::Pointer paux = typename BinsObjectDynamic<Configure>::Pointer(new BinsObjectDynamic<Configure>(it_begin_topo, it_end_topo ) );
401  paux.swap(mpTopographicBinsObjectDynamic);
402 
403 
404  std::cout << "Gathering Information From Topographic Domain for the first time" << std::endl;
405  if(ovewrite_particle_data==false)
406  {
407  std::cout << "Not overwriting particle data (assuming correct initial conditions in calculation domain)" << std::endl;
408  }
409  else
410  {
411  std::cout << "Replacing particle information using the Topographic domain" << std::endl;
412  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
413  //KRATOS_WATCH(offset) //(flag managed only by MoveParticles
414  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
415  vector<unsigned int> element_partition;
416  #ifdef _OPENMP
417  int number_of_threads = omp_get_max_threads();
418  #else
419  int number_of_threads = 1;
420  #endif
421  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
422 
423  #pragma omp parallel for
424  for(int kkk=0; kkk<number_of_threads; kkk++)
425  {
426  ResultContainerType results(max_results);
427  ResultIteratorType result_begin = results.begin();
428 
429  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
430  {
431  if (results.size()!=max_results)
432  results.resize(max_results);
433  //const int & elem_id = ielem->Id();
434  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
435  Element::Pointer pelement(*it_begin_topo); //we have no idea in which element it might be from the topographic domain, so we just set it in the first element.
436 
437  //Geometry<Node >& geom = ielem->GetGeometry();
438  //array_1d<double,TDim+1> N;
439 
440 
441  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
442  int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
443  //std::cout << "elem " << ii << " with " << (unsigned int)number_of_particles_in_elem << " particles" << std::endl;
444 
445  for (int iii=0; iii<number_of_particles_in_elem ; iii++ )
446  {
447  //KRATOS_WATCH(iii)
448  if (iii>mmaximum_number_of_particles) //it means we are out of our portion of the array, abort loop!
449  break;
450 
451  PFEM_Particle_Fluid & pparticle = element_particle_pointers[offset+iii];
452 
453 
454  bool erase_flag= pparticle.GetEraseFlag();
455  if (erase_flag==false)
456  {
457  OverwriteParticleDataUsingTopographicDomain(pparticle,pelement,mcalculation_domain_complete_displacement,result_begin, max_results);
458 
459  }
460 
461 
462  }
463 
464  }
465  }
466  }
467  KRATOS_CATCH("")
468 
469  }
470 
471  //TOOL TO TRANSFER INFORMATION FROM ONE DOMAIN TO OTHER when necessary. to be don
472  void PreReseedUsingTopographicDomain(const int minimum_number_of_particles, array_1d<double, 3 > domains_added_displacement)
473  //mtopographic_model_part(topographic_model_part)
474  {
475  KRATOS_TRY
476 
477  if(mintialized_transfer_tool==false)
478  KRATOS_THROW_ERROR(std::logic_error, "TRANSFER TOOL NOT INITIALIZED!", "");
479  const unsigned int max_results = 1000;
480  std::cout << "executing transfer tool" << std::endl;
481  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
482  mcalculation_domain_added_displacement = domains_added_displacement;
483  mcalculation_domain_complete_displacement += domains_added_displacement;
484 
485  ContainerType& rElements_topo = mtopographic_model_part_pointer->ElementsArray();
486  IteratorType it_begin_topo = rElements_topo.begin();
487 
488  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
489  //KRATOS_WATCH(offset) //(flag managed only by MoveParticles
490  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
491  vector<unsigned int> element_partition;
492  #ifdef _OPENMP
493  int number_of_threads = omp_get_max_threads();
494  #else
495  int number_of_threads = 1;
496  #endif
497  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
498 
499  #pragma omp parallel for
500  for(int kkk=0; kkk<number_of_threads; kkk++)
501  {
502  ResultContainerType results(max_results);
503  ResultIteratorType result_begin = results.begin();
504 
505  Element::Pointer pelement(*it_begin_topo); //we have no idea in which element it might be from the topographic domain, so we just set it in the first element.
506 
507  BoundedMatrix<double, (TDim+1), 3 > pos;
508  BoundedMatrix<double, (TDim+1) , (TDim+1) > N;
509  unsigned int freeparticle=0; //we start with the first position in the particles array
510 
511  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
512  {
513  if (results.size()!=max_results)
514  results.resize(max_results);
515  //const int & elem_id = ielem->Id();
516  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
517 
518  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
519  int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
520  if (number_of_particles_in_elem<(minimum_number_of_particles))// && (ielem->GetGeometry())[0].Y()<0.10 )
521  {
522  //KRATOS_WATCH("elem with little particles")
523  Geometry< Node >& geom = ielem->GetGeometry();
524  ComputeGaussPointPositionsForPreReseed(geom, pos, N);
525  //double conductivity = ielem->GetProperties()[CONDUCTIVITY];
526  //KRATOS_WATCH(conductivity);
527  for (unsigned int j = 0; j < (pos.size1()); j++) //i am dropping the last one, the one in the middle of the element
528  {
529  bool keep_looking = true;
530  while(keep_looking)
531  {
532  if (mparticles_vector[freeparticle].GetEraseFlag()==true)
533  {
534  #pragma omp critical
535  {
536  if (mparticles_vector[freeparticle].GetEraseFlag()==true)
537  {
538  mparticles_vector[freeparticle].GetEraseFlag()=false;
539  keep_looking=false;
540  }
541  }
542  if (keep_looking==false)
543  break;
544  /*
545  else if (freeparticle<(it_end_particle_model_part-1))
546  freeparticle++;
547  */
548  else
549  freeparticle++;
550  //break;
551  }
552  else
553  {
554  //if (freeparticle<(it_end_particle_model_part-1))
555  freeparticle++;
556  //else
557  //break; //we finished the list and we couldnt find a free space
558  }
559  }
560 
561  PFEM_Particle_Fluid pparticle(pos(j,0),pos(j,1),pos(j,2));
562  /*
563  PFEM_Particle_Fluid & pparticle = mparticles_vector[freeparticle];
564  pparticle.X() = pos(j,0);
565  pparticle.Y() = pos(j,1);
566  pparticle.Z() = pos(j,2);
567  */
568  array_1d<double,TDim+1>aux2_N;
569  bool is_found = CalculatePosition(geom,pos(j,0),pos(j,1),pos(j,2),aux2_N);
570  if (is_found==false)
571  {
572  KRATOS_WATCH(aux2_N);
573  }
574 
575  pparticle.GetEraseFlag()=false;
576  OverwriteParticleDataUsingTopographicDomain(pparticle,pelement,mcalculation_domain_complete_displacement,result_begin, max_results);
577  //and we copy it to the array:
578  mparticles_vector[freeparticle] = pparticle;
579  element_particle_pointers(offset+number_of_particles_in_elem) = &mparticles_vector[freeparticle];
580  number_of_particles_in_elem++;
581 
582  }
583  }
584  }
585  }
586 
587  KRATOS_CATCH("")
588 
589  }
590 
591  void CalculateVelOverElemSize()
592  {
593  KRATOS_TRY
594 
595  //ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
596 
597  const double nodal_weight = 1.0/ (1.0 + double (TDim) );
598 
599  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
600  vector<unsigned int> element_partition;
601  #ifdef _OPENMP
602  int number_of_threads = omp_get_max_threads();
603  #else
604  int number_of_threads = 1;
605  #endif
606  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
607 
608  if (muse_mesh_velocity_to_convect==false)
609  {
610  #pragma omp parallel for
611  for(int kkk=0; kkk<number_of_threads; kkk++)
612  {
613  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
614  {
615  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
616  Geometry<Node >& geom = ielem->GetGeometry();
617 
618  array_1d<double, 3 >vector_mean_velocity=ZeroVector(3);
619 
620  for (unsigned int i=0; i != (TDim+1) ; i++)
621  vector_mean_velocity += geom[i].FastGetSolutionStepValue(VELOCITY);
622  vector_mean_velocity *= nodal_weight;
623 
624  const double mean_velocity = sqrt ( pow(vector_mean_velocity[0],2) + pow(vector_mean_velocity[1],2) + pow(vector_mean_velocity[2],2) );
625  ielem->SetValue(VELOCITY_OVER_ELEM_SIZE, mean_velocity / ( ielem->GetValue(MEAN_SIZE) ) );
626  }
627  }
628  }
629  else
630  {
631  #pragma omp parallel for
632  for(int kkk=0; kkk<number_of_threads; kkk++)
633  {
634  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
635  {
636  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
637  Geometry<Node >& geom = ielem->GetGeometry();
638 
639  array_1d<double, 3 >vector_mean_velocity=ZeroVector(3);
640 
641  for (unsigned int i=0; i != (TDim+1) ; i++)
642  vector_mean_velocity += geom[i].FastGetSolutionStepValue(VELOCITY)-geom[i].FastGetSolutionStepValue(MESH_VELOCITY);
643  vector_mean_velocity *= nodal_weight;
644 
645  const double mean_velocity = sqrt ( pow(vector_mean_velocity[0],2) + pow(vector_mean_velocity[1],2) + pow(vector_mean_velocity[2],2) );
646  ielem->SetValue(VELOCITY_OVER_ELEM_SIZE, mean_velocity / ( ielem->GetValue(MEAN_SIZE) ) );
647  }
648  }
649  }
650  KRATOS_CATCH("")
651  }
652 
653 
654 
655  //name self explained
656  void ResetBoundaryConditions(bool fully_reset_nodes)
657  {
658  KRATOS_TRY
659 
660  if (fully_reset_nodes)
661  {
662  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
663  vector<unsigned int> node_partition;
664  #ifdef _OPENMP
665  int number_of_threads = omp_get_max_threads();
666  #else
667  int number_of_threads = 1;
668  #endif
669  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
670 
671  #pragma omp parallel for
672  for(int kkk=0; kkk<number_of_threads; kkk++)
673  {
674  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
675  {
676  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
677 
678  if (inode->IsFixed(VELOCITY_X))
679  {
680  inode->FastGetSolutionStepValue(VELOCITY_X)=inode->GetSolutionStepValue(VELOCITY_X,1);
681  }
682  if (inode->IsFixed(VELOCITY_Y))
683  {
684  inode->FastGetSolutionStepValue(VELOCITY_Y)=inode->GetSolutionStepValue(VELOCITY_Y,1);
685  }
686  if constexpr (TDim==3)
687  if (inode->IsFixed(VELOCITY_Z))
688  {
689  inode->FastGetSolutionStepValue(VELOCITY_Z)=inode->GetSolutionStepValue(VELOCITY_Z,1);
690  }
691 
692  if (inode->IsFixed(PRESSURE))
693  inode->FastGetSolutionStepValue(PRESSURE)=inode->GetSolutionStepValue(PRESSURE,1);
694  inode->GetSolutionStepValue(PRESSURE,1)=inode->FastGetSolutionStepValue(PRESSURE);
695  }
696  }
697  }
698  else //for fractional step only!
699  {
700  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
701  vector<unsigned int> node_partition;
702  #ifdef _OPENMP
703  int number_of_threads = omp_get_max_threads();
704  #else
705  int number_of_threads = 1;
706  #endif
707  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
708 
709  #pragma omp parallel for
710  for(int kkk=0; kkk<number_of_threads; kkk++)
711  {
712  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
713  {
714  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
715 
716  const array_1d<double, 3 > original_velocity = inode->FastGetSolutionStepValue(VELOCITY);
717 
718  if (inode->IsFixed(VELOCITY_X) || inode->IsFixed(VELOCITY_Y) || inode->IsFixed(VELOCITY_Z) )
719  {
720 
721  const array_1d<double, 3 > & normal = inode->FastGetSolutionStepValue(NORMAL);
722  const double normal_scalar_sq = normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2];
723  const array_1d<double, 3 > normal_adimensionalized = normal / sqrt(normal_scalar_sq);
724  array_1d<double, 3 > & velocity = inode->FastGetSolutionStepValue(VELOCITY);
725 
726  array_1d<double, 3 > normal_velocity;
727  for (unsigned int j=0; j!=3; j++)
728  normal_velocity[j] = fabs(normal_adimensionalized[j])*original_velocity[j];
729 
730  if (inode->IsFixed(VELOCITY_X))
731  {
732  velocity[0] = original_velocity[0] - normal_velocity[0];
733  }
734  if (inode->IsFixed(VELOCITY_Y))
735  {
736  velocity[1] = original_velocity[1] - normal_velocity[1];
737  }
738  if constexpr (TDim==3)
739  if (inode->IsFixed(VELOCITY_Z))
740  {
741  velocity[2] = original_velocity[2] - normal_velocity[2];
742  }
743 
744  }
745 
746  if (inode->IsFixed(PRESSURE))
747  inode->FastGetSolutionStepValue(PRESSURE)=inode->GetSolutionStepValue(PRESSURE,1);
748  }
749  }
750  }
751  KRATOS_CATCH("")
752  }
753 
754  //setting the normal component of the velocity to zero
755  void ResetBoundaryConditionsSlip()
756  {
757  KRATOS_TRY
758 
759  {
760  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
761  vector<unsigned int> node_partition;
762  #ifdef _OPENMP
763  int number_of_threads = omp_get_max_threads();
764  #else
765  int number_of_threads = 1;
766  #endif
767  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
768 
769  #pragma omp parallel for
770  for(int kkk=0; kkk<number_of_threads; kkk++)
771  {
772  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
773  {
774  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
775 
776  if(inode->Is(SLIP))
777  {
778 
779  array_1d<double, 3 >& velocity = inode->FastGetSolutionStepValue(VELOCITY);
780  const array_1d<double, 3 > & normal = inode->FastGetSolutionStepValue(NORMAL);
781  const double normal_scalar_sq = normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2];
782  const array_1d<double, 3 > normal_adimensionalized = normal / sqrt(normal_scalar_sq);
783  //calculating the normal component of the velocity
784  array_1d<double, 3 > normal_velocity;
785  for (unsigned int j=0; j!=3; j++)
786  normal_velocity[j] = normal_adimensionalized[j]*velocity[j];
787 
788  const double dot_prod = normal_velocity[0]*velocity[0] + normal_velocity[1]*velocity[1] + normal_velocity[2]*velocity[2];
789  //if the dot product of velocity * normal velocity is lower than zero, then they have opposite signs and we must invert the direction:
790  if (dot_prod<0.0)
791  normal_velocity*= -1.0;
792 
793  velocity -= normal_velocity; //substracting the normal component
794  }
795  else if (inode->IsFixed(VELOCITY_X) && inode->IsFixed(VELOCITY_Y) )
796  {
797  inode->FastGetSolutionStepValue(VELOCITY) = inode->GetSolutionStepValue(VELOCITY,1);
798  }
799 
800  }
801  }
802  }
803  KRATOS_CATCH("")
804  }
805 
806 
807  void CalculateDeltaVelocity()
808  {
809  KRATOS_TRY
810  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
811  vector<unsigned int> node_partition;
812  #ifdef _OPENMP
813  int number_of_threads = omp_get_max_threads();
814  #else
815  int number_of_threads = 1;
816  #endif
817  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
818 
819  #pragma omp parallel for
820  for(int kkk=0; kkk<number_of_threads; kkk++)
821  {
822  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
823  {
824  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
825  inode->FastGetSolutionStepValue(DELTA_VELOCITY) = inode->FastGetSolutionStepValue(VELOCITY) - inode->FastGetSolutionStepValue(PROJECTED_VELOCITY) ;
826  }
827  }
828 
829  KRATOS_CATCH("")
830  }
831 
832  void CopyVectorVarToPreviousTimeStep(const Variable< array_1d<double, 3 > >& OriginVariable,
833  ModelPart::NodesContainerType& rNodes)
834  {
835  KRATOS_TRY
836  ModelPart::NodesContainerType::iterator inodebegin = rNodes.begin();
837  vector<unsigned int> node_partition;
838  #ifdef _OPENMP
839  int number_of_threads = omp_get_max_threads();
840  #else
841  int number_of_threads = 1;
842  #endif
843  OpenMPUtils::CreatePartition(number_of_threads, rNodes.size(), node_partition);
844 
845  #pragma omp parallel for
846  for(int kkk=0; kkk<number_of_threads; kkk++)
847  {
848  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
849  {
850  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
851  noalias(inode->GetSolutionStepValue(OriginVariable,1)) = inode->FastGetSolutionStepValue(OriginVariable);
852  }
853  }
854  KRATOS_CATCH("")
855  }
856 
857  void CopyScalarVarToPreviousTimeStep(const Variable<double>& OriginVariable,
858  ModelPart::NodesContainerType& rNodes)
859  {
860  KRATOS_TRY
861  ModelPart::NodesContainerType::iterator inodebegin = rNodes.begin();
862  vector<unsigned int> node_partition;
863  #ifdef _OPENMP
864  int number_of_threads = omp_get_max_threads();
865  #else
866  int number_of_threads = 1;
867  #endif
868  OpenMPUtils::CreatePartition(number_of_threads, rNodes.size(), node_partition);
869 
870  #pragma omp parallel for
871  for(int kkk=0; kkk<number_of_threads; kkk++)
872  {
873  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
874  {
875  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
876  inode->GetSolutionStepValue(OriginVariable,1) = inode->FastGetSolutionStepValue(OriginVariable);
877  }
878  }
879  KRATOS_CATCH("")
880  }
881 
882 
883  //to move all the particles across the streamlines. heavy task!
884  void MoveParticles(const bool discriminate_streamlines) //,const bool pressure_gradient_integrate)
885  {
886 
887  KRATOS_TRY
888 
889  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
890 
891  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
892  //moveparticlesdiff reads from the pointers of one part (ie odd) and saves into the other part (ie even part)
893  //since it is the only function in the whole procedure that does this, it must use alternatively one part and the other.
894  //KRATOS_WATCH(offset)
895 
896  bool even_timestep;
897  if (offset!=0) even_timestep=false;
898  else even_timestep=true;
899 
900  const int post_offset = mmaximum_number_of_particles*int(even_timestep); //and we also save the offset to know the location in which we will save the pointers after we've moved the particles
901  //KRATOS_WATCH(post_offset)
902 
903 
904  double delta_t = CurrentProcessInfo[DELTA_TIME];
905 
906  const array_1d<double,3> gravity= CurrentProcessInfo[GRAVITY];
907 
908  array_1d<double,TDim+1> N;
909  const unsigned int max_results = 10000;
910 
911  //double integration_distance= 2.0;
912 
913  max_nsubsteps = 10;
914  max_substep_dt=delta_t/double(max_nsubsteps);
915 
916  vector<unsigned int> element_partition;
917  #ifdef _OPENMP
918  int number_of_threads = omp_get_max_threads();
919  #else
920  int number_of_threads = 1;
921  #endif
922  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
923 
924  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
925 
926 
927 
928  //before doing anything we must reset the vector of nodes contained by each element (particles that are inside each element.
929  #pragma omp parallel for
930  for(int kkk=0; kkk<number_of_threads; kkk++)
931  {
932  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
933  {
934  ModelPart::ElementsContainerType::iterator old_element = ielembegin+ii;
935 
936  int & number_of_particles = old_element->GetValue(NUMBER_OF_FLUID_PARTICLES);
937 
938  mnumber_of_particles_in_elems_aux(ii)=number_of_particles;
939  mnumber_of_particles_in_elems(ii)=0;
940  //we reset the local vectors for a faster access;
941  }
942  }
943 
944 
945  bool nonzero_mesh_velocity = false;
946  //seeing if we have to use the mesh_velocity or not
947  for(ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
948  inode!=mr_model_part.NodesEnd(); inode++)
949  {
950  const array_1d<double, 3 > velocity = inode->FastGetSolutionStepValue(MESH_VELOCITY);
951  for(unsigned int i = 0; i!=3; i++)
952  {
953  if (fabs(velocity[i])>1.0e-9)
954  nonzero_mesh_velocity=true;
955  }
956  if( nonzero_mesh_velocity==true)
957  break;
958  }
959 
960  if ( nonzero_mesh_velocity==true)
961  muse_mesh_velocity_to_convect = true; // if there is mesh velocity, then we have to take it into account when moving the particles
962  else
963  muse_mesh_velocity_to_convect = false; //otherwise, we can avoid reading the values since we know it is zero everywhere (to save time!)
964 
965  KRATOS_INFO("MoveParticleUtilityPfem2") << "Convecting particles" << std::endl;
966  //We move the particles across the fixed mesh and saving change data into them (using the function MoveParticle)
967 
968  const bool local_use_mesh_velocity_to_convect = muse_mesh_velocity_to_convect;
969 
970  #pragma omp parallel for
971  for(int kkk=0; kkk<number_of_threads; kkk++)
972  {
973 
974  const array_1d<double,3> mesh_displacement = mcalculation_domain_added_displacement; //if it is a standard problem, displacements are zero and therefore nothing is added.
975  ResultContainerType results(max_results);
976 
977  GlobalPointersVector< Element > elements_in_trajectory;
978  elements_in_trajectory.resize(20);
979 
980  for(unsigned int ielem=element_partition[kkk]; ielem<element_partition[kkk+1]; ielem++)
981  {
982  //for(unsigned int ielem=0; ielem<mr_model_part.Elements().size(); ielem++)
983  //{
984 
985  ModelPart::ElementsContainerType::iterator old_element = ielembegin+ielem;
986  const int old_element_id = old_element->Id();
987 
988  ParticlePointerVector& old_element_particle_pointers = *mpointers_to_particle_pointers_vectors(old_element_id-1);
989 
990  if ( (results.size()) !=max_results)
991  results.resize(max_results);
992 
993 
994 
995  unsigned int number_of_elements_in_trajectory=0; //excluding the origin one (current one, ielem)
996 
997  for(int ii=0; ii<(mnumber_of_particles_in_elems_aux(ielem)); ii++)
998  {
999 
1000  PFEM_Particle_Fluid & pparticle = old_element_particle_pointers[offset+ii];
1001 
1002  Element::Pointer pcurrent_element( *old_element.base() );
1003  ResultIteratorType result_begin = results.begin();
1004  bool & erase_flag=pparticle.GetEraseFlag();
1005  if (erase_flag==false){
1006  MoveParticle(pparticle,pcurrent_element,elements_in_trajectory,number_of_elements_in_trajectory,result_begin,max_results, mesh_displacement, discriminate_streamlines, local_use_mesh_velocity_to_convect); //saqué N de los argumentos, no lo necesito ya q empieza SIEMPRE en un nodo y no me importa donde termina
1007 
1008  const int current_element_id = pcurrent_element->Id();
1009 
1010  int & number_of_particles_in_current_elem = mnumber_of_particles_in_elems(current_element_id-1);
1011  //int & number_of_water_particles_in_current_elem = mnumber_of_water_particles_in_elems(current_element_id-1);
1012 
1013  if (number_of_particles_in_current_elem<mmaximum_number_of_particles && erase_flag==false)
1014  {
1015  {
1016 
1017  ParticlePointerVector& current_element_particle_pointers = *mpointers_to_particle_pointers_vectors(current_element_id-1);
1018 
1019  #pragma omp critical
1020  {
1021  if (number_of_particles_in_current_elem<mmaximum_number_of_particles) // we cant go over this node, there's no room. otherwise we would be in the position of the first particle of the next element!!
1022  {
1023 
1024  current_element_particle_pointers(post_offset+number_of_particles_in_current_elem) = &pparticle;
1025 
1026  number_of_particles_in_current_elem++ ;
1027  if (number_of_particles_in_current_elem>mmaximum_number_of_particles)
1028  KRATOS_WATCH("MAL");
1029 
1030  }
1031  else
1032  pparticle.GetEraseFlag()=true; //so we just delete it!
1033  }
1034  }
1035  }
1036  else
1037  pparticle.GetEraseFlag()=true; //so we just delete it!
1038 
1039 
1040 
1041  }
1042  }
1043  }
1044  }
1045 
1046  //now we pass info from the local vector to the elements:
1047  #pragma omp parallel for
1048  for(int kkk=0; kkk<number_of_threads; kkk++)
1049  {
1050  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
1051  {
1052  ModelPart::ElementsContainerType::iterator old_element = ielembegin+ii;
1053 
1054  old_element->GetValue(NUMBER_OF_FLUID_PARTICLES) = mnumber_of_particles_in_elems(ii);
1055  //old_element->GetValue(NUMBER_OF_WATER_PARTICLES) = mnumber_of_water_particles_in_elems(ii);
1056  }
1057 
1058  }
1059 
1060 
1061  //after having changed everything we change the status of the modd_timestep flag:
1062  CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET] = post_offset;; //
1063 
1064  KRATOS_CATCH("")
1065  }
1066 
1067  void TransferLagrangianToEulerian() //explicit
1068  {
1069  KRATOS_TRY
1070 
1071  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
1072  //const double delta_t =CurrentProcessInfo[DELTA_TIME];
1073  const double threshold= 0.0/(double(TDim)+1.0);
1074 
1075 
1076  KRATOS_INFO("MoveParticleUtilityPfem2") << "Projecting info to mesh" << std::endl;
1077 
1078 
1079  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
1080  //KRATOS_WATCH(offset) //(flag managed only by MoveParticles
1081 
1082  //we must project data from the particles (lagrangian) into the eulerian mesh
1083  //ValuesVectorType eulerian_nodes_old_temperature;
1084  //int nnodes = mr_model_part.Nodes().size();
1085  //array_1d<double,(n_nodes)> eulerian_nodes_sumweights;
1086 
1087  //we save data from previous time step of the eulerian mesh in case we must reuse it later cos no particle was found around the nodes
1088  //though we could've use a bigger buffer, to be changed later!
1089  //after having saved data, we reset them to zero, this way it's easier to add the contribution of the surrounding particles.
1090  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
1091  vector<unsigned int> node_partition;
1092  #ifdef _OPENMP
1093  int number_of_threads = omp_get_max_threads();
1094  #else
1095  int number_of_threads = 1;
1096  #endif
1097  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
1098 
1099  #pragma omp parallel for
1100  for(int kkk=0; kkk<number_of_threads; kkk++)
1101  {
1102  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
1103  {
1104  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
1105  inode->FastGetSolutionStepValue(DISTANCE)=0.0;
1106  inode->FastGetSolutionStepValue(PROJECTED_VELOCITY)=ZeroVector(3);
1107  inode->FastGetSolutionStepValue(YP)=0.0;
1108  }
1109 
1110  }
1111 
1112  //adding contribution, loop on elements, since each element has stored the particles found inside of it
1113  vector<unsigned int> element_partition;
1114  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
1115 
1116  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
1117  #pragma omp parallel for
1118  for(int kkk=0; kkk<number_of_threads; kkk++)
1119  {
1120  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
1121  {
1122  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
1123 
1124  array_1d<double,3*(TDim+1)> nodes_positions;
1125  array_1d<double,3*(TDim+1)> nodes_addedvel = ZeroVector(3*(TDim+1));
1126 
1127  array_1d<double,(TDim+1)> nodes_added_distance = ZeroVector((TDim+1));
1128  array_1d<double,(TDim+1)> nodes_addedweights = ZeroVector((TDim+1));
1129  //array_1d<double,(TDim+1)> weighting_inverse_divisor;
1130 
1131  Geometry<Node >& geom = ielem->GetGeometry();
1132 
1133  for (int i=0 ; i!=(TDim+1) ; ++i)
1134  {
1135  nodes_positions[i*3+0]=geom[i].X();
1136  nodes_positions[i*3+1]=geom[i].Y();
1137  nodes_positions[i*3+2]=geom[i].Z();
1138  //weighting_inverse_divisor[i]=1.0/((geom[i].FastGetSolutionStepValue(MEAN_SIZE))*1.01);
1139  }
1142 
1143  int & number_of_particles_in_elem= ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
1144  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
1145 
1146  for (int iii=0; iii<number_of_particles_in_elem ; iii++ )
1147  {
1148  if (iii==mmaximum_number_of_particles) //it means we are out of our portion of the array, abort loop!
1149  break;
1150 
1151  PFEM_Particle_Fluid & pparticle = element_particle_pointers[offset+iii];
1152 
1153  if (pparticle.GetEraseFlag()==false)
1154  {
1155 
1156  array_1d<double,3> & position = pparticle.Coordinates();
1157 
1158  const array_1d<float,3>& velocity = pparticle.GetVelocity();
1159 
1160  const float& particle_distance = pparticle.GetDistance(); // -1 if water, +1 if air
1161 
1162  array_1d<double,TDim+1> N;
1163  bool is_found = CalculatePosition(nodes_positions,position[0],position[1],position[2],N);
1164  if (is_found==false) //something went wrong. if it was close enough to the edge we simply send it inside the element.
1165  {
1166  KRATOS_WATCH(N);
1167  for (int j=0 ; j!=(TDim+1); j++)
1168  if (N[j]<0.0 && N[j]> -1e-5)
1169  N[j]=1e-10;
1170 
1171  }
1172 
1173  for (int j=0 ; j!=(TDim+1); j++) //going through the 3/4 nodes of the element
1174  {
1175  //double sq_dist = 0;
1176  //these lines for a weighting function based on the distance (or square distance) from the node insteadof the shape functions
1177  //for (int k=0 ; k!=(TDim); k++) sq_dist += ((position[k] - nodes_positions[j*3+k])*(position[k] - nodes_positions[j*3+k]));
1178  //double weight = (1.0 - (sqrt(sq_dist)*weighting_inverse_divisor[j] ) );
1179 
1180  double weight=N(j);
1181  //weight=N(j)*N(j)*N(j);
1182  if (weight<threshold) weight=1e-10;
1183  if (weight<0.0) {KRATOS_WATCH(weight)}//;weight=0.0;KRATOS_WATCH(velocity);KRATOS_WATCH(N);KRATOS_WATCH(number_of_particles_in_elem);}//{KRATOS_WATCH(weight); KRATOS_WATCH(geom[j].Id()); KRATOS_WATCH(position);}
1184  else
1185  {
1186  nodes_addedweights[j]+= weight;
1187  //nodes_addedtemp[j] += weight * particle_temp;
1188 
1189  nodes_added_distance[j] += weight*particle_distance;
1190 
1191  //nodes_added_oxygen[j] += weight*particle_oxygen;
1192 
1193  for (int k=0 ; k!=(TDim); k++) //x,y,(z)
1194  {
1195  nodes_addedvel[j*3+k] += weight * double(velocity[k]);
1196  }
1197 
1198  }//
1199  }
1200  }
1201  }
1202 
1203  for (int i=0 ; i!=(TDim+1) ; ++i) {
1204  geom[i].SetLock();
1205  geom[i].FastGetSolutionStepValue(DISTANCE) +=nodes_added_distance[i];
1206  geom[i].FastGetSolutionStepValue(PROJECTED_VELOCITY_X) +=nodes_addedvel[3*i+0];
1207  geom[i].FastGetSolutionStepValue(PROJECTED_VELOCITY_Y) +=nodes_addedvel[3*i+1];
1208  geom[i].FastGetSolutionStepValue(PROJECTED_VELOCITY_Z) +=nodes_addedvel[3*i+2]; //we are updating info to the previous time step!!
1209 
1210  geom[i].FastGetSolutionStepValue(YP) +=nodes_addedweights[i];
1211  geom[i].UnSetLock();
1212  }
1213  }
1214  }
1215 
1216  #pragma omp parallel for
1217  for(int kkk=0; kkk<number_of_threads; kkk++)
1218  {
1219  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
1220  {
1221  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
1222  double sum_weights = inode->FastGetSolutionStepValue(YP);
1223  if (sum_weights>0.00001)
1224  {
1225  //inode->FastGetSolutionStepValue(TEMPERATURE_OLD_IT)=(inode->FastGetSolutionStepValue(TEMPERATURE_OLD_IT))/sum_weights; //resetting the temperature
1226  double & dist = inode->FastGetSolutionStepValue(DISTANCE);
1227  dist /=sum_weights; //resetting the density
1228  inode->FastGetSolutionStepValue(PROJECTED_VELOCITY)=(inode->FastGetSolutionStepValue(PROJECTED_VELOCITY))/sum_weights; //resetting the velocity
1229 
1230  }
1231 
1232  else //this should never happen because other ways to recover the information have been executed before, but leaving it just in case..
1233  {
1234  inode->FastGetSolutionStepValue(DISTANCE)=3.0; //resetting the temperature
1235  //inode->FastGetSolutionStepValue(DISTANCE)=inode->GetSolutionStepValue(DISTANCE,1); //resetting the temperature
1236  inode->FastGetSolutionStepValue(PROJECTED_VELOCITY)=inode->GetSolutionStepValue(VELOCITY,1);
1237 
1238  }
1240  if (inode->IsFixed(DISTANCE))
1241  inode->FastGetSolutionStepValue(DISTANCE)=inode->GetSolutionStepValue(DISTANCE,1);
1242  }
1243  }
1244 
1245 
1246  KRATOS_CATCH("")
1247  }
1248 
1249 
1250 
1251  void TransferLagrangianToEulerianImp() //semi implicit
1252  {
1253  KRATOS_TRY
1254 
1255  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
1256 
1257  std::cout << "projecting info to mesh (semi implicit)" << std::endl;
1258 
1259 
1260  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
1261  //KRATOS_WATCH(offset) //(flag managed only by MoveParticles
1262 
1263  //we must project data from the particles (lagrangian) into the eulerian mesh
1264  //ValuesVectorType eulerian_nodes_old_temperature;
1265  //int nnodes = mr_model_part.Nodes().size();
1266  //array_1d<double,(n_nodes)> eulerian_nodes_sumweights;
1267 
1268  //we save data from previous time step of the eulerian mesh in case we must reuse it later cos no particle was found around the nodes
1269  //though we could've use a bigger buffer, to be changed later!
1270  //after having saved data, we reset them to zero, this way it's easier to add the contribution of the surrounding particles.
1271  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
1272  vector<unsigned int> node_partition;
1273  #ifdef _OPENMP
1274  int number_of_threads = omp_get_max_threads();
1275  #else
1276  int number_of_threads = 1;
1277  #endif
1278  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
1279 
1280  #pragma omp parallel for
1281  for(int kkk=0; kkk<number_of_threads; kkk++)
1282  {
1283  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
1284  {
1285  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
1286  inode->FastGetSolutionStepValue(DISTANCE)=0.0;
1287  inode->FastGetSolutionStepValue(PROJECTED_VELOCITY)=ZeroVector(3);
1288  inode->FastGetSolutionStepValue(YP)=0.0;
1289  }
1290 
1291  }
1292 
1293  //adding contribution, loop on elements, since each element has stored the particles found inside of it
1294  vector<unsigned int> element_partition;
1295  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
1296 
1297  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
1298  #pragma omp parallel for
1299  for(int kkk=0; kkk<number_of_threads; kkk++)
1300  {
1301 
1302  //creating a matrix for each of the problems.
1303  BoundedMatrix<double, TDim+1 , TDim+1 > mass_matrix; // WE ONLY NEED ONE! they are the same for all the variables! //_x,mass_matrix_y,mass_matrix_z,mass_matrix_d; //mass matrices for the projected vel (x,y,z) and the distance
1304  array_1d<double,(TDim+1)> rhs_x,rhs_y,rhs_z,rhs_d;
1305 
1306  array_1d<double,3*(TDim+1)> nodes_positions;
1307  array_1d<double,3*(TDim+1)> nodes_addedvel = ZeroVector(3*(TDim+1));
1308 
1309  array_1d<double,(TDim+1)> nodes_added_distance = ZeroVector((TDim+1));
1310  array_1d<double,(TDim+1)> nodes_addedweights = ZeroVector((TDim+1));
1311 
1312  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
1313  {
1314  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
1315 
1316  nodes_addedvel = ZeroVector(3*(TDim+1)); //resetting vectors
1317  nodes_added_distance = ZeroVector((TDim+1)); //resetting vectors
1318  nodes_addedweights = ZeroVector((TDim+1)); //resetting vectors
1319  mass_matrix = ZeroMatrix(TDim+1 , TDim+1 ); //resetting matrices. WE ONLY NEED ONE! they are the same for all the variable. only the rhs changes.
1320  //mass_matrix_y = ZeroMatrix(TDim+1 , TDim+1 ); //resetting matrices
1321  //mass_matrix_z = ZeroMatrix(TDim+1 , TDim+1 ); //resetting matrices
1322  //mass_matrix_d = ZeroMatrix(TDim+1 , TDim+1 ); //resetting matrices
1323  rhs_x = ZeroVector((TDim+1)); //resetting vectors
1324  rhs_y = ZeroVector((TDim+1)); //resetting vectors
1325  rhs_z = ZeroVector((TDim+1)); //resetting vectors
1326  rhs_d = ZeroVector((TDim+1)); //resetting vectors
1327 
1328  Geometry<Node >& geom = ielem->GetGeometry();
1329  const double elem_volume = geom.Area();
1330 
1331  for (int i=0 ; i!=(TDim+1) ; ++i) //saving the nodal positions for faster access
1332  {
1333  nodes_positions[i*3+0]=geom[i].X();
1334  nodes_positions[i*3+1]=geom[i].Y();
1335  nodes_positions[i*3+2]=geom[i].Z();
1336  }
1339 
1340  int & number_of_particles_in_elem= ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
1341  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
1342 
1343  for (int iii=0; iii<number_of_particles_in_elem ; iii++ )
1344  {
1345  if (iii==mmaximum_number_of_particles) //it means we are out of our portion of the array, abort loop!
1346  break;
1347 
1348  PFEM_Particle_Fluid & pparticle = element_particle_pointers[offset+iii];
1349 
1350  if (pparticle.GetEraseFlag()==false)
1351  {
1352 
1353  array_1d<double,3> & position = pparticle.Coordinates();
1354 
1355  const array_1d<float,3>& velocity = pparticle.GetVelocity();
1356 
1357  const float& particle_distance = pparticle.GetDistance(); // -1 if water, +1 if air
1358 
1359  array_1d<double,TDim+1> N;
1360  bool is_found = CalculatePosition(nodes_positions,position[0],position[1],position[2],N);
1361  if (is_found==false) //something went wrong. if it was close enough to the edge we simply send it inside the element.
1362  {
1363  KRATOS_WATCH(N);
1364  for (int j=0 ; j!=(TDim+1); j++)
1365  if (N[j]<0.0 && N[j]> -1e-5)
1366  N[j]=1e-10;
1367 
1368  }
1369 
1370  for (int j=0 ; j!=(TDim+1); j++) //going through the 3/4 nodes of the element
1371  {
1372  double weight=N(j);
1373  for (int k=0 ; k!=(TDim+1); k++) //building the mass matrix
1374  mass_matrix(j,k) += weight*N(k);
1375 
1376  rhs_x[j] += weight * double(velocity[0]);
1377  rhs_y[j] += weight * double(velocity[1]);
1378  rhs_z[j] += weight * double(velocity[2]);
1379  rhs_d[j] += weight * double(particle_distance);
1380 
1381  //adding also a part with the lumped mass matrix to reduce overshoots and undershoots
1382  if(true)
1383  {
1384  double this_particle_weight = weight*elem_volume/(double(number_of_particles_in_elem))*0.1; //can be increased or reduced to change the lumped mass contrubtion
1385  nodes_addedweights[j]+= this_particle_weight;
1386  nodes_added_distance[j] += this_particle_weight*particle_distance;
1387  for (int k=0 ; k!=(TDim); k++) //x,y,(z)
1388  {
1389  nodes_addedvel[j*3+k] += this_particle_weight * double(velocity[k]);
1390  }
1391  }
1392  }
1393  }
1394  }
1395 
1396  //now we invert the matrix
1397  BoundedMatrix<double, TDim+1 , TDim+1 > inverse_mass_matrix=ZeroMatrix(TDim+1 , TDim+1);
1398  if constexpr (TDim==3)
1399  InvertMatrix( mass_matrix, inverse_mass_matrix);
1400  else
1401  InvertMatrix3x3( mass_matrix, inverse_mass_matrix);
1402  //and now compute the elemental contribution to the gobal system:
1403 
1404  if(number_of_particles_in_elem>(TDim*3)) //otherwise it's impossible to define a correctly the gradients, therefore the results inside the element are useless.
1405  {
1406  for (int i=0 ; i!=(TDim+1); i++)
1407  {
1408  for (int j=0 ; j!=(TDim+1); j++)
1409  {
1410  nodes_addedvel[3*i+0] += inverse_mass_matrix(i,j)*rhs_x[j]*elem_volume*(1.0/(double(1+TDim)));
1411  nodes_addedvel[3*i+1] += inverse_mass_matrix(i,j)*rhs_y[j]*elem_volume*(1.0/(double(1+TDim)));
1412  nodes_addedvel[3*i+2] += inverse_mass_matrix(i,j)*rhs_z[j]*elem_volume*(1.0/(double(1+TDim)));
1413  nodes_added_distance[i] += inverse_mass_matrix(i,j)*rhs_d[j]*elem_volume*(1.0/(double(1+TDim)));
1414 
1415  }
1416  }
1417  //and also to the mass matrix. LUMPED (but for the contribution of the grandient at elemental level.
1418  for (int i=0 ; i!=(TDim+1); i++)
1419  nodes_addedweights[i] += elem_volume*(1.0/(double(1+TDim)));
1420  }
1421 
1422 
1423  for (int i=0 ; i!=(TDim+1) ; ++i) {
1424  geom[i].SetLock();
1425  geom[i].FastGetSolutionStepValue(DISTANCE) +=nodes_added_distance[i];
1426  geom[i].FastGetSolutionStepValue(PROJECTED_VELOCITY_X) +=nodes_addedvel[3*i+0];
1427  geom[i].FastGetSolutionStepValue(PROJECTED_VELOCITY_Y) +=nodes_addedvel[3*i+1];
1428  geom[i].FastGetSolutionStepValue(PROJECTED_VELOCITY_Z) +=nodes_addedvel[3*i+2]; //we are updating info to the previous time step!!
1429 
1430  geom[i].FastGetSolutionStepValue(YP) +=nodes_addedweights[i];
1431  geom[i].UnSetLock();
1432  }
1433  }
1434  }
1435 
1436  #pragma omp parallel for
1437  for(int kkk=0; kkk<number_of_threads; kkk++)
1438  {
1439  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
1440  {
1441  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
1442  double sum_weights = inode->FastGetSolutionStepValue(YP);
1443  if (sum_weights>0.00001)
1444  {
1445  //inode->FastGetSolutionStepValue(TEMPERATURE_OLD_IT)=(inode->FastGetSolutionStepValue(TEMPERATURE_OLD_IT))/sum_weights; //resetting the temperature
1446  double & dist = inode->FastGetSolutionStepValue(DISTANCE);
1447  dist /=sum_weights; //resetting the density
1448  inode->FastGetSolutionStepValue(PROJECTED_VELOCITY)=(inode->FastGetSolutionStepValue(PROJECTED_VELOCITY))/sum_weights; //resetting the velocity
1449 
1450  }
1451 
1452  else //this should never happen because other ways to recover the information have been executed before, but leaving it just in case..
1453  {
1454  inode->FastGetSolutionStepValue(DISTANCE)=3.0; //resetting the temperature
1455  //inode->FastGetSolutionStepValue(DISTANCE)=inode->GetSolutionStepValue(DISTANCE,1); //resetting the temperature
1456  inode->FastGetSolutionStepValue(PROJECTED_VELOCITY)=inode->GetSolutionStepValue(VELOCITY,1);
1457 
1458  }
1460  if (inode->IsFixed(DISTANCE))
1461  inode->FastGetSolutionStepValue(DISTANCE)=inode->GetSolutionStepValue(DISTANCE,1);
1462  }
1463  }
1464 
1465 
1466  KRATOS_CATCH("")
1467  }
1468 
1469  void AccelerateParticlesWithoutMovingUsingDeltaVelocity()
1470  {
1471  KRATOS_TRY
1472  //std::cout << "updating particles" << std::endl;
1473  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
1474 
1475  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
1476  //(flag managed only by MoveParticles
1477  //KRATOS_WATCH(offset)
1478  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
1479 
1480 
1481  vector<unsigned int> element_partition;
1482  #ifdef _OPENMP
1483  int number_of_threads = omp_get_max_threads();
1484  #else
1485  int number_of_threads = 1;
1486  #endif
1487  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
1488 
1489  #pragma omp parallel for
1490  for(int kkk=0; kkk<number_of_threads; kkk++)
1491  {
1492  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
1493  {
1494  //const int & elem_id = ielem->Id();
1495  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
1496  Element::Pointer pelement(*ielem.base());
1497  Geometry<Node >& geom = ielem->GetGeometry();
1498 
1499  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
1500  int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
1501  //std::cout << "elem " << ii << " with " << (unsigned int)number_of_particles_in_elem << " particles" << std::endl;
1502 
1503  for (int iii=0; iii<number_of_particles_in_elem ; iii++ )
1504  {
1505  //KRATOS_WATCH(iii)
1506  if (iii>mmaximum_number_of_particles) //it means we are out of our portion of the array, abort loop!
1507  break;
1508 
1509  PFEM_Particle_Fluid & pparticle = element_particle_pointers[offset+iii];
1510 
1511 
1512  bool erase_flag= pparticle.GetEraseFlag();
1513  if (erase_flag==false)
1514  {
1515  AccelerateParticleUsingDeltaVelocity(pparticle,pelement,geom); //'lite' version, we pass by reference the geometry, so much cheaper
1516  }
1517 
1518 
1519  }
1520  }
1521  }
1522  KRATOS_CATCH("")
1523  }
1524 
1525  //**************************************************************************************************************
1526  //**************************************************************************************************************
1527 
1528 
1529  template< class TDataType > void AddUniqueWeakPointer
1530  (GlobalPointersVector< TDataType >& v, const typename TDataType::WeakPointer candidate)
1531  {
1532  typename GlobalPointersVector< TDataType >::iterator i = v.begin();
1533  typename GlobalPointersVector< TDataType >::iterator endit = v.end();
1534  while ( i != endit && (i)->Id() != (candidate.lock())->Id())
1535  {
1536  i++;
1537  }
1538  if( i == endit )
1539  {
1540  v.push_back(candidate);
1541  }
1542 
1543  }
1544 
1545  //**************************************************************************************************************
1546  //**************************************************************************************************************
1547 
1548  void PreReseed(int minimum_number_of_particles)
1549  {
1550  KRATOS_TRY
1551 
1552 
1553  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
1554  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET];
1555  const int max_results = 1000;
1556 
1557  //tools for the paralelization
1558  unsigned int number_of_threads = OpenMPUtils::GetNumThreads();
1559  vector<unsigned int> elem_partition;
1560  int number_of_rows=mr_model_part.Elements().size();
1561  elem_partition.resize(number_of_threads + 1);
1562  int elem_partition_size = number_of_rows / number_of_threads;
1563  elem_partition[0] = 0;
1564  elem_partition[number_of_threads] = number_of_rows;
1565  //KRATOS_WATCH(elem_partition_size);
1566  for (unsigned int i = 1; i < number_of_threads; i++)
1567  elem_partition[i] = elem_partition[i - 1] + elem_partition_size;
1568 
1569  const bool local_use_mesh_velocity_to_convect = muse_mesh_velocity_to_convect;
1570  #pragma omp parallel firstprivate(elem_partition)
1571  {
1572  ResultContainerType results(max_results);
1573  int k = OpenMPUtils::ThisThread();
1574  ModelPart::ElementsContainerType::iterator it_begin = mr_model_part.ElementsBegin() + elem_partition[k];
1575  ModelPart::ElementsContainerType::iterator it_end = mr_model_part.ElementsBegin() + elem_partition[k+1] ;
1576  //ModelPart::NodesContainerType local_list=aux[k];
1577  //PointerVectorSet<PFEM_Particle_Fluid, IndexedObject> & list=aux[k];
1578  //KRATOS_WATCH(k);
1579  BoundedMatrix<double, (TDim+1), 3 > pos;
1580  BoundedMatrix<double, (TDim+1) , (TDim+1) > N;
1581  unsigned int freeparticle=0; //we start with the first position in the particles array
1582 
1583  //int local_id=1;
1584  for (ModelPart::ElementsContainerType::iterator ielem = it_begin; ielem != it_end; ielem++)
1585  {
1586  results.resize(max_results);
1587  //const int & elem_id = ielem->Id();
1588  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
1589  int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
1590  if (number_of_particles_in_elem<(minimum_number_of_particles))// && (ielem->GetGeometry())[0].Y()<0.10 )
1591  {
1592  //KRATOS_WATCH("elem with little particles")
1593  Geometry< Node >& geom = ielem->GetGeometry();
1594  ComputeGaussPointPositionsForPreReseed(geom, pos, N);
1595  //double conductivity = ielem->GetProperties()[CONDUCTIVITY];
1596  //KRATOS_WATCH(conductivity);
1597  for (unsigned int j = 0; j < (pos.size1()); j++) //i am dropping the last one, the one in the middle of the element
1598  {
1599  bool keep_looking = true;
1600  while(keep_looking)
1601  {
1602  if (mparticles_vector[freeparticle].GetEraseFlag()==true)
1603  {
1604  #pragma omp critical
1605  {
1606  if (mparticles_vector[freeparticle].GetEraseFlag()==true)
1607  {
1608  mparticles_vector[freeparticle].GetEraseFlag()=false;
1609  keep_looking=false;
1610  }
1611  }
1612  if (keep_looking==false)
1613  break;
1614  else
1615  freeparticle++;
1616  }
1617  else
1618  {
1619  freeparticle++;
1620  }
1621  }
1622 
1623  PFEM_Particle_Fluid pparticle(pos(j,0),pos(j,1),pos(j,2));
1624 
1625  array_1d<double,TDim+1>aux2_N;
1626  bool is_found = CalculatePosition(geom,pos(j,0),pos(j,1),pos(j,2),aux2_N);
1627  if (is_found==false)
1628  {
1629  KRATOS_WATCH(aux2_N);
1630  }
1631 
1632 
1633  pparticle.GetEraseFlag()=false;
1634 
1635  ResultIteratorType result_begin = results.begin();
1636  Element::Pointer pelement( *ielem.base() );
1637  MoveParticle_inverse_way(pparticle, pelement, result_begin, max_results, local_use_mesh_velocity_to_convect);
1638 
1639  //and we copy it to the array:
1640  mparticles_vector[freeparticle] = pparticle;
1641 
1642  element_particle_pointers(offset+number_of_particles_in_elem) = &mparticles_vector[freeparticle];
1643  pparticle.GetEraseFlag()=false;
1644 
1645  number_of_particles_in_elem++;
1646 
1647 
1648  }
1649  }
1650  }
1651  }
1652 
1653 
1654 
1655 
1656  KRATOS_CATCH("")
1657  }
1658 
1659 
1660  //**************************************************************************************************************
1661  //**************************************************************************************************************
1662 
1663 
1664  void PostReseed(int minimum_number_of_particles, double mass_correction_factor ) //pooyan's way
1665  {
1666  KRATOS_TRY
1667 
1668  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
1669  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET];
1670 
1671  if (mass_correction_factor>0.5) mass_correction_factor=0.5;
1672  if (mass_correction_factor<-0.5) mass_correction_factor=-0.5;
1673  //mass_correction_factor=0.0;
1674 
1675  //ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
1676  //const double delta_t = CurrentProcessInfo[DELTA_TIME];
1677  //array_1d<double,3> & gravity= CurrentProcessInfo[GRAVITY];
1678  //const int max_results = 1000;
1679 
1680  const double threshold = mass_correction_factor*0.5;
1681 
1682  //TOOLS FOR THE PARALELIZATION
1683  //int last_id= (mr_linea_model_part.NodesEnd()-1)->Id();
1684  unsigned int number_of_threads = OpenMPUtils::GetNumThreads();
1685  //KRATOS_WATCH(number_of_threads);
1686  vector<unsigned int> elem_partition;
1687  int number_of_rows=mr_model_part.Elements().size();
1688  //KRATOS_WATCH(number_of_threads);
1689  //KRATOS_THROW_ERROR(std::logic_error, "Add ----NODAL_H---- variable!!!!!! ERROR", "");
1690  elem_partition.resize(number_of_threads + 1);
1691  int elem_partition_size = number_of_rows / number_of_threads;
1692  elem_partition[0] = 0;
1693  elem_partition[number_of_threads] = number_of_rows;
1694  //KRATOS_WATCH(elem_partition_size);
1695  for (unsigned int i = 1; i < number_of_threads; i++)
1696  elem_partition[i] = elem_partition[i - 1] + elem_partition_size;
1697  //typedef Node PointType;
1698  //std::vector<ModelPart::NodesContainerType> aux;// aux;
1699  //aux.resize(number_of_threads);
1700 
1701  //ModelPart::NodesContainerType::iterator it_begin_particle_model_part = mr_linea_model_part.NodesBegin();
1702  //ModelPart::NodesContainerType::iterator it_end_particle_model_part = mr_linea_model_part.NodesEnd();
1703 
1704  #pragma omp parallel firstprivate(elem_partition) // firstprivate(results)//we will add the nodes in different parts of aux and later assemple everything toghether, remaming particles ids to get consecutive ids
1705  {
1706  unsigned int reused_particles=0;
1707 
1708  unsigned int freeparticle = 0; //we start by the first position;
1709 
1710  int k = OpenMPUtils::ThisThread();
1711  ModelPart::ElementsContainerType::iterator it_begin = mr_model_part.ElementsBegin() + elem_partition[k];
1712  ModelPart::ElementsContainerType::iterator it_end = mr_model_part.ElementsBegin() + elem_partition[k+1] ;
1713 
1714  BoundedMatrix<double, (3+2*TDim), 3 > pos; //7 particles (2D) or 9 particles (3D)
1715  BoundedMatrix<double, (3+2*TDim), (TDim+1) > N;
1716 
1717  array_1d<double, 3 > vel_complete, vel_without_air_nodes;
1718  double sum_Ns_without_air_nodes;
1719  double mesh_distance;
1720 
1721  array_1d<double, (3+2*TDim) > distances;
1722  array_1d<int, (3+2*TDim) > positions;
1723  array_1d<bool, (3+2*TDim) > is_water_particle; //for both
1724 
1725 
1726  unsigned int number_of_reseeded_particles;
1727  //unsigned int number_of_water_reseeded_particles;
1728 
1729  //array_1d<double, 3 > nodes_distances;
1730 
1731  //int local_id=1;
1732  for (ModelPart::ElementsContainerType::iterator ielem = it_begin; ielem != it_end; ielem++)
1733  {
1734  //results.resize(max_results);
1735 
1736  int & number_of_particles_in_elem= ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
1737  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
1738 
1739 
1740  Geometry< Node >& geom = ielem->GetGeometry();
1741  if ( (number_of_particles_in_elem<(minimum_number_of_particles)))// && (geom[0].Y()<0.10) ) || (number_of_water_particles_in_elem>2 && number_of_particles_in_elem<(minimum_number_of_particles) ) )
1742  {
1743 
1744  //bool reseed_more=false;
1745  number_of_reseeded_particles=0;
1746 
1747 
1748  //reseed_more=true;
1749  number_of_reseeded_particles= 3+2*TDim;
1750  ComputeGaussPointPositionsForPostReseed(geom, pos, N);
1751 
1752  distances = ZeroVector(3+2*TDim);
1753 
1754  bool has_water_node=false;
1755  bool has_air_node=false;
1756  double mean_element_distance = 0.0;
1757 
1758 
1759  for (unsigned int j = 0; j < (TDim+1); j++)
1760  {
1761  mean_element_distance += (1.0/double(TDim+1))*(geom[j].FastGetSolutionStepValue(DISTANCE));
1762  if ((geom[j].FastGetSolutionStepValue(DISTANCE))<0.0)
1763  has_water_node=true;
1764  else
1765  has_air_node=true;
1766  }
1767 
1768  //first we check the particle distance according to the nodal values
1769  for (unsigned int j = 0; j < number_of_reseeded_particles; j++) //first we order particles
1770  {
1771  positions[j]=j+1; //just creating a vector from 1 to 7 or whathever our lenght is (7 for 2d, 9 for 3d)
1772  for (unsigned int l = 0; l < (TDim+1); l++)
1773  {
1774  distances[j] += N(j, l) * geom[l].FastGetSolutionStepValue(DISTANCE);
1775  }
1776  }
1777 
1778 
1779  if ( (has_air_node && has_water_node) ) //for slit elements we use the distance function
1780  {
1781  for (unsigned int j = 0; j < number_of_reseeded_particles ; j++) //first we order particles
1782  {
1783  if (distances[j]>threshold)
1784  is_water_particle[j]=false;
1785  else
1786  is_water_particle[j]=true;
1787  }
1788  }
1789  else if (has_air_node)
1790  {
1791  double water_fraction = 0.5 - 0.5*(mean_element_distance);
1792  if (water_fraction>0.9 && mass_correction_factor<0.0) //to avoid seeding air particles when we are in a pure water element
1793  mass_correction_factor = 0.0;
1794  unsigned int number_of_water_reseeded_particles = double(number_of_reseeded_particles)*(1.01+mass_correction_factor*1.0)*water_fraction;
1795 
1796  BubbleSort(distances, positions, number_of_reseeded_particles); //ok. now we have the particles ordered from the "watermost" to "airmost". therefore we will fill the water particles and later the air ones using that order
1797 
1798  for (unsigned int j = 0; j < number_of_reseeded_particles ; j++) //first we order particles
1799  {
1800  int array_position = positions[j]-1;
1801  if (array_position>3 && number_of_reseeded_particles==4)
1802  {
1803  KRATOS_WATCH("error in reseeding")
1804  }
1805 
1806  if ( (j+1) <= number_of_water_reseeded_particles ) //means it is a water particle
1807  is_water_particle[array_position]=true;
1808  else
1809  is_water_particle[array_position]=false;
1810  }
1811  }
1812  else //only water particles
1813  {
1814  for (unsigned int j = 0; j < number_of_reseeded_particles ; j++) //first we order particles
1815  is_water_particle[j]=true;
1816  }
1817 
1818  bool fix_distance = false;
1819  unsigned int node_with_fixed_distance = 0;
1820  for (unsigned int j = 0; j < (TDim+1) ; j++) //we go over the 3/4 nodes:
1821  {
1822  if ((geom[j].IsFixed(DISTANCE)))
1823  {
1824  fix_distance = true;
1825  node_with_fixed_distance = j;
1826  }
1827  }
1828  // so now if the 3 were fixed, we assign the sign of the first node to all the particles:
1829  if (fix_distance)
1830  {
1831  bool is_water_for_all_particles=true;
1832  if ((geom[node_with_fixed_distance].FastGetSolutionStepValue(DISTANCE))>0.0)
1833  is_water_for_all_particles=false;
1834 
1835  for (unsigned int j = 0; j < number_of_reseeded_particles ; j++) //first we order particles
1836  is_water_particle[j]=is_water_for_all_particles;
1837  }
1838 
1839 
1840  for (unsigned int j = 0; j < number_of_reseeded_particles; j++)
1841  {
1842  //now we have to find an empty space ( a particle that was about to be deleted) in the particles model part. once found. there will be our renewed particle:
1843  bool keep_looking = true;
1844  while(keep_looking)
1845  {
1846  if (mparticles_vector[freeparticle].GetEraseFlag()==true)
1847  {
1848  #pragma omp critical
1849  {
1850  if (mparticles_vector[freeparticle].GetEraseFlag()==true)
1851  {
1852  mparticles_vector[freeparticle].GetEraseFlag()=false;
1853  keep_looking=false;
1854  }
1855  }
1856  if (keep_looking==false)
1857  break;
1858 
1859  else
1860  freeparticle++;
1861  }
1862  else
1863  {
1864  freeparticle++;
1865  }
1866  }
1867 
1868 
1869  PFEM_Particle_Fluid pparticle(pos(j,0),pos(j,1),pos(j,2));
1870 
1871 
1872  array_1d<float, 3 > & vel = pparticle.GetVelocity();
1873  float& distance= pparticle.GetDistance();
1874 
1875  array_1d<double,TDim+1>aux_N;
1876  bool is_found = CalculatePosition(geom,pos(j,0),pos(j,1),pos(j,2),aux_N);
1877  if (is_found==false)
1878  {
1879  KRATOS_WATCH(aux_N);
1880  KRATOS_WATCH(j)
1881  KRATOS_WATCH(ielem->Id())
1882  }
1883 
1884 
1885  noalias(vel_complete)=ZeroVector(3);
1886  noalias(vel_without_air_nodes)=ZeroVector(3);
1887  sum_Ns_without_air_nodes=0.0;
1888 
1889  noalias(vel) = ZeroVector(3);
1890  distance=0.0;
1891  mesh_distance = 0.0;
1892  //oxygen = 0.0;
1893 
1894  for (unsigned int l = 0; l < (TDim+1); l++)
1895  {
1896  noalias(vel_complete) += N(j, l) * geom[l].FastGetSolutionStepValue(VELOCITY);
1897  mesh_distance += N(j,l) * geom[l].FastGetSolutionStepValue(DISTANCE);
1898  if ((geom[l].FastGetSolutionStepValue(DISTANCE))<0.0)
1899  {
1900  sum_Ns_without_air_nodes+=N(j, l);
1901  noalias(vel_without_air_nodes) += N(j, l) * geom[l].FastGetSolutionStepValue(VELOCITY);
1902  }
1903  }
1904 
1906  if (is_water_particle[j])
1907  {
1908  distance=-1.0;
1909  }
1910  else
1911  {
1912  //if (mesh_distance<2.0)
1913  distance=1.0;
1914  //else
1915  // distance=3.0;
1916  }
1917 
1918  if (distance<0.0 && sum_Ns_without_air_nodes>0.01)
1919  vel = vel_without_air_nodes / sum_Ns_without_air_nodes ;
1920  else
1921  vel = vel_complete;
1922 
1923  pparticle.GetEraseFlag()=false;
1924 
1925 
1926  mparticles_vector[freeparticle]=pparticle;
1927  element_particle_pointers(offset+number_of_particles_in_elem) = &mparticles_vector[freeparticle];
1928  number_of_particles_in_elem++;
1929 
1930 
1931  if (keep_looking)
1932  {
1933  KRATOS_THROW_ERROR(std::logic_error, "FINISHED THE LIST AND COULDNT FIND A FREE CELL FOR THE NEW PARTICLE!", "");
1934  }
1935  else
1936  {
1937  reused_particles++;
1938  }
1939 
1940  }
1941  }
1942  }
1943  }
1944 
1945  KRATOS_CATCH("")
1946  }
1947 
1948  void ExecuteParticlesPritingTool( ModelPart& lagrangian_model_part, int input_filter_factor )
1949  {
1950  KRATOS_TRY
1951  //mfilter_factor; //we will only print one out of every "filter_factor" particles of the total particle list
1952 
1953  if(mparticle_printing_tool_initialized==false)
1954  {
1955  mfilter_factor=input_filter_factor;
1956 
1957  if(lagrangian_model_part.NodesBegin()-lagrangian_model_part.NodesEnd()>0)
1958  KRATOS_THROW_ERROR(std::logic_error, "AN EMPTY MODEL PART IS REQUIRED FOR THE PRINTING OF PARTICLES", "");
1959 
1960  lagrangian_model_part.AddNodalSolutionStepVariable(VELOCITY);
1961  lagrangian_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);
1962  lagrangian_model_part.AddNodalSolutionStepVariable(DISTANCE);
1963 
1964  for (unsigned int i=0; i!=((mmaximum_number_of_particles*mnelems)/mfilter_factor)+mfilter_factor; i++)
1965  {
1966  Node ::Pointer pnode = lagrangian_model_part.CreateNewNode( i+mlast_node_id+1 , 0.0, 0.0, 0.0); //recordar que es el nueevo model part!!
1967  //pnode->SetBufferSize(mr_model_part.NodesBegin()->GetBufferSize());
1968  pnode->SetBufferSize(1);
1969  }
1970  mparticle_printing_tool_initialized=true;
1971  }
1972 
1973  //resetting data of the unused particles
1974  const double inactive_particle_position= -10.0;
1975  array_1d<double,3>inactive_particle_position_vector;
1976  inactive_particle_position_vector(0)=inactive_particle_position;
1977  inactive_particle_position_vector(1)=inactive_particle_position;
1978  inactive_particle_position_vector(2)=inactive_particle_position;
1979  ModelPart::NodesContainerType::iterator inodebegin = lagrangian_model_part.NodesBegin();
1980  for(unsigned int ii=0; ii<lagrangian_model_part.Nodes().size(); ii++)
1981  {
1982  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
1983  inode->FastGetSolutionStepValue(DISTANCE) = 0.0;
1984  inode->FastGetSolutionStepValue(VELOCITY) = ZeroVector(3);
1985  inode->FastGetSolutionStepValue(DISPLACEMENT) = inactive_particle_position_vector;
1986  }
1987 
1988 
1989  int counter=0;
1990  //ModelPart::NodesContainerType::iterator it_begin = lagrangian_model_part.NodesBegin();
1991  for (int i=0; i!=mmaximum_number_of_particles*mnelems; i++)
1992  {
1993  PFEM_Particle_Fluid& pparticle =mparticles_vector[i];
1994  if(pparticle.GetEraseFlag()==false && i%mfilter_factor==0)
1995  {
1996  ModelPart::NodesContainerType::iterator inode = inodebegin+counter; //copying info from the particle to the (printing) node.
1997  inode->FastGetSolutionStepValue(DISTANCE) = pparticle.GetDistance();
1998  inode->FastGetSolutionStepValue(VELOCITY) = pparticle.GetVelocity();
1999  inode->FastGetSolutionStepValue(DISPLACEMENT) = pparticle.Coordinates();
2000  counter++;
2001  }
2002  }
2003 
2004  KRATOS_CATCH("")
2005 
2006  }
2007 
2008  void ExecuteParticlesPritingToolForDroppletsOnly( ModelPart& lagrangian_model_part, int input_filter_factor )
2009  {
2010  KRATOS_TRY
2011  //mfilter_factor; //we will only print one out of every "filter_factor" particles of the total particle list
2012  const int first_particle_id=1000000;
2013  if(mparticle_printing_tool_initialized==false)
2014  {
2015  mfilter_factor=input_filter_factor;
2016 
2017  if(lagrangian_model_part.NodesBegin()-lagrangian_model_part.NodesEnd()>0)
2018  KRATOS_THROW_ERROR(std::logic_error, "AN EMPTY MODEL PART IS REQUIRED FOR THE PRINTING OF PARTICLES", "");
2019 
2020  lagrangian_model_part.AddNodalSolutionStepVariable(VELOCITY);
2021  lagrangian_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);
2022  lagrangian_model_part.AddNodalSolutionStepVariable(DISTANCE);
2023 
2024  for (unsigned int i=0; i!=((mmaximum_number_of_particles*mnelems)/mfilter_factor)+mfilter_factor; i++)
2025  {
2026  Node ::Pointer pnode = lagrangian_model_part.CreateNewNode( i+first_particle_id+1 , 0.0, 0.0, 0.0); //recordar que es el nueevo model part!!
2027  //pnode->SetBufferSize(mr_model_part.NodesBegin()->GetBufferSize());
2028  pnode->SetBufferSize(1);
2029  }
2030  mparticle_printing_tool_initialized=true;
2031  }
2032 
2033  //resetting data of the unused particles
2034  const double inactive_particle_position= -10.0;
2035  array_1d<double,3>inactive_particle_position_vector;
2036  inactive_particle_position_vector(0)=inactive_particle_position;
2037  inactive_particle_position_vector(1)=inactive_particle_position;
2038  inactive_particle_position_vector(2)=inactive_particle_position;
2039  ModelPart::NodesContainerType::iterator inodebegin = lagrangian_model_part.NodesBegin();
2040  for(unsigned int ii=0; ii<lagrangian_model_part.Nodes().size(); ii++)
2041  {
2042  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
2043  inode->FastGetSolutionStepValue(DISTANCE) = 0.0;
2044  inode->FastGetSolutionStepValue(VELOCITY) = ZeroVector(3);
2045  inode->FastGetSolutionStepValue(DISPLACEMENT) = inactive_particle_position_vector;
2046  }
2047  const int max_number_of_printed_particles=lagrangian_model_part.Nodes().size();
2048 
2049 
2050  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
2051  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
2052  //(flag managed only by MoveParticles
2053  //KRATOS_WATCH(offset)
2054  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
2055 
2056  int counter=0;
2057  for(unsigned int ii=0; ii<mr_model_part.Elements().size(); ii++)
2058  {
2059  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
2060  Element::Pointer pelement(*ielem.base());
2061  Geometry<Node >& geom = ielem->GetGeometry();
2062  //double mean_elem_dist=0.0;
2063  bool pure_air_elem=true;
2064  for(unsigned int j=0; j<(TDim+1); j++)
2065  {
2066  if (geom[j].FastGetSolutionStepValue(DISTANCE)<0.0)
2067  pure_air_elem=false;
2068  //mean_elem_dist += geom[j].FastGetSolutionStepValue(DISTANCE);
2069  }
2070  //if (mean_elem_dist>0.0) //only air elements
2071  if (pure_air_elem==true)
2072  {
2073  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
2074  int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
2075  //std::cout << "elem " << ii << " with " << (unsigned int)number_of_particles_in_elem << " particles" << std::endl;
2076 
2077  for (int iii=0; iii<number_of_particles_in_elem ; iii++ )
2078  {
2079  //KRATOS_WATCH(iii)
2080  if (iii>mmaximum_number_of_particles) //it means we are out of our portion of the array, abort loop!
2081  break;
2082 
2083  PFEM_Particle_Fluid & pparticle = element_particle_pointers[offset+iii];
2084 
2085 
2086  bool erase_flag= pparticle.GetEraseFlag();
2087  if (erase_flag==false && pparticle.GetDistance()<0.0)
2088  {
2089  ModelPart::NodesContainerType::iterator inode = inodebegin+counter; //copying info from the particle to the (printing) node.
2090  inode->FastGetSolutionStepValue(DISTANCE) = pparticle.GetDistance();
2091  inode->FastGetSolutionStepValue(VELOCITY) = pparticle.GetVelocity();
2092  inode->FastGetSolutionStepValue(DISPLACEMENT) = pparticle.Coordinates();
2093  counter++;
2094  }
2095  }
2096  }
2097 
2098  if (counter>(max_number_of_printed_particles-30)) //we are approaching the end of the model part. so we stop before it's too late
2099  break;
2100  }
2101 
2102  KRATOS_CATCH("")
2103 
2104  }
2105 
2106  void AssignNodalVelocityUsingInletConditions(const double inlet_vel)
2107  {
2108  KRATOS_TRY
2109 
2110  //first we are going to delete all the velocities!
2111  ModelPart::ConditionsContainerType::iterator iconditionbegin = mr_model_part.ConditionsBegin();
2112  vector<unsigned int> condition_partition;
2113  #ifdef _OPENMP
2114  int number_of_threads = omp_get_max_threads();
2115  #else
2116  int number_of_threads = 1;
2117  #endif
2118 
2119  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Conditions().size(), condition_partition);
2120 
2121  #pragma omp parallel for
2122  for(int kkk=0; kkk<number_of_threads; kkk++)
2123  {
2124  for(unsigned int ii=condition_partition[kkk]; ii<condition_partition[kkk+1]; ii++)
2125  {
2126  ModelPart::ConditionsContainerType::iterator icondition = iconditionbegin+ii;
2127  if ( icondition->GetValue(IS_INLET) > 0.5 )
2128  {
2129  Geometry<Node >& geom = icondition->GetGeometry();
2130  array_1d<double,3> normal = ZeroVector(3);
2131  this->CalculateNormal(geom,normal);
2132  const double normal_lenght = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
2133  const array_1d<double,3> velocity = - inlet_vel/normal_lenght * normal;
2134  for (unsigned int l = 0; l < (TDim); l++)
2135  {
2136  geom[l].SetLock();
2137  geom[l].FastGetSolutionStepValue(VELOCITY) = velocity;
2138  geom[l].UnSetLock();
2139  }
2140  }
2141  }
2142  }
2143 
2144 
2145  KRATOS_CATCH("")
2146  }
2147 
2148  void RotateParticlesAndDomainVelocities(array_1d<double, 3 > rotations)
2149  {
2150  KRATOS_TRY
2151 
2152 
2153 
2154  if(fabs(rotations[0])>0.000000001 || fabs(rotations[1])>0.000000001)
2155  KRATOS_THROW_ERROR(std::invalid_argument,"ROTATIONS ONLY IMPLEMENTED AROUND Z AXIS! (xy plane) ","");
2156 
2157  const double cosinus_theta = cos(rotations[2]);
2158  const double sinus_theta = sin(rotations[2]);
2159 
2160  //std::cout << "updating particles" << std::endl;
2161  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
2162 
2163  const int offset = CurrentProcessInfo[WATER_PARTICLE_POINTERS_OFFSET]; //the array of pointers for each element has twice the required size so that we use a part in odd timesteps and the other in even ones.
2164  //(flag managed only by MoveParticles
2165  //KRATOS_WATCH(offset)
2166  ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.ElementsBegin();
2167 
2168 
2169  vector<unsigned int> element_partition;
2170  #ifdef _OPENMP
2171  int number_of_threads = omp_get_max_threads();
2172  #else
2173  int number_of_threads = 1;
2174  #endif
2175  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Elements().size(), element_partition);
2176 
2177  #pragma omp parallel for
2178  for(int kkk=0; kkk<number_of_threads; kkk++)
2179  {
2180  for(unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
2181  {
2182  //const int & elem_id = ielem->Id();
2183  ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
2184  Element::Pointer pelement(*ielem.base());
2185 
2186  ParticlePointerVector& element_particle_pointers = (ielem->GetValue(FLUID_PARTICLE_POINTERS));
2187  int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_FLUID_PARTICLES);
2188  //std::cout << "elem " << ii << " with " << (unsigned int)number_of_particles_in_elem << " particles" << std::endl;
2189 
2190  for (int iii=0; iii<number_of_particles_in_elem ; iii++ )
2191  {
2192  //KRATOS_WATCH(iii)
2193  if (iii>mmaximum_number_of_particles) //it means we are out of our portion of the array, abort loop!
2194  break;
2195 
2196  PFEM_Particle_Fluid & pparticle = element_particle_pointers[offset+iii];
2197 
2198 
2199  bool erase_flag= pparticle.GetEraseFlag();
2200  if (erase_flag==false)
2201  {
2202  array_1d<float, 3 > & vel = pparticle.GetVelocity();
2203  const float vel_x = vel[0];
2204  const float vel_y = vel[1];
2205  vel[0] = cosinus_theta*vel_x + sinus_theta*vel_y;
2206  vel[1] = cosinus_theta*vel_y - sinus_theta*vel_x;
2207  }
2208 
2209 
2210  }
2211  }
2212  }
2213 
2214 
2215  ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.NodesBegin();
2216  vector<unsigned int> node_partition;
2217  OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.Nodes().size(), node_partition);
2218 
2219  #pragma omp parallel for
2220  for(int kkk=0; kkk<number_of_threads; kkk++)
2221  {
2222  for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
2223  {
2224  ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
2225  if (inode->IsFixed(VELOCITY_X)==false)
2226  {
2227  array_1d<double, 3 > & vel = inode->FastGetSolutionStepValue(VELOCITY);
2228  const double vel_x = vel[0];
2229  const double vel_y = vel[1];
2230  vel[0] = cosinus_theta*vel_x + sinus_theta*vel_y;
2231  vel[1] = cosinus_theta*vel_y - sinus_theta*vel_x;
2232  }
2233  }
2234  }
2235  KRATOS_CATCH("")
2236  }
2237 
2238 
2239  protected:
2240 
2241 
2242  private:
2243 
2244 
2245  void Check()
2246  {
2247  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(DISTANCE) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing DISTANCE variable on solution step data","");
2248  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(VELOCITY) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing VELOCITY variable on solution step data","");
2249  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(PRESSURE) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing PRESSURE variable on solution step data","");
2250  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(PROJECTED_VELOCITY) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing PROJECTED_VELOCITY variable on solution step data","");
2251  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(DELTA_VELOCITY) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing DELTA_VELOCITY variable on solution step data","");
2252  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(MESH_VELOCITY) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing MESH_VELOCITY variable on solution step data","");
2253  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(YP) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing YP variable on solution step data","");
2254  if(mr_model_part.NodesBegin()->SolutionStepsDataHas(NORMAL) == false) KRATOS_THROW_ERROR(std::invalid_argument,"missing NORMAL variable on solution step data","");
2255  }
2256 
2257 
2261  void MoveParticle( PFEM_Particle_Fluid & pparticle,
2262  Element::Pointer & pelement,
2263  GlobalPointersVector< Element >& elements_in_trajectory,
2264  unsigned int & number_of_elements_in_trajectory,
2265  ResultIteratorType result_begin,
2266  const unsigned int MaxNumberOfResults,
2267  const array_1d<double,3> mesh_displacement,
2268  const bool discriminate_streamlines,
2269  const bool use_mesh_velocity_to_convect)
2270  {
2271 
2272  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
2273  double delta_t = CurrentProcessInfo[DELTA_TIME];
2274  array_1d<double,3> & gravity= CurrentProcessInfo[GRAVITY];
2275  unsigned int nsubsteps;
2276  double substep_dt;
2277 
2278 
2279  bool KEEP_INTEGRATING=false;
2280  bool is_found;
2281  //bool have_air_node;
2282  //bool have_water_node;
2283 
2284  array_1d<double,3> vel;
2285  array_1d<double,3> vel_without_other_phase_nodes=ZeroVector(3);
2286  array_1d<double,3> position;
2287  array_1d<double,3> mid_position;
2288  array_1d<double,TDim+1> N;
2289 
2290  //we start with the first position, then it will enter the loop.
2291  position = pparticle.Coordinates(); //initial coordinates
2292 
2293  const float particle_distance = pparticle.GetDistance();
2294  array_1d<float,3> particle_velocity = pparticle.GetVelocity();
2295  //double distance=0.0;
2296  array_1d<double,3> last_useful_vel;
2297  double sum_Ns_without_other_phase_nodes;
2298  //double pressure=0.0;
2300  //bool flying_water_particle=true; //if a water particle does not find a water element in its whole path, then we add the gravity*dt
2301  double only_integral = 0.0 ;
2302 
2303  is_found = FindNodeOnMesh(position, N ,pelement,result_begin,MaxNumberOfResults); //good, now we know where this point is:
2304  if(is_found == true)
2305  {
2306  KEEP_INTEGRATING=true;
2307  Geometry< Node >& geom = pelement->GetGeometry();//the element we're in
2308  vel=ZeroVector(3);
2309  vel_without_other_phase_nodes = ZeroVector(3);
2310  sum_Ns_without_other_phase_nodes=0.0;
2311  //distance=0.0;
2312 
2313  if (particle_distance<0.0 && discriminate_streamlines==true)
2314  {
2315  for(unsigned int j=0; j<(TDim+1); j++)
2316  {
2317  if ((geom[j].FastGetSolutionStepValue(DISTANCE))<0.0) //ok. useful info!
2318  {
2319  sum_Ns_without_other_phase_nodes += N[j];
2320  noalias(vel_without_other_phase_nodes) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2321  if (use_mesh_velocity_to_convect)
2322  noalias(vel_without_other_phase_nodes) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2323  }
2324 
2325  noalias(vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2326  if (use_mesh_velocity_to_convect)
2327  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2328  }
2329 
2330  if (sum_Ns_without_other_phase_nodes>0.01)
2331  {
2332  vel = vel_without_other_phase_nodes / sum_Ns_without_other_phase_nodes;
2333  //flying_water_particle=false;
2334  }
2335  else
2336  {
2337  vel = particle_velocity;
2338  if (use_mesh_velocity_to_convect)
2339  {
2340  for(unsigned int j=0; j<(TDim+1); j++)
2341  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2342  }
2343  }
2344  }
2345  else // air particle or we are not following streamlines
2346  {
2347  for(unsigned int j=0; j<(TDim+1); j++)
2348  {
2349  noalias(vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2350  if (use_mesh_velocity_to_convect)
2351  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2352  }
2353  //flying_water_particle=false;
2354  }
2355 
2356 
2357 
2358  //calculating substep to get +- courant(substep) = 0.1
2359  nsubsteps = 10.0 * (delta_t * pelement->GetValue(VELOCITY_OVER_ELEM_SIZE));
2360  if (nsubsteps<1)
2361  nsubsteps=1;
2362  substep_dt = delta_t / double(nsubsteps);
2363 
2364  only_integral = 1.0;// weight;//*double(nsubsteps);
2365 
2366 
2367  position += vel*substep_dt;//weight;
2368 
2370  last_useful_vel=vel;
2372 
2373  //DONE THE FIRST LOCATION OF THE PARTICLE, NOW WE PROCEED TO STREAMLINE INTEGRATION USING THE MESH VELOCITY
2375  unsigned int check_from_element_number=0;
2376 
2377 
2378  for(unsigned int i=0; i<(nsubsteps-1); i++)// this is for the substeps n+1. in the first one we already knew the position of the particle.
2379  {
2380  if (KEEP_INTEGRATING==true)
2381  {
2382  is_found = FindNodeOnMesh(position, N ,pelement,elements_in_trajectory,number_of_elements_in_trajectory,check_from_element_number,result_begin,MaxNumberOfResults); //good, now we know where this point is:
2383  if(is_found == true)
2384  {
2385  Geometry< Node >& geom = pelement->GetGeometry();//the element we're in
2386  sum_Ns_without_other_phase_nodes=0.0;
2387 
2388  if (particle_distance<0.0 && discriminate_streamlines==true)
2389  {
2390  vel_without_other_phase_nodes = ZeroVector(3);
2391 
2392  for(unsigned int j=0; j<TDim+1; j++)
2393  {
2394  if ((geom[j].FastGetSolutionStepValue(DISTANCE))<0.0) //ok. useful info!
2395  {
2396  sum_Ns_without_other_phase_nodes += N[j];
2397  noalias(vel_without_other_phase_nodes) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2398  if (use_mesh_velocity_to_convect)
2399  noalias(vel_without_other_phase_nodes) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2400  }
2401  noalias(vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2402  if (use_mesh_velocity_to_convect)
2403  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2404  }
2405 
2406  //if (have_water_node)
2407  //if (distance<0.0)
2408  if (sum_Ns_without_other_phase_nodes>0.01)
2409  {
2410  vel = vel_without_other_phase_nodes / sum_Ns_without_other_phase_nodes;
2411  //flying_water_particle=false;
2412  }
2413  else
2414  {
2415  particle_velocity += substep_dt * gravity;
2416  vel = particle_velocity;
2417  if (use_mesh_velocity_to_convect)
2418  {
2419  for(unsigned int j=0; j<(TDim+1); j++)
2420  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2421  }
2422  }
2423  }
2424  else //air particle or we are not discriminating streamlines
2425  {
2426  vel_without_other_phase_nodes = ZeroVector(3);
2427  vel = ZeroVector(3);
2428  for(unsigned int j=0; j<(TDim+1); j++)
2429  {
2430  noalias(vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2431  if (use_mesh_velocity_to_convect)
2432  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2433  }
2434 
2435  //flying_water_particle=false;
2436  }
2437 
2438  only_integral += 1.0; //values saved for the current time step
2439 
2440  position+=vel*substep_dt;//weight;
2441 
2442  }
2443  else
2444  {
2445  KEEP_INTEGRATING=false;
2446  break;
2447  }
2448  }
2449  else
2450  break;
2451 
2452 
2453  }
2454 
2455  }
2456 
2457  //if there's a mesh velocity, we add it at the end in a single step:
2458  position-=mesh_displacement;
2459 
2460 
2461  if (KEEP_INTEGRATING==false) (pparticle.GetEraseFlag()=true);
2462  else is_found = FindNodeOnMesh(position, N ,pelement,result_begin,MaxNumberOfResults); //we must save the pointer of the last element that we're in (inside the pointervector pelement)
2463 
2464  if (is_found==false) ( pparticle.GetEraseFlag()=true);
2465 
2466  pparticle.Coordinates() = position;
2467  }
2468 
2469 
2470  void AccelerateParticleUsingDeltaVelocity(
2471  PFEM_Particle_Fluid & pparticle,
2472  Element::Pointer & pelement,
2473  Geometry< Node >& geom)
2474  {
2475  array_1d<double,TDim+1> N;
2476 
2477  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
2478  const double delta_t = CurrentProcessInfo[DELTA_TIME];
2479  array_1d<double,3> gravity = CurrentProcessInfo[GRAVITY];
2480 
2481 
2482  //we start with the first position, then it will enter the loop.
2483  array_1d<double,3> coords = pparticle.Coordinates();
2484  float & particle_distance = pparticle.GetDistance();
2485  //double distance=0.0;
2486  array_1d<double,3> delta_velocity = ZeroVector(3);
2487 
2488  array_1d<double,3> delta_velocity_without_air = ZeroVector(3);
2489  array_1d<double,3> delta_velocity_without_water = ZeroVector(3);
2490  double sum_Ns_without_water_nodes = 0.0;
2491  double sum_Ns_without_air_nodes = 0.0;
2492 
2493  bool is_found = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2494  if(is_found == false)
2495  {
2496  KRATOS_WATCH(N)
2497  for (int j=0 ; j!=(TDim+1); j++)
2498  if (N[j]<0.0 )
2499  N[j]=1e-10;
2500  }
2501 
2502 
2503  if (particle_distance>0.0) //no problem. air
2504  {
2505  for(unsigned int j=0; j<(TDim+1); j++)
2506  {
2507  //just for air
2508  if ((geom[j].FastGetSolutionStepValue(DISTANCE))>0.0)
2509  {
2510  noalias(delta_velocity_without_water) += geom[j].FastGetSolutionStepValue(DELTA_VELOCITY)*N[j];
2511  sum_Ns_without_air_nodes += N[j];
2512 
2513  }
2514  //both air and water
2515  noalias(delta_velocity) += geom[j].FastGetSolutionStepValue(DELTA_VELOCITY)*N[j];
2516  }
2517 
2518  if (sum_Ns_without_water_nodes>0.01)
2519  {
2520  //delta_velocity = delta_velocity_without_water/sum_Ns_without_water_nodes ; //commented = using all the velocities always!
2521  }
2522  //else we use the complete field
2523 
2524  }
2525  else //water particle
2526  {
2527 
2528  for(unsigned int j=0; j<(TDim+1); j++)
2529  {
2530  if ((geom[j].FastGetSolutionStepValue(DISTANCE))<0.0)
2531  {
2532  noalias(delta_velocity_without_air) += geom[j].FastGetSolutionStepValue(DELTA_VELOCITY)*N[j];
2533  sum_Ns_without_air_nodes += N[j];
2534 
2535  }
2536  noalias(delta_velocity) += geom[j].FastGetSolutionStepValue(DELTA_VELOCITY)*N[j];
2537  }
2538 
2539  if (sum_Ns_without_air_nodes>0.01)
2540  {
2541  delta_velocity = delta_velocity_without_air/sum_Ns_without_air_nodes ;
2542  }
2543  else
2544  {
2545  if (mDENSITY_WATER>(10.0*mDENSITY_AIR))
2546  {
2547  delta_velocity=gravity*(1.0-mDENSITY_AIR/mDENSITY_WATER)*delta_t;
2548  }
2549  }
2550 
2551  }
2552  pparticle.GetVelocity() = pparticle.GetVelocity() + delta_velocity;
2553  }
2554 
2555  void MoveParticle_inverse_way(
2556  PFEM_Particle_Fluid & pparticle,
2557  Element::Pointer & pelement, //NOT A REFERENCE!! WE SHALL NOT OVERWRITE THE ELEMENT IT BELONGS TO!
2558  ResultIteratorType result_begin,
2559  const unsigned int MaxNumberOfResults,
2560  const bool use_mesh_velocity_to_convect)
2561  {
2562 
2563  ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo();
2564  double delta_t = CurrentProcessInfo[DELTA_TIME];
2565  unsigned int nsubsteps;
2566  double substep_dt;
2567 
2568 
2569  bool KEEP_INTEGRATING=false;
2570  bool is_found;
2571 
2572  array_1d<double,3> vel;
2573  array_1d<double,3> particle_vel;
2574  array_1d<double,3> position;
2575  array_1d<double,3> mid_position;
2576  array_1d<double,TDim+1> N;
2577 
2578 
2579  //we start with the first position, then it will enter the loop.
2580  position = pparticle.Coordinates(); // + (pparticle)->FastGetSolutionStepValue(DISPLACEMENT); //initial coordinates
2581 
2582 
2583  float & distance = pparticle.GetDistance();
2584  double only_integral = 0.0 ;
2585 
2586  is_found = FindNodeOnMesh(position, N ,pelement,result_begin,MaxNumberOfResults); //good, now we know where this point is:
2587  if(is_found == true)
2588  {
2589  KEEP_INTEGRATING=true;
2590  Geometry< Node >& geom = pelement->GetGeometry();//the element we're in
2591  vel=ZeroVector(3);
2592  particle_vel=ZeroVector(3);
2593  distance=0.0;
2594 
2595  for(unsigned int j=0; j<(TDim+1); j++)
2596  {
2597  distance += geom[j].FastGetSolutionStepValue(DISTANCE)*N(j);
2598  noalias(particle_vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2599  noalias(vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j];
2600  if (use_mesh_velocity_to_convect)
2601  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2602  }
2603  //calculating substep to get +- courant(substep) = 1/4
2604  nsubsteps = 10.0 * (delta_t * pelement->GetValue(VELOCITY_OVER_ELEM_SIZE));
2605  if (nsubsteps<1)
2606  nsubsteps=1;
2607  substep_dt = delta_t / double(nsubsteps);
2608 
2609  only_integral = 1.0;// weight;//*double(nsubsteps);
2610  position -= vel*substep_dt;//weight;
2611 
2612  for(unsigned int i=0; i<(nsubsteps-1); i++)// this is for the substeps n+1. in the first one we already knew the position of the particle.
2613  { if (KEEP_INTEGRATING==true) {
2614  is_found = FindNodeOnMesh(position, N ,pelement,result_begin,MaxNumberOfResults); //good, now we know where this point is:
2615  if(is_found == true)
2616  {
2617  Geometry< Node >& geom = pelement->GetGeometry();//the element we're in
2618 
2619  vel=ZeroVector(3);
2620  particle_vel=ZeroVector(3);
2621  distance=0.0;
2622 
2623 
2624  for(unsigned int j=0; j<(TDim+1); j++)
2625  {
2626  noalias(particle_vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j] ;
2627  noalias(vel) += geom[j].FastGetSolutionStepValue(VELOCITY)*N[j] ;
2628  distance += geom[j].FastGetSolutionStepValue(DISTANCE)*N(j);
2629  if (use_mesh_velocity_to_convect)
2630  noalias(vel) -= geom[j].FastGetSolutionStepValue(MESH_VELOCITY)*N[j];
2631  }
2632 
2633 
2634  only_integral += 1.0;//weight ; //values saved for the current time step
2635  position-=vel*substep_dt;//weight;
2636  }
2637  else KEEP_INTEGRATING=false;
2638  }
2639 
2640 
2641  }
2642 
2644  if(distance>0.0)
2645  {
2646  //if(distance<2.0)
2647  distance=1.0;
2648  //else
2649  // distance=3.0;
2650  }
2651  else
2652  distance=-1.0;
2653 
2654 
2655 
2656  pparticle.GetVelocity()=particle_vel;
2657  }
2658  //else {KRATOS_WATCH(position); }
2659 
2660  }
2661 
2662  void OverwriteParticleDataUsingTopographicDomain(
2663  PFEM_Particle_Fluid & pparticle,
2664  Element::Pointer & pelement,
2665  array_1d<double,3> domains_offset,
2666  ResultIteratorType result_begin,
2667  const unsigned int MaxNumberOfResults)
2668  {
2669  array_1d<double,TDim+1> N;
2670 
2671  //we start with the first position, then it will enter the loop.
2672  array_1d<double,3> coords = pparticle.Coordinates()+domains_offset;
2673  float & particle_distance = pparticle.GetDistance();
2674  bool is_found = FindNodeOnTopographicMesh(coords, N ,pelement,result_begin,MaxNumberOfResults); //good, now we know where this point is:
2675 
2676  if (is_found) //it is part of the solid topographic domain
2677  {
2678  particle_distance= -1.0;
2679  }
2680  else //it is outside the topographic domain, therefore it is air or whatever it means
2681  {
2682  particle_distance= 1.0;
2683  }
2684 
2685  pparticle.GetVelocity() = ZeroVector(3);
2686  }
2687 
2688 
2689 
2690 
2695  bool FindNodeOnMesh( array_1d<double,3>& position,
2696  array_1d<double,TDim+1>& N,
2697  Element::Pointer & pelement,
2698  ResultIteratorType result_begin,
2699  const unsigned int MaxNumberOfResults)
2700  {
2701  typedef std::size_t SizeType;
2702 
2703  const array_1d<double,3>& coords = position;
2704  array_1d<double,TDim+1> aux_N;
2705  //before using the bin to search for possible elements we check first the last element in which the particle was.
2706  Geometry<Node >& geom_default = pelement->GetGeometry(); //(*(i))->GetGeometry();
2707  bool is_found_1 = CalculatePosition(geom_default,coords[0],coords[1],coords[2],N);
2708  if(is_found_1 == true) //that was easy!
2709  {
2710  return true;
2711  }
2712 
2713  //to begin with we check the neighbour elements; it is a bit more expensive
2714  GlobalPointersVector< Element >& neighb_elems = pelement->GetValue(NEIGHBOUR_ELEMENTS);
2715  //the first we check is the one that has negative shape function, because it means it went outside in this direction:
2716  //commented, it is not faster than simply checking all the neighbours (branching)
2717  /*
2718  unsigned int checked_element=0;
2719  for (unsigned int i=0;i!=(TDim+1);i++)
2720  {
2721  if (N[i]<0.0)
2722  {
2723  checked_element=i;
2724  Geometry<Node >& geom = neighb_elems[i].GetGeometry();
2725  bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],aux_N);
2726  if (is_found_2)
2727  {
2728  pelement=Element::Pointer(((neighb_elems(i))));
2729  N=aux_N;
2730  return true;
2731  }
2732  break;
2733  }
2734  }
2735  */
2736  //we check all the neighbour elements
2737  for (unsigned int i=0;i!=(neighb_elems.size());i++)
2738  {
2739  if(neighb_elems(i).get()!=nullptr)
2740  {
2741  Geometry<Node >& geom = neighb_elems[i].GetGeometry();
2742  bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2743  if (is_found_2)
2744  {
2745  pelement = neighb_elems[i].shared_from_this();
2746  return true;
2747  }
2748  }
2749  }
2750 
2751  //if checking all the neighbour elements did not work, we have to use the bins
2752  //ask to the container for the list of candidate elements
2753  SizeType results_found = mpBinsObjectDynamic->SearchObjectsInCell(Point{coords}, result_begin, MaxNumberOfResults );
2754 
2755  if(results_found>0){
2756  //loop over the candidate elements and check if the particle falls within
2757  for(SizeType i = 0; i< results_found; i++)
2758  {
2759  Geometry<Node >& geom = (*(result_begin+i))->GetGeometry();
2760 
2761  //find local position
2762  bool is_found = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2763 
2764  if(is_found == true)
2765  {
2766  pelement=Element::Pointer((*(result_begin+i)));
2767  return true;
2768  }
2769  }
2770  }
2771  //if nothing worked, then:
2772  //not found case
2773  return false;
2774  }
2775 
2776 
2777  // VERSION INCLUDING PREDEFINED ELEMENTS FOLLOWING A TRAJECTORY
2778  bool FindNodeOnMesh( array_1d<double,3>& position,
2779  array_1d<double,TDim+1>& N,
2780  Element::Pointer & pelement,
2781  GlobalPointersVector< Element >& elements_in_trajectory,
2782  unsigned int & number_of_elements_in_trajectory,
2783  unsigned int & check_from_element_number,
2784  ResultIteratorType result_begin,
2785  const unsigned int MaxNumberOfResults)
2786  {
2787  typedef std::size_t SizeType;
2788 
2789  const array_1d<double,3>& coords = position;
2790  array_1d<double,TDim+1> aux_N;
2791  //before using the bin to search for possible elements we check first the last element in which the particle was.
2792  Geometry<Node >& geom_default = pelement->GetGeometry(); //(*(i))->GetGeometry();
2793  bool is_found_1 = CalculatePosition(geom_default,coords[0],coords[1],coords[2],N);
2794  if(is_found_1 == true)
2795  {
2796  return true; //that was easy!
2797  }
2798 
2799  //if it was not found in the first element, we can proceed to check in the following elements (in the trajectory defined by previous particles that started from the same element.
2800  for (unsigned int i=(check_from_element_number);i!=number_of_elements_in_trajectory;i++)
2801  {
2802  Geometry<Node >& geom = elements_in_trajectory[i].GetGeometry();
2803  bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],aux_N);
2804  if (is_found_2)
2805  {
2806  pelement = elements_in_trajectory[i].shared_from_this();
2807  N=aux_N;
2808  check_from_element_number = i+1 ; //now i element matches pelement, so to avoid cheching twice the same element we send the counter to the following element.
2809  return true;
2810  }
2811 
2812  }
2813 
2814  //now we check the neighbour elements:
2815  GlobalPointersVector< Element >& neighb_elems = pelement->GetValue(NEIGHBOUR_ELEMENTS);
2816  //the first we check is the one that has negative shape function, because it means it went outside in this direction:
2817  //commented, it is not faster than simply checking all the neighbours (branching)
2818  /*
2819  unsigned int checked_element=0;
2820  for (unsigned int i=0;i!=(TDim+1);i++)
2821  {
2822  if (N[i]<0.0)
2823  {
2824  checked_element=i;
2825  Geometry<Node >& geom = neighb_elems[i].GetGeometry();
2826  bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],aux_N);
2827  if (is_found_2)
2828  {
2829  pelement=Element::Pointer(((neighb_elems(i))));
2830  N=aux_N;
2831  return true;
2832  }
2833  break;
2834  }
2835  }
2836  */
2837  //we check all the neighbour elements
2838  for (unsigned int i=0;i!=(neighb_elems.size());i++)
2839  {
2840  if(neighb_elems(i).get()!=nullptr)
2841  {
2842  Geometry<Node >& geom = neighb_elems[i].GetGeometry();
2843  bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2844  if (is_found_2)
2845  {
2846  pelement = neighb_elems[i].shared_from_this();
2847  if (number_of_elements_in_trajectory<20)
2848  {
2849  elements_in_trajectory(number_of_elements_in_trajectory)=pelement;
2850  number_of_elements_in_trajectory++;
2851  check_from_element_number = number_of_elements_in_trajectory; //we do it after doing the ++ to the counter, so we woudlnt enter the loop that searches in the elements_in_trajectory list. we are the particle that is adding elements to the list
2852  }
2853  return true;
2854  }
2855  }
2856  }
2857 
2858 
2859  //if checking all the neighbour elements did not work, we have to use the bins
2860  //ask to the container for the list of candidate elements
2861  SizeType results_found = mpBinsObjectDynamic->SearchObjectsInCell(Point{coords}, result_begin, MaxNumberOfResults );
2862 
2863  if(results_found>0)
2864  {
2865  //loop over the candidate elements and check if the particle falls within
2866  for(SizeType i = 0; i< results_found; i++)
2867  {
2868  Geometry<Node >& geom = (*(result_begin+i))->GetGeometry();
2869 
2870  //find local position
2871  bool is_found = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2872 
2873  if(is_found == true)
2874  {
2875  pelement=Element::Pointer((*(result_begin+i)));
2876  if (number_of_elements_in_trajectory<20)
2877  {
2878  elements_in_trajectory(number_of_elements_in_trajectory)=pelement;
2879  number_of_elements_in_trajectory++;
2880  check_from_element_number = number_of_elements_in_trajectory; //we do it after doing the ++ to the counter, so we woudlnt enter the loop that searches in the elements_in_trajectory list. we are the particle that is adding elements to the list
2881  }
2882  return true;
2883  }
2884  }
2885  }
2886 
2887  //not found case
2888  return false;
2889  }
2890 
2891 
2896  bool FindNodeOnTopographicMesh( array_1d<double,3>& position,
2897  array_1d<double,TDim+1>& N,
2898  Element::Pointer & pelement,
2899  ResultIteratorType result_begin,
2900  const unsigned int MaxNumberOfResults)
2901  {
2902  typedef std::size_t SizeType;
2903 
2904  const array_1d<double,3>& coords = position;
2905  array_1d<double,TDim+1> aux_N;
2906  //before using the bin to search for possible elements we check first the last element in which the particle was.
2907 
2908  //ModelPart::ElementsContainerType::iterator i = mr_model_part.ElementsBegin()+last_element;
2909  Geometry<Node >& geom_default = pelement->GetGeometry(); //(*(i))->GetGeometry();
2910  bool is_found_1 = CalculatePosition(geom_default,coords[0],coords[1],coords[2],N);
2911  if(is_found_1 == true)
2912  {
2913  //pelement = (*(i));
2914  return true;
2915  }
2916 
2917  //to begin with we check the neighbour elements:
2918  GlobalPointersVector< Element >& neighb_elems = pelement->GetValue(NEIGHBOUR_ELEMENTS);
2919  for (unsigned int i=0;i!=(neighb_elems.size());i++)
2920  {
2921  if(neighb_elems(i).get()!=nullptr)
2922  {
2923  Geometry<Node >& geom = neighb_elems[i].GetGeometry();
2924  bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2925  if (is_found_2)
2926  {
2927  pelement = neighb_elems[i].shared_from_this();
2928  return true;
2929  }
2930  }
2931  }
2932 
2933 
2934  //ask to the container for the list of candidate elements
2935  SizeType results_found = mpTopographicBinsObjectDynamic->SearchObjectsInCell(Point{coords}, result_begin, MaxNumberOfResults );
2936  //KRATOS_WATCH(results_found)
2937 
2938  if(results_found>0){
2939  //loop over the candidate elements and check if the particle falls within
2940  for(SizeType i = 0; i< results_found; i++)
2941  {
2942  Geometry<Node >& geom = (*(result_begin+i))->GetGeometry();
2943 
2944  //find local position
2945  bool is_found = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
2946 
2947  if(is_found == true)
2948  {
2949  pelement=Element::Pointer((*(result_begin+i)));
2950  return true;
2951  }
2952  }
2953  }
2954 
2955  //not found case
2956  return false;
2957  }
2958 
2959  //***************************************
2960  //***************************************
2961 
2962  inline bool CalculatePosition(Geometry<Node >&geom,
2963  const double xc, const double yc, const double zc,
2964  array_1d<double, 3 > & N
2965  )
2966  {
2967  double x0 = geom[0].X();
2968  double y0 = geom[0].Y();
2969  double x1 = geom[1].X();
2970  double y1 = geom[1].Y();
2971  double x2 = geom[2].X();
2972  double y2 = geom[2].Y();
2973 
2974  double area = CalculateVol(x0, y0, x1, y1, x2, y2);
2975  double inv_area = 0.0;
2976  if (area == 0.0)
2977  {
2978  KRATOS_THROW_ERROR(std::logic_error, "element with zero area found", "");
2979  } else
2980  {
2981  inv_area = 1.0 / area;
2982  }
2983 
2984 
2985  N[0] = CalculateVol(x1, y1, x2, y2, xc, yc) * inv_area;
2986  N[1] = CalculateVol(x2, y2, x0, y0, xc, yc) * inv_area;
2987  N[2] = CalculateVol(x0, y0, x1, y1, xc, yc) * inv_area;
2988  //KRATOS_WATCH(N);
2989 
2990  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
2991  return true;
2992 
2993  return false;
2994  }
2996  //using the pre loaded nodal coordinates
2997  inline bool CalculatePosition(const array_1d<double,3*(TDim+1)>& nodes_positions,
2998  const double xc, const double yc, const double zc,
2999  array_1d<double, 3 > & N
3000  )
3001  {
3002  const double& x0 = nodes_positions[0];
3003  const double& y0 = nodes_positions[1];
3004  const double& x1 = nodes_positions[3];
3005  const double& y1 = nodes_positions[4];
3006  const double& x2 = nodes_positions[6];
3007  const double& y2 = nodes_positions[7];
3008 
3009  double area = CalculateVol(x0, y0, x1, y1, x2, y2);
3010  double inv_area = 0.0;
3011  if (area == 0.0)
3012  {
3013  KRATOS_THROW_ERROR(std::logic_error, "element with zero area found", "");
3014  } else
3015  {
3016  inv_area = 1.0 / area;
3017  }
3018 
3019 
3020  N[0] = CalculateVol(x1, y1, x2, y2, xc, yc) * inv_area;
3021  N[1] = CalculateVol(x2, y2, x0, y0, xc, yc) * inv_area;
3022  N[2] = CalculateVol(x0, y0, x1, y1, xc, yc) * inv_area;
3023  //KRATOS_WATCH(N);
3024 
3025  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
3026  return true;
3027 
3028  return false;
3029  }
3030 
3031 
3032  //***************************************
3033  //***************************************
3034 
3035  inline bool CalculatePosition(Geometry<Node >&geom,
3036  const double xc, const double yc, const double zc,
3037  array_1d<double, 4 > & N
3038  )
3039  {
3040 
3041  double x0 = geom[0].X();
3042  double y0 = geom[0].Y();
3043  double z0 = geom[0].Z();
3044  double x1 = geom[1].X();
3045  double y1 = geom[1].Y();
3046  double z1 = geom[1].Z();
3047  double x2 = geom[2].X();
3048  double y2 = geom[2].Y();
3049  double z2 = geom[2].Z();
3050  double x3 = geom[3].X();
3051  double y3 = geom[3].Y();
3052  double z3 = geom[3].Z();
3053 
3054  double vol = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3);
3055 
3056  double inv_vol = 0.0;
3057  if (vol < 0.000000000000000000000000000001)
3058  {
3059  KRATOS_THROW_ERROR(std::logic_error, "element with zero vol found", "");
3060  } else
3061  {
3062  inv_vol = 1.0 / vol;
3063  }
3064 
3065  N[0] = CalculateVol(x1, y1, z1, x3, y3, z3, x2, y2, z2, xc, yc, zc) * inv_vol;
3066  N[1] = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, xc, yc, zc) * inv_vol;
3067  N[2] = CalculateVol(x3, y3, z3, x1, y1, z1, x0, y0, z0, xc, yc, zc) * inv_vol;
3068  N[3] = CalculateVol(x3, y3, z3, x0, y0, z0, x2, y2, z2, xc, yc, zc) * inv_vol;
3069 
3070 
3071  if (N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[3] >= 0.0 &&
3072  N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0 && N[3] <= 1.0)
3073  //if the xc yc zc is inside the tetrahedron return true
3074  return true;
3075 
3076  return false;
3077  }
3079  //using the pre loaded nodal coordinates
3080  inline bool CalculatePosition(const array_1d<double,3*(TDim+1)>& nodes_positions,
3081  const double xc, const double yc, const double zc,
3082  array_1d<double, 4 > & N
3083  )
3084  {
3085 
3086  const double& x0 = nodes_positions[0];
3087  const double& y0 = nodes_positions[1];
3088  const double& z0 = nodes_positions[2];
3089  const double& x1 = nodes_positions[3];
3090  const double& y1 = nodes_positions[4];
3091  const double& z1 = nodes_positions[5];
3092  const double& x2 = nodes_positions[6];
3093  const double& y2 = nodes_positions[7];
3094  const double& z2 = nodes_positions[8];
3095  const double& x3 = nodes_positions[9];
3096  const double& y3 = nodes_positions[10];
3097  const double& z3 = nodes_positions[11];
3098 
3099  double vol = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3);
3100 
3101  double inv_vol = 0.0;
3102  if (vol < 0.000000000000000000000000000001)
3103  {
3104  KRATOS_THROW_ERROR(std::logic_error, "element with zero vol found", "");
3105  } else
3106  {
3107  inv_vol = 1.0 / vol;
3108  }
3109 
3110  N[0] = CalculateVol(x1, y1, z1, x3, y3, z3, x2, y2, z2, xc, yc, zc) * inv_vol;
3111  N[1] = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, xc, yc, zc) * inv_vol;
3112  N[2] = CalculateVol(x3, y3, z3, x1, y1, z1, x0, y0, z0, xc, yc, zc) * inv_vol;
3113  N[3] = CalculateVol(x3, y3, z3, x0, y0, z0, x2, y2, z2, xc, yc, zc) * inv_vol;
3114 
3115 
3116  if (N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[3] >= 0.0 &&
3117  N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0 && N[3] <= 1.0)
3118  //if the xc yc zc is inside the tetrahedron return true
3119  return true;
3120 
3121  return false;
3122  }
3123 
3124  inline double CalculateVol(const double x0, const double y0,
3125  const double x1, const double y1,
3126  const double x2, const double y2
3127  )
3128  {
3129  return 0.5 * ((x1 - x0)*(y2 - y0)- (y1 - y0)*(x2 - x0));
3130  }
3131  //***************************************
3132  //***************************************
3133 
3134  inline double CalculateVol(const double x0, const double y0, const double z0,
3135  const double x1, const double y1, const double z1,
3136  const double x2, const double y2, const double z2,
3137  const double x3, const double y3, const double z3
3138  )
3139  {
3140  double x10 = x1 - x0;
3141  double y10 = y1 - y0;
3142  double z10 = z1 - z0;
3143 
3144  double x20 = x2 - x0;
3145  double y20 = y2 - y0;
3146  double z20 = z2 - z0;
3147 
3148  double x30 = x3 - x0;
3149  double y30 = y3 - y0;
3150  double z30 = z3 - z0;
3151 
3152  double detJ = x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30;
3153  return detJ * 0.1666666666666666666667;
3154  }
3155 
3156 
3157 
3158  void ComputeGaussPointPositions_4(Geometry< Node >& geom, BoundedMatrix<double, 7, 3 > & pos,BoundedMatrix<double, 7, 3 > & N)
3159  {
3160  double one_third = 1.0 / 3.0;
3161  double one_sixt = 0.15; //1.0 / 6.0;
3162  double two_third = 0.7; //2.0 * one_third;
3163 
3164  N(0, 0) = one_sixt;
3165  N(0, 1) = one_sixt;
3166  N(0, 2) = two_third;
3167  N(1, 0) = two_third;
3168  N(1, 1) = one_sixt;
3169  N(1, 2) = one_sixt;
3170  N(2, 0) = one_sixt;
3171  N(2, 1) = two_third;
3172  N(2, 2) = one_sixt;
3173  N(3, 0) = one_third;
3174  N(3, 1) = one_third;
3175  N(3, 2) = one_third;
3176 
3177  //first
3178  pos(0, 0) = one_sixt * geom[0].X() + one_sixt * geom[1].X() + two_third * geom[2].X();
3179  pos(0, 1) = one_sixt * geom[0].Y() + one_sixt * geom[1].Y() + two_third * geom[2].Y();
3180  pos(0, 2) = one_sixt * geom[0].Z() + one_sixt * geom[1].Z() + two_third * geom[2].Z();
3181 
3182  //second
3183  pos(1, 0) = two_third * geom[0].X() + one_sixt * geom[1].X() + one_sixt * geom[2].X();
3184  pos(1, 1) = two_third * geom[0].Y() + one_sixt * geom[1].Y() + one_sixt * geom[2].Y();
3185  pos(1, 2) = two_third * geom[0].Z() + one_sixt * geom[1].Z() + one_sixt * geom[2].Z();
3186 
3187  //third
3188  pos(2, 0) = one_sixt * geom[0].X() + two_third * geom[1].X() + one_sixt * geom[2].X();
3189  pos(2, 1) = one_sixt * geom[0].Y() + two_third * geom[1].Y() + one_sixt * geom[2].Y();
3190  pos(2, 2) = one_sixt * geom[0].Z() + two_third * geom[1].Z() + one_sixt * geom[2].Z();
3191 
3192  //fourth
3193  pos(3, 0) = one_third * geom[0].X() + one_third * geom[1].X() + one_third * geom[2].X();
3194  pos(3, 1) = one_third * geom[0].Y() + one_third * geom[1].Y() + one_third * geom[2].Y();
3195  pos(3, 2) = one_third * geom[0].Z() + one_third * geom[1].Z() + one_third * geom[2].Z();
3196 
3197  }
3198 
3199 
3200  void ComputeGaussPointPositionsForPostReseed(Geometry< Node >& geom, BoundedMatrix<double, 7, 3 > & pos,BoundedMatrix<double, 7, 3 > & N) //2d
3201  {
3202  double one_third = 1.0 / 3.0;
3203  double one_eight = 0.12; //1.0 / 6.0;
3204  double three_quarters = 0.76; //2.0 * one_third;
3205 
3206  N(0, 0) = one_eight;
3207  N(0, 1) = one_eight;
3208  N(0, 2) = three_quarters;
3209 
3210  N(1, 0) = three_quarters;
3211  N(1, 1) = one_eight;
3212  N(1, 2) = one_eight;
3213 
3214  N(2, 0) = one_eight;
3215  N(2, 1) = three_quarters;
3216  N(2, 2) = one_eight;
3217 
3218  N(3, 0) = one_third;
3219  N(3, 1) = one_third;
3220  N(3, 2) = one_third;
3221 
3222  N(4, 0) = one_eight;
3223  N(4, 1) = 0.44;
3224  N(4, 2) = 0.44;
3225 
3226  N(5, 0) = 0.44;
3227  N(5, 1) = one_eight;
3228  N(5, 2) = 0.44;
3229 
3230  N(6, 0) = 0.44;
3231  N(6, 1) = 0.44;
3232  N(6, 2) = one_eight;
3233 
3234 
3235  //first
3236  pos(0, 0) = one_eight * geom[0].X() + one_eight * geom[1].X() + three_quarters * geom[2].X();
3237  pos(0, 1) = one_eight * geom[0].Y() + one_eight * geom[1].Y() + three_quarters * geom[2].Y();
3238  pos(0, 2) = one_eight * geom[0].Z() + one_eight * geom[1].Z() + three_quarters * geom[2].Z();
3239 
3240  //second
3241  pos(1, 0) = three_quarters * geom[0].X() + one_eight * geom[1].X() + one_eight * geom[2].X();
3242  pos(1, 1) = three_quarters * geom[0].Y() + one_eight * geom[1].Y() + one_eight * geom[2].Y();
3243  pos(1, 2) = three_quarters * geom[0].Z() + one_eight * geom[1].Z() + one_eight * geom[2].Z();
3244 
3245  //third
3246  pos(2, 0) = one_eight * geom[0].X() + three_quarters * geom[1].X() + one_eight * geom[2].X();
3247  pos(2, 1) = one_eight * geom[0].Y() + three_quarters * geom[1].Y() + one_eight * geom[2].Y();
3248  pos(2, 2) = one_eight * geom[0].Z() + three_quarters * geom[1].Z() + one_eight * geom[2].Z();
3249 
3250  //fourth
3251  pos(3, 0) = one_third * geom[0].X() + one_third * geom[1].X() + one_third * geom[2].X();
3252  pos(3, 1) = one_third * geom[0].Y() + one_third * geom[1].Y() + one_third * geom[2].Y();
3253  pos(3, 2) = one_third * geom[0].Z() + one_third * geom[1].Z() + one_third * geom[2].Z();
3254 
3255  //fifth
3256  pos(4, 0) = one_eight * geom[0].X() + 0.44 * geom[1].X() + 0.44 * geom[2].X();
3257  pos(4, 1) = one_eight * geom[0].Y() + 0.44 * geom[1].Y() + 0.44 * geom[2].Y();
3258  pos(4, 2) = one_eight * geom[0].Z() + 0.44 * geom[1].Z() + 0.44 * geom[2].Z();
3259 
3260  //sixth
3261  pos(5, 0) = 0.44 * geom[0].X() + one_eight * geom[1].X() + 0.44 * geom[2].X();
3262  pos(5, 1) = 0.44 * geom[0].Y() + one_eight * geom[1].Y() + 0.44 * geom[2].Y();
3263  pos(5, 2) = 0.44 * geom[0].Z() + one_eight * geom[1].Z() + 0.44 * geom[2].Z();
3264 
3265  //seventh
3266  pos(6, 0) = 0.44 * geom[0].X() + 0.44 * geom[1].X() + one_eight * geom[2].X();
3267  pos(6, 1) = 0.44 * geom[0].Y() + 0.44 * geom[1].Y() + one_eight * geom[2].Y();
3268  pos(6, 2) = 0.44 * geom[0].Z() + 0.44 * geom[1].Z() + one_eight * geom[2].Z();
3269 
3270 
3271 
3272 
3273  }
3274 
3275  void ComputeGaussPointPositionsForPostReseed(Geometry< Node >& geom, BoundedMatrix<double, 9, 3 > & pos,BoundedMatrix<double, 9, 4 > & N) //3D
3276  {
3277  double one_quarter = 0.25;
3278  double small_fraction = 0.1; //1.0 / 6.0;
3279  double big_fraction = 0.7; //2.0 * one_third;
3280  double mid_fraction = 0.3; //2.0 * one_third;
3281 
3282  N(0, 0) = big_fraction;
3283  N(0, 1) = small_fraction;
3284  N(0, 2) = small_fraction;
3285  N(0, 3) = small_fraction;
3286 
3287  N(1, 0) = small_fraction;
3288  N(1, 1) = big_fraction;
3289  N(1, 2) = small_fraction;
3290  N(1, 3) = small_fraction;
3291 
3292  N(2, 0) = small_fraction;
3293  N(2, 1) = small_fraction;
3294  N(2, 2) = big_fraction;
3295  N(2, 3) = small_fraction;
3296 
3297  N(3, 0) = small_fraction;
3298  N(3, 1) = small_fraction;
3299  N(3, 2) = small_fraction;
3300  N(3, 3) = big_fraction;
3301 
3302  N(4, 0) = one_quarter;
3303  N(4, 1) = one_quarter;
3304  N(4, 2) = one_quarter;
3305  N(4, 3) = one_quarter;
3306 
3307  N(5, 0) = small_fraction;
3308  N(5, 1) = mid_fraction;
3309  N(5, 2) = mid_fraction;
3310  N(5, 3) = mid_fraction;
3311 
3312  N(6, 0) = mid_fraction;
3313  N(6, 1) = small_fraction;
3314  N(6, 2) = mid_fraction;
3315  N(6, 3) = mid_fraction;
3316 
3317  N(7, 0) = mid_fraction;
3318  N(7, 1) = mid_fraction;
3319  N(7, 2) = small_fraction;
3320  N(7, 3) = mid_fraction;
3321 
3322  N(8, 0) = mid_fraction;
3323  N(8, 1) = mid_fraction;
3324  N(8, 2) = mid_fraction;
3325  N(8, 3) = small_fraction;
3326 
3327  pos=ZeroMatrix(9,3);
3328  for (unsigned int i=0; i!=4; i++) //going through the 4 nodes
3329  {
3330  array_1d<double, 3 > & coordinates = geom[i].Coordinates();
3331  for (unsigned int j=0; j!=9; j++) //going through the 9 particles
3332  {
3333  for (unsigned int k=0; k!=3; k++) //x,y,z
3334  pos(j,k) += N(j,i) * coordinates[k];
3335  }
3336  }
3337 
3338 
3339  }
3340 
3341 
3342 
3343  void ComputeGaussPointPositionsForPreReseed(Geometry< Node >& geom, BoundedMatrix<double, 3, 3 > & pos,BoundedMatrix<double, 3, 3 > & N) //2D
3344  {
3345 
3346  N(0, 0) = 0.5;
3347  N(0, 1) = 0.25;
3348  N(0, 2) = 0.25;
3349 
3350  N(1, 0) = 0.25;
3351  N(1, 1) = 0.5;
3352  N(1, 2) = 0.25;
3353 
3354  N(2, 0) = 0.25;
3355  N(2, 1) = 0.25;
3356  N(2, 2) = 0.5;
3357 
3358  //first
3359  pos(0, 0) = 0.5 * geom[0].X() + 0.25 * geom[1].X() + 0.25 * geom[2].X();
3360  pos(0, 1) = 0.5 * geom[0].Y() + 0.25 * geom[1].Y() + 0.25 * geom[2].Y();
3361  pos(0, 2) = 0.5 * geom[0].Z() + 0.25 * geom[1].Z() + 0.25 * geom[2].Z();
3362 
3363  //second
3364  pos(1, 0) = 0.25 * geom[0].X() + 0.5 * geom[1].X() + 0.25 * geom[2].X();
3365  pos(1, 1) = 0.25 * geom[0].Y() + 0.5 * geom[1].Y() + 0.25 * geom[2].Y();
3366  pos(1, 2) = 0.25 * geom[0].Z() + 0.5 * geom[1].Z() + 0.25 * geom[2].Z();
3367 
3368  //third
3369  pos(2, 0) = 0.25 * geom[0].X() + 0.25 * geom[1].X() + 0.5 * geom[2].X();
3370  pos(2, 1) = 0.25 * geom[0].Y() + 0.25 * geom[1].Y() + 0.5 * geom[2].Y();
3371  pos(2, 2) = 0.25 * geom[0].Z() + 0.25 * geom[1].Z() + 0.5 * geom[2].Z();
3372 
3373  }
3374 
3375  void ComputeGaussPointPositionsForPreReseed(Geometry< Node >& geom, BoundedMatrix<double, 4, 3 > & pos,BoundedMatrix<double, 4, 4 > & N) //3D
3376  {
3377  //creating 4 particles, each will be closer to a node and equidistant to the other nodes
3378 
3379 
3380  N(0, 0) = 0.4;
3381  N(0, 1) = 0.2;
3382  N(0, 2) = 0.2;
3383  N(0, 3) = 0.2;
3384 
3385  N(1, 0) = 0.2;
3386  N(1, 1) = 0.4;
3387  N(1, 2) = 0.2;
3388  N(1, 3) = 0.2;
3389 
3390  N(2, 0) = 0.2;
3391  N(2, 1) = 0.2;
3392  N(2, 2) = 0.4;
3393  N(2, 3) = 0.2;
3394 
3395  N(3, 0) = 0.2;
3396  N(3, 1) = 0.2;
3397  N(3, 2) = 0.2;
3398  N(3, 3) = 0.4;
3399 
3400  pos=ZeroMatrix(4,3);
3401  for (unsigned int i=0; i!=4; i++) //going through the 4 nodes
3402  {
3403  array_1d<double, 3 > & coordinates = geom[i].Coordinates();
3404  for (unsigned int j=0; j!=4; j++) //going through the 4 particles
3405  {
3406  for (unsigned int k=0; k!=3; k++) //x,y,z
3407  pos(j,k) += N(j,i) * coordinates[k];
3408  }
3409  }
3410 
3411  }
3412 
3413 
3414 
3415  void ComputeGaussPointPositions_45(Geometry< Node >& geom, BoundedMatrix<double, 45, 3 > & pos,BoundedMatrix<double, 45, 3 > & N)
3416  {
3417  //std::cout << "NEW ELEMENT" << std::endl;
3418  unsigned int counter=0;
3419  for (unsigned int i=0; i!=9;i++)
3420  {
3421  for (unsigned int j=0; j!=(9-i);j++)
3422  {
3423  N(counter,0)=0.05+double(i)*0.1;
3424  N(counter,1)=0.05+double(j)*0.1;
3425  N(counter,2)=1.0 - ( N(counter,1)+ N(counter,0) ) ;
3426  pos(counter, 0) = N(counter,0) * geom[0].X() + N(counter,1) * geom[1].X() + N(counter,2) * geom[2].X();
3427  pos(counter, 1) = N(counter,0) * geom[0].Y() + N(counter,1) * geom[1].Y() + N(counter,2) * geom[2].Y();
3428  pos(counter, 2) = N(counter,0) * geom[0].Z() + N(counter,1) * geom[1].Z() + N(counter,2) * geom[2].Z();
3429  //std::cout << N(counter,0) << " " << N(counter,1) << " " << N(counter,2) << " " << std::endl;
3430  counter++;
3431 
3432  }
3433  }
3434 
3435  }
3436 
3437  void ComputeGaussPointPositions_initial(Geometry< Node >& geom, BoundedMatrix<double, 15, 3 > & pos,BoundedMatrix<double, 15, 3 > & N) //2D
3438  {
3439  //std::cout << "NEW ELEMENT" << std::endl;
3440  unsigned int counter=0;
3441  for (unsigned int i=0; i!=5;i++)
3442  {
3443  for (unsigned int j=0; j!=(5-i);j++)
3444  {
3445  N(counter,0)=0.05+double(i)*0.2;
3446  N(counter,1)=0.05+double(j)*0.2;
3447  N(counter,2)=1.0 - ( N(counter,1)+ N(counter,0) ) ;
3448  pos(counter, 0) = N(counter,0) * geom[0].X() + N(counter,1) * geom[1].X() + N(counter,2) * geom[2].X();
3449  pos(counter, 1) = N(counter,0) * geom[0].Y() + N(counter,1) * geom[1].Y() + N(counter,2) * geom[2].Y();
3450  pos(counter, 2) = N(counter,0) * geom[0].Z() + N(counter,1) * geom[1].Z() + N(counter,2) * geom[2].Z();
3451  //std::cout << N(counter,0) << " " << N(counter,1) << " " << N(counter,2) << " " << std::endl;
3452  counter++;
3453 
3454  }
3455  }
3456 
3457  }
3458 
3459  void ComputeGaussPointPositions_initial(Geometry< Node >& geom, BoundedMatrix<double, 20, 3 > & pos,BoundedMatrix<double, 20, 4 > & N) //3D
3460  {
3461  //std::cout << "NEW ELEMENT" << std::endl;
3462  //double total;
3463  double fraction_increment;
3464  unsigned int counter=0;
3465  for (unsigned int i=0; i!=4;i++) //going to build a particle "pyramid"(tetrahedra) by layers. the first layer will be made by a triangle of 4 base X 4 height. since it is a triangle, it means it will have 10 particles
3466  {
3467  //std::cout << "inside i" << i << std::endl;
3468  for (unsigned int j=0; j!=(4-i);j++)
3469  {
3470  //std::cout << "inside j" << j << std::endl;
3471  for (unsigned int k=0; k!=(4-i-j);k++)
3472  {
3473  //std::cout << "inside k" << k << std::endl;
3474  N(counter,0)= 0.27 * ( 0.175 + double(i) ) ; //this is our "surface" in which we will build each layer, so we must construct a triangle using what's left of the shape functions total (a total of 1)
3475 
3476  //total = 1.0 - N(counter,0);
3477  fraction_increment = 0.27; //
3478 
3479  N(counter,1)=fraction_increment * (0.175 + double(j));
3480  N(counter,2)=fraction_increment * (0.175 + double(k));
3481  N(counter,3)=1.0 - ( N(counter,0)+ N(counter,1) + N(counter,2) ) ;
3482  pos(counter, 0) = N(counter,0) * geom[0].X() + N(counter,1) * geom[1].X() + N(counter,2) * geom[2].X() + N(counter,3) * geom[3].X();
3483  pos(counter, 1) = N(counter,0) * geom[0].Y() + N(counter,1) * geom[1].Y() + N(counter,2) * geom[2].Y() + N(counter,3) * geom[3].Y();
3484  pos(counter, 2) = N(counter,0) * geom[0].Z() + N(counter,1) * geom[1].Z() + N(counter,2) * geom[2].Z() + N(counter,3) * geom[3].Z();
3485  //std::cout << N(counter,0) << " " << N(counter,1) << " " << N(counter,2) << " " << std::endl;
3486  counter++;
3487  }
3488 
3489  }
3490  }
3491 
3492  }
3493 
3494 
3495  // Bubble Sort Function for Descending Order
3496  void BubbleSort(array_1d<double,7> &distances , array_1d<int,7 > &positions, unsigned int & arrange_number)
3497  {
3498  int i, j;
3499  bool flag = true; // set flag to 1 to start first pass
3500  double temp; // holding variable
3501  int temp_position;
3502  int numLength = arrange_number;
3503  for(i = 1; (i <= numLength) && flag; i++)
3504  {
3505  flag = false;
3506  for (j=0; j < (numLength -1); j++)
3507  {
3508  if (distances[j+1] < distances[j]) // descending order simply changes to >
3509  {
3510  temp = distances[j]; // swap elements
3511  distances[j] = distances[j+1];
3512  distances[j+1] = temp;
3513 
3514  temp_position = positions[j]; //swap positions
3515  positions[j] = positions[j+1];
3516  positions[j+1] = temp_position;
3517 
3518  flag = true; // indicates that a swap occurred.
3519  }
3520  }
3521  }
3522  return; //arrays are passed to functions by address; nothing is returned
3523  }
3524 
3525 
3526 
3527  void BubbleSort(array_1d<double,9> &distances , array_1d<int,9 > &positions, unsigned int & arrange_number)
3528  {
3529  int i, j;
3530  bool flag = true; // set flag to 1 to start first pass
3531  double temp; // holding variable
3532  int temp_position;
3533  int numLength = arrange_number;
3534  for(i = 1; (i <= numLength) && flag; i++)
3535  {
3536  flag = false;
3537  for (j=0; j < (numLength -1); j++)
3538  {
3539  if (distances[j+1] < distances[j]) // descending order simply changes to >
3540  {
3541  temp = distances[j]; // swap elements
3542  distances[j] = distances[j+1];
3543  distances[j+1] = temp;
3544 
3545  temp_position = positions[j]; //swap positions
3546  positions[j] = positions[j+1];
3547  positions[j+1] = temp_position;
3548 
3549  flag = true; // indicates that a swap occurred.
3550  }
3551  }
3552  }
3553  return; //arrays are passed to functions by address; nothing is returned
3554  }
3555 
3556  template<class T>
3557  bool InvertMatrix(const T& input, T& inverse)
3558  {
3559  typedef permutation_matrix<std::size_t> pmatrix;
3560 
3561  // create a working copy of the input
3562  T A(input);
3563 
3564  // create a permutation matrix for the LU-factorization
3565  pmatrix pm(A.size1());
3566 
3567  // perform LU-factorization
3568  int res = lu_factorize(A, pm);
3569  if (res != 0)
3570  return false;
3571 
3572  // create identity matrix of "inverse"
3573  inverse.assign(identity_matrix<double> (A.size1()));
3574 
3575  // backsubstitute to get the inverse
3576  lu_substitute(A, pm, inverse);
3577 
3578  return true;
3579  }
3580 
3581  bool InvertMatrix3x3(const BoundedMatrix<double, TDim+1 , TDim+1 >& A, BoundedMatrix<double, TDim+1 , TDim+1 >& result)
3582  {
3583  double determinant = +A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
3584  -A(0,1)*(A(1,0)*A(2,2)-A(1,2)*A(2,0))
3585  +A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
3586  double invdet = 1/determinant;
3587  result(0,0) = (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
3588  result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
3589  result(2,0) = (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
3590  result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
3591  result(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
3592  result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
3593  result(0,2) = (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
3594  result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
3595  result(2,2) = (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
3596 
3597  return true;
3598  }
3599 
3600 
3601  ModelPart& mr_model_part;
3602  ModelPart* mtopographic_model_part_pointer;
3603  array_1d<double, 3 > mcalculation_domain_complete_displacement;
3604  array_1d<double, 3 > mcalculation_domain_added_displacement;
3605  bool mintialized_transfer_tool;
3606  bool muse_mesh_velocity_to_convect;
3607  int m_nparticles;
3608  int mnelems;
3609  double mDENSITY_WATER;
3610  double mDENSITY_AIR;
3611 
3612  //vector<double> mareas_vector; UNUSED SO COMMENTED
3613  int max_nsubsteps;
3614  double max_substep_dt;
3615  int mmaximum_number_of_particles;
3616  std::vector< PFEM_Particle_Fluid > mparticles_vector; //Point<3>
3617  int mlast_elem_id;
3618  bool modd_timestep;
3619  bool mparticle_printing_tool_initialized;
3620  unsigned int mfilter_factor;
3621  unsigned int mlast_node_id;
3622  //ModelPart& mr_particle_model_part;
3623 
3624  vector<int> mnumber_of_particles_in_elems;
3625  vector<int> mnumber_of_particles_in_elems_aux;
3626  vector<ParticlePointerVector*> mpointers_to_particle_pointers_vectors;
3627 
3628  typename BinsObjectDynamic<Configure>::Pointer mpBinsObjectDynamic;
3629  typename BinsObjectDynamic<Configure>::Pointer mpTopographicBinsObjectDynamic;
3630 
3631 
3632  void CalculateNormal(Geometry<Node >& pGeometry, array_1d<double,3>& An );
3633 
3634  };
3635 
3636  template<>
3637  void MoveParticleUtilityPFEM2<2>::CalculateNormal(Geometry<Node >& pGeometry, array_1d<double,3>& An )
3638  {
3639  array_1d<double,2> v1;
3640  v1[0] = pGeometry[1].X() - pGeometry[0].X();
3641  v1[1] = pGeometry[1].Y() - pGeometry[0].Y();
3642 
3643  An[0] = -v1[1];
3644  An[1] = v1[0];
3645  An[2] = 0.0;
3646 
3647  //now checking orientation using the normal:
3648  const unsigned int NumNodes = 2;
3649  array_1d<double,3> nodal_normal = ZeroVector(3);
3650  for (unsigned int iNode = 0; iNode < NumNodes; ++iNode)
3651  nodal_normal += pGeometry[iNode].FastGetSolutionStepValue(NORMAL);
3652 
3653  double dot_prod = nodal_normal[0]*An[0] + nodal_normal[1]*An[1];
3654  if (dot_prod<0.0)
3655  {
3656  //std::cout << "inverting the normal" << std::endl;
3657  An *= -1.0; // inverting the direction of the normal!!!
3658  }
3659 
3660  }
3661 
3662  template<>
3663  void MoveParticleUtilityPFEM2<3>::CalculateNormal(Geometry<Node >& pGeometry, array_1d<double,3>& An )
3664  {
3665  array_1d<double,3> v1,v2;
3666  v1[0] = pGeometry[1].X() - pGeometry[0].X();
3667  v1[1] = pGeometry[1].Y() - pGeometry[0].Y();
3668  v1[2] = pGeometry[1].Z() - pGeometry[0].Z();
3669 
3670  v2[0] = pGeometry[2].X() - pGeometry[0].X();
3671  v2[1] = pGeometry[2].Y() - pGeometry[0].Y();
3672  v2[2] = pGeometry[2].Z() - pGeometry[0].Z();
3673 
3674  MathUtils<double>::CrossProduct(An,v1,v2);
3675  An *= 0.5;
3676 
3677  //now checking orientation using the normal:
3678  const unsigned int NumNodes = 3;
3679  array_1d<double,3> nodal_normal = ZeroVector(3);
3680  for (unsigned int iNode = 0; iNode < NumNodes; ++iNode)
3681  nodal_normal += pGeometry[iNode].FastGetSolutionStepValue(NORMAL);
3682 
3683  double dot_prod = nodal_normal[0]*An[0] + nodal_normal[1]*An[1] + nodal_normal[2]*An[2];
3684  if (dot_prod<0.0)
3685  {
3686  //std::cout << "inverting the normal!!" << std::endl;
3687  An *= -1.0; // inverting the direction of the normal!!!
3688  }
3689 
3690  }
3691 
3692 } // namespace Kratos.
3693 
3694 #endif // KRATOS_MOVE_PART_UTILITY_DIFF2_INCLUDED defined
array_1d.h
bins_dynamic_objects.h
cfd_variables.h
Kratos::BinsObjectDynamic
Short class definition.
Definition: bins_dynamic_objects.h:57
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::Geometry::GetValue
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometry.h:627
Kratos::Geometry::Area
virtual double Area() const
This method calculate and return area or surface area of this geometry depending to it's dimension.
Definition: geometry.h:1345
Kratos::GlobalPointersVector
This class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
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::size
size_type size() const
Definition: global_pointers_vector.h:307
Kratos::GlobalPointersVector::resize
void resize(size_type new_dim) const
Definition: global_pointers_vector.h:366
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::ConditionsBegin
ConditionIterator ConditionsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1361
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::Conditions
ConditionsContainerType & Conditions(IndexType ThisIndex=0)
Definition: model_part.h:1381
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
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::ModelPart::ElementsArray
ElementsContainerType::ContainerType & ElementsArray(IndexType ThisIndex=0)
Definition: model_part.h:1209
Kratos::MoveParticleUtilityPFEM2
Definition: move_particle_utility_pfem2.h:112
Kratos::MoveParticleUtilityPFEM2::AccelerateParticlesWithoutMovingUsingDeltaVelocity
void AccelerateParticlesWithoutMovingUsingDeltaVelocity()
Definition: move_particle_utility_pfem2.h:1469
Kratos::MoveParticleUtilityPFEM2::ResetBoundaryConditionsSlip
void ResetBoundaryConditionsSlip()
Definition: move_particle_utility_pfem2.h:755
Kratos::MoveParticleUtilityPFEM2::MountBin
void MountBin()
Definition: move_particle_utility_pfem2.h:361
Kratos::MoveParticleUtilityPFEM2::PreReseedUsingTopographicDomain
void PreReseedUsingTopographicDomain(const int minimum_number_of_particles, array_1d< double, 3 > domains_added_displacement)
Definition: move_particle_utility_pfem2.h:472
Kratos::MoveParticleUtilityPFEM2::AssignNodalVelocityUsingInletConditions
void AssignNodalVelocityUsingInletConditions(const double inlet_vel)
Definition: move_particle_utility_pfem2.h:2106
Kratos::MoveParticleUtilityPFEM2::TransferLagrangianToEulerian
void TransferLagrangianToEulerian()
Definition: move_particle_utility_pfem2.h:1067
Kratos::MoveParticleUtilityPFEM2::CopyVectorVarToPreviousTimeStep
void CopyVectorVarToPreviousTimeStep(const Variable< array_1d< double, 3 > > &OriginVariable, ModelPart::NodesContainerType &rNodes)
Definition: move_particle_utility_pfem2.h:832
Kratos::MoveParticleUtilityPFEM2::CalculateDeltaVelocity
void CalculateDeltaVelocity()
Definition: move_particle_utility_pfem2.h:807
Kratos::MoveParticleUtilityPFEM2::ContainerType
Configure::ContainerType ContainerType
Definition: move_particle_utility_pfem2.h:118
Kratos::MoveParticleUtilityPFEM2::~MoveParticleUtilityPFEM2
~MoveParticleUtilityPFEM2()
Definition: move_particle_utility_pfem2.h:358
Kratos::MoveParticleUtilityPFEM2::Configure
SpatialContainersConfigure< TDim > Configure
Definition: move_particle_utility_pfem2.h:115
Kratos::MoveParticleUtilityPFEM2::IntializeTransferTool
void IntializeTransferTool(ModelPart *topographic_model_part, array_1d< double, 3 > initial_domains_offset, bool ovewrite_particle_data)
Definition: move_particle_utility_pfem2.h:383
Kratos::MoveParticleUtilityPFEM2::ResetBoundaryConditions
void ResetBoundaryConditions(bool fully_reset_nodes)
Definition: move_particle_utility_pfem2.h:656
Kratos::MoveParticleUtilityPFEM2::ExecuteParticlesPritingToolForDroppletsOnly
void ExecuteParticlesPritingToolForDroppletsOnly(ModelPart &lagrangian_model_part, int input_filter_factor)
Definition: move_particle_utility_pfem2.h:2008
Kratos::MoveParticleUtilityPFEM2::TransferLagrangianToEulerianImp
void TransferLagrangianToEulerianImp()
Definition: move_particle_utility_pfem2.h:1251
Kratos::MoveParticleUtilityPFEM2::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(MoveParticleUtilityPFEM2)
Kratos::MoveParticleUtilityPFEM2::PointType
Configure::PointType PointType
Definition: move_particle_utility_pfem2.h:116
Kratos::MoveParticleUtilityPFEM2::IteratorType
Configure::IteratorType IteratorType
Definition: move_particle_utility_pfem2.h:120
Kratos::MoveParticleUtilityPFEM2::ExecuteParticlesPritingTool
void ExecuteParticlesPritingTool(ModelPart &lagrangian_model_part, int input_filter_factor)
Definition: move_particle_utility_pfem2.h:1948
Kratos::MoveParticleUtilityPFEM2::ParticlePointerVector
PointerVector< PFEM_Particle_Fluid, PFEM_Particle_Fluid *, std::vector< PFEM_Particle_Fluid * > > ParticlePointerVector
Definition: move_particle_utility_pfem2.h:124
Kratos::MoveParticleUtilityPFEM2::PostReseed
void PostReseed(int minimum_number_of_particles, double mass_correction_factor)
Definition: move_particle_utility_pfem2.h:1664
Kratos::MoveParticleUtilityPFEM2::CopyScalarVarToPreviousTimeStep
void CopyScalarVarToPreviousTimeStep(const Variable< double > &OriginVariable, ModelPart::NodesContainerType &rNodes)
Definition: move_particle_utility_pfem2.h:857
Kratos::MoveParticleUtilityPFEM2::RotateParticlesAndDomainVelocities
void RotateParticlesAndDomainVelocities(array_1d< double, 3 > rotations)
Definition: move_particle_utility_pfem2.h:2148
Kratos::MoveParticleUtilityPFEM2::CalculateVelOverElemSize
void CalculateVelOverElemSize()
Definition: move_particle_utility_pfem2.h:591
Kratos::MoveParticleUtilityPFEM2::ResultContainerType
Configure::ResultContainerType ResultContainerType
Definition: move_particle_utility_pfem2.h:121
Kratos::MoveParticleUtilityPFEM2::PreReseed
void PreReseed(int minimum_number_of_particles)
Definition: move_particle_utility_pfem2.h:1548
Kratos::MoveParticleUtilityPFEM2::MoveParticleUtilityPFEM2
MoveParticleUtilityPFEM2(ModelPart &model_part, int maximum_number_of_particles)
Definition: move_particle_utility_pfem2.h:134
Kratos::MoveParticleUtilityPFEM2::MoveParticles
void MoveParticles(const bool discriminate_streamlines)
Definition: move_particle_utility_pfem2.h:884
Kratos::MoveParticleUtilityPFEM2::ResultIteratorType
Configure::ResultIteratorType ResultIteratorType
Definition: move_particle_utility_pfem2.h:123
Kratos::MoveParticleUtilityPFEM2::AddUniqueWeakPointer
void AddUniqueWeakPointer(GlobalPointersVector< TDataType > &v, const typename TDataType::WeakPointer candidate)
Definition: move_particle_utility_pfem2.h:1530
Kratos::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::PFEM_Particle_Fluid
PFEM Particle class.
Definition: pfem_particle_fluidonly.h:103
Kratos::PFEM_Particle_Fluid::GetEraseFlag
bool & GetEraseFlag()
Definition: pfem_particle_fluidonly.h:191
Kratos::PFEM_Particle_Fluid::GetVelocity
array_1d< float, 3 > & GetVelocity()
Definition: pfem_particle_fluidonly.h:151
Kratos::PFEM_Particle_Fluid::GetDistance
float & GetDistance()
Definition: pfem_particle_fluidonly.h:161
Kratos::Point
Point class.
Definition: point.h:59
Kratos::Point::Coordinates
CoordinatesArrayType const & Coordinates() const
Definition: point.h:215
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::PointerVector
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::SpatialContainersConfigure
Thhis class is a container for spatial search.
Definition: spatial_containers_configure.h:59
Kratos::SpatialContainersConfigure::IteratorType
typename ContainerType::iterator IteratorType
Definition: spatial_containers_configure.h:82
Kratos::SpatialContainersConfigure::ResultIteratorType
typename ResultContainerType::iterator ResultIteratorType
Definition: spatial_containers_configure.h:85
Kratos::SpatialContainersConfigure::ContainerType
typename PointerVectorSet< TEntity, IndexedObject >::ContainerType ContainerType
Container definition.
Definition: spatial_containers_configure.h:80
Kratos::SpatialContainersConfigure::ResultContainerType
typename PointerVectorSet< TEntity, IndexedObject >::ContainerType ResultContainerType
Definition: spatial_containers_configure.h:83
Kratos::Variable
Variable class contains all information needed to store and retrive data from a data container.
Definition: variable.h:63
Kratos::array_1d< double, 3 >
data_value_container.h
KRATOS_THROW_ERROR
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
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
deprecated_variables.h
enrichment_utilities.h
global_pointer_variables.h
KRATOS_INFO
#define KRATOS_INFO(label)
Definition: logger.h:250
line_2d_2.h
math_utils.h
model_part.h
CSharpKratosWrapper::ModelPart
Kratos::ModelPart ModelPart
Definition: kratos_wrapper.h:31
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::bool
REACTION_CHECK_STIFFNESS_FACTOR INNER_LOOP_ITERATION DISTANCE_THRESHOLD ACTIVE_CHECK_FACTOR AUXILIAR_COORDINATES NORMAL_GAP WEIGHTED_GAP WEIGHTED_SCALAR_RESIDUAL bool
Definition: contact_structural_mechanics_application_variables.h:93
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::SizeType
std::size_t SizeType
The definition of the size type.
Definition: mortar_classes.h:43
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::int
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
Kratos::double
TABLE_NUMBER_ANGULAR_VELOCITY TABLE_NUMBER_MOMENT I33 BEAM_INERTIA_ROT_UNIT_LENGHT_Y KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BEAM_INERTIA_ROT_UNIT_LENGHT_Z) typedef std double
Definition: DEM_application_variables.h:182
PecletTest.velocity
float velocity
Definition: PecletTest.py:54
edgebased_PureConvection.ProcessInfo
ProcessInfo
Definition: edgebased_PureConvection.py:116
edgebased_PureConvection.dist
float dist
Definition: edgebased_PureConvection.py:89
face_heat.model_part
model_part
Definition: face_heat.py:14
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_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
isotropic_damage_automatic_differentiation.threshold
threshold
Definition: isotropic_damage_automatic_differentiation.py:135
lagrangian_droplet_test.lagrangian_model_part
lagrangian_model_part
Definition: lagrangian_droplet_test.py:45
ode_solve.d
int d
Definition: ode_solve.py:397
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
read_stl.node_id
int node_id
Definition: read_stl.py:12
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
rotatingcone_PureConvectionBenchmarking.delta_t
float delta_t
Definition: rotatingcone_PureConvectionBenchmarking.py:129
script_THERMAL_CORRECT.counter
int counter
Definition: script_THERMAL_CORRECT.py:218
script.inlet_vel
inlet_vel
Definition: script.py:181
sensitivityMatrix.A
A
Definition: sensitivityMatrix.py:70
sensitivityMatrix.N
N
Definition: sensitivityMatrix.py:29
tensors::i
integer i
Definition: TensorModule.f:17
tensors::l
integer l
Definition: TensorModule.f:17
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
openmp_utils.h
pfem_2_application_variables.h
pfem_particle_fluidonly.h
spatial_containers.h
spatial_containers_configure.h
ResultIteratorType
Configure::ResultIteratorType ResultIteratorType
Definition: transfer_utility.h:252
triangle_2d_3.h
triangle_3d_3.h
Generated on Sat Feb 3 2024 02:06:39 for KratosMultiphysics by doxygen 1.9.1