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.
volume_averaging_tool.h
Go to the documentation of this file.
1 //
2 // Project Name: Kratos
3 // Last Modified by: $Author: Guillermo Casas $
4 // Date: $Date: 2012-03-08 08:56:42 $
5 //
6 //
7 
8 #if !defined(KRATOS_VOLUME_AVERAGING_TOOL )
9 #define KRATOS_VOLUME_AVERAGING_TOOL
10 
11 // /* External includes */
12 
13 // System includes
14 #include <string>
15 #include <iostream>
16 #include <cstdlib>
17 
18 
19 // Project includes
20 #include "includes/define.h"
21 #include "includes/model_part.h"
22 #include "geometries/triangle_2d_3.h"
23 #include "utilities/timer.h"
24 
25 // #include "geometries/tetrahedra_3d_4.h"
26 
27 
28 //Database includes
29 #include "spatial_containers/spatial_containers.h"
30 #include "utilities/binbased_fast_point_locator.h"
31 #include "utilities/binbased_nodes_in_element_locator.h"
32 
33 // #include "custom_utilities/drag_utilities.h"
34 
35 namespace Kratos
36 {
37 
40 
44 
48 
52 
56 
58 
62 //class VolumeAveragingTool
63 template<std::size_t TDim >
64 class VolumeAveragingTool
65 {
66 public:
69 
71  KRATOS_CLASS_POINTER_DEFINITION(VolumeAveragingTool<TDim >);
72 
76 
78  VolumeAveragingTool() {} //
79 
81  virtual ~VolumeAveragingTool() {}
82 
83 
87 
88  void CalculateSimpleCellBasedPorosity(
89  ModelPart& rFluid_ModelPart ,
90  ModelPart& rDEM_ModelPart,
91  BinBasedFastPointLocator<TDim>& bin_of_objects_fluid
92  )
93  {
94 
95 
96 
97 
98  }
99 
103 
104 
105  //**********************************************************************
106  //**********************************************************************
107 
111 
112 
116 
117 
121 
123  virtual std::string Info() const
124  {
125  return "";
126  }
127 
129  virtual void PrintInfo(std::ostream& rOStream) const {}
130 
132  virtual void PrintData(std::ostream& rOStream) const {}
133 
134 
138 
140 
141 protected:
144 
145 
149 
150 
154 
155 
159 
160 
164 
165 
169 
170 
174 
175 
177 
178 private:
181 
182 
186 
187 
188 
189  inline void CalculateCenterAndSearchRadius(Geometry<Node >&geom,
190  double& xc, double& yc, double& zc, double& R, array_1d<double,3>& N
191  )
192  {
193  double x0 = geom[0].X();
194  double y0 = geom[0].Y();
195  double x1 = geom[1].X();
196  double y1 = geom[1].Y();
197  double x2 = geom[2].X();
198  double y2 = geom[2].Y();
199 
200 
201  xc = 0.3333333333333333333*(x0+x1+x2);
202  yc = 0.3333333333333333333*(y0+y1+y2);
203  zc = 0.0;
204 
205  double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0);
206  double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1);
207  double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2);
208 
209  R = R1;
210  if (R2 > R) R = R2;
211  if (R3 > R) R = R3;
212 
213  R = 1.01 * sqrt(R);
214  }
215  //***************************************
216  //***************************************
217  inline void CalculateCenterAndSearchRadius(Geometry<Node >&geom,
218  double& xc, double& yc, double& zc, double& R, array_1d<double,4>& N
219 
220  )
221  {
222  double x0 = geom[0].X();
223  double y0 = geom[0].Y();
224  double z0 = geom[0].Z();
225  double x1 = geom[1].X();
226  double y1 = geom[1].Y();
227  double z1 = geom[1].Z();
228  double x2 = geom[2].X();
229  double y2 = geom[2].Y();
230  double z2 = geom[2].Z();
231  double x3 = geom[3].X();
232  double y3 = geom[3].Y();
233  double z3 = geom[3].Z();
234 
235 
236  xc = 0.25*(x0+x1+x2+x3);
237  yc = 0.25*(y0+y1+y2+y3);
238  zc = 0.25*(z0+z1+z2+z3);
239 
240  double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0) + (zc-z0)*(zc-z0);
241  double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1) + (zc-z1)*(zc-z1);
242  double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2) + (zc-z2)*(zc-z2);
243  double R4 = (xc-x3)*(xc-x3) + (yc-y3)*(yc-y3) + (zc-z3)*(zc-z3);
244 
245  R = R1;
246  if (R2 > R) R = R2;
247  if (R3 > R) R = R3;
248  if (R4 > R) R = R4;
249 
250  R = sqrt(R);
251  }
252  //***************************************
253  //***************************************
254  inline double CalculateVol( const double x0, const double y0,
255  const double x1, const double y1,
256  const double x2, const double y2
257  )
258  {
259  return 0.5*( (x1-x0)*(y2-y0)- (y1-y0)*(x2-x0) );
260  }
261  //***************************************
262  //***************************************
263  inline double CalculateVol( const double x0, const double y0, const double z0,
264  const double x1, const double y1, const double z1,
265  const double x2, const double y2, const double z2,
266  const double x3, const double y3, const double z3
267  )
268  {
269  double x10 = x1 - x0;
270  double y10 = y1 - y0;
271  double z10 = z1 - z0;
272 
273  double x20 = x2 - x0;
274  double y20 = y2 - y0;
275  double z20 = z2 - z0;
276 
277  double x30 = x3 - x0;
278  double y30 = y3 - y0;
279  double z30 = z3 - z0;
280 
281  double detJ = x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30;
282  return detJ*0.1666666666666666666667;
283 
284  //return 0.5*( (x1-x0)*(y2-y0)- (y1-y0)*(x2-x0) );
285  }
286  //***************************************
287  //***************************************
288  inline bool CalculatePosition( Geometry<Node >&geom,
289  const double xc, const double yc, const double zc,
290  array_1d<double,3>& N
291  )
292  {
293  double x0 = geom[0].X();
294  double y0 = geom[0].Y();
295  double x1 = geom[1].X();
296  double y1 = geom[1].Y();
297  double x2 = geom[2].X();
298  double y2 = geom[2].Y();
299 
300  double area = CalculateVol(x0,y0,x1,y1,x2,y2);
301  double inv_area = 0.0;
302  if (std::abs(area) < std::numeric_limits<double>::epsilon()) {
303 
304  //The interpolated node will not be inside an elemente with zero area
305  return false;
306  }
307  else {
308  inv_area = 1.0 / area;
309  }
310 
311  N[0] = CalculateVol(x1,y1,x2,y2,xc,yc) * inv_area;
312  N[1] = CalculateVol(x2,y2,x0,y0,xc,yc) * inv_area;
313  N[2] = CalculateVol(x0,y0,x1,y1,xc,yc) * inv_area;
314 
315 
316  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
317  return true;
318 
319  return false;
320  }
321 
322  //***************************************
323  //***************************************
324 
325  inline bool CalculatePosition( Geometry<Node >&geom,
326  const double xc, const double yc, const double zc,
327  array_1d<double,4>& N
328  )
329  {
330 
331  double x0 = geom[0].X();
332  double y0 = geom[0].Y();
333  double z0 = geom[0].Z();
334  double x1 = geom[1].X();
335  double y1 = geom[1].Y();
336  double z1 = geom[1].Z();
337  double x2 = geom[2].X();
338  double y2 = geom[2].Y();
339  double z2 = geom[2].Z();
340  double x3 = geom[3].X();
341  double y3 = geom[3].Y();
342  double z3 = geom[3].Z();
343 
344  double vol = CalculateVol(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3);
345 
346  double inv_vol = 0.0;
347  if (vol < 0.0000000000001)
348  {
349 
350  //The interpolated node will not be inside an elemente with zero volume
351  return false;
352  }
353  else
354  {
355  inv_vol = 1.0 / vol;
356  }
357 
358  N[0] = CalculateVol(x1,y1,z1,x3,y3,z3,x2,y2,z2,xc,yc,zc) * inv_vol;
359  N[1] = CalculateVol(x3,y3,z3,x0,y0,z0,x2,y2,z2,xc,yc,zc) * inv_vol;
360  N[2] = CalculateVol(x3,y3,z3,x1,y1,z1,x0,y0,z0,xc,yc,zc) * inv_vol;
361  N[3] = CalculateVol(x0,y0,z0,x1,y1,z1,x2,y2,z2,xc,yc,zc) * inv_vol;
362 
363 
364  if (N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[3] >=0.0 &&
365  N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0 && N[3] <=1.0)
366  //if the xc yc zc is inside the tetrahedron return true
367  return true;
368 
369  return false;
370  }
371 //el_it Element iterator
372 //N Shape functions
373 //step_data_size
374 //pnode pointer to the node
375 
376 
377  //projecting an array1D 2Dversion
378  void Interpolate(
379  Element::Pointer el_it,
380  const array_1d<double,3>& N,
381  Node::Pointer pnode,
382  Variable<array_1d<double,3> >& rOriginVariable,
383  Variable<array_1d<double,3> >& rDestinationVariable)
384  {
385 
386  //Geometry element of the rOrigin_ModelPart
387  Geometry< Node >& geom = el_it->GetGeometry();
388 
389  //unsigned int buffer_size = pnode->GetBufferSize();
390 
391  //for (unsigned int step = 0; step < buffer_size; step++)
392  //{
393  //getting the data of the solution step
394  array_1d<double,3>& step_data = (pnode)->FastGetSolutionStepValue(rDestinationVariable, 0);
395  //const array_1d<double,3>& velocity = (pnode)->FastGetSolutionStepValue(VELOCITY, 0);
396  const array_1d<double,3>& node0_data = geom[0].FastGetSolutionStepValue(rOriginVariable, 0);
397  const array_1d<double,3>& node1_data = geom[1].FastGetSolutionStepValue(rOriginVariable, 0);
398  const array_1d<double,3>& node2_data = geom[2].FastGetSolutionStepValue(rOriginVariable, 0);
399 
400  //copying this data in the position of the vector we are interested in
401  for (unsigned int j= 0; j< TDim; j++)
402  {
403  step_data[j] = N[0] * node0_data[j] + N[1] * node1_data[j] + N[2] * node2_data[j];
404 
405  }
406 
407 
408  //}
409 // pnode->GetValue(IS_VISITED) = 1.0;
410 
411  }
412 
413  //projecting an array1D 3Dversion
414  void Interpolate(
415  Element::Pointer el_it,
416  const array_1d<double,4>& N,
417  Node::Pointer pnode,
418  Variable<array_1d<double,3> >& rOriginVariable,
419  Variable<array_1d<double,3> >& rDestinationVariable) {
420 
421  //Geometry element of the rOrigin_ModelPart
422  Geometry< Node >& geom = el_it->GetGeometry();
423 
424  //getting the data of the solution step
425  array_1d<double, 3 > & step_data = (pnode)->FastGetSolutionStepValue(rDestinationVariable, 0);
426  const array_1d<double, 3 > & node0_data = geom[0].FastGetSolutionStepValue(rOriginVariable, 0);
427  const array_1d<double, 3 > & node1_data = geom[1].FastGetSolutionStepValue(rOriginVariable, 0);
428  const array_1d<double, 3 > & node2_data = geom[2].FastGetSolutionStepValue(rOriginVariable, 0);
429  const array_1d<double, 3 > & node3_data = geom[3].FastGetSolutionStepValue(rOriginVariable, 0);
430 
431  //copying this data in the position of the vector we are interested in
432  for (unsigned int j = 0; j < TDim; j++) {
433  step_data[j] = N[0] * node0_data[j] + N[1] * node1_data[j] + N[2] * node2_data[j] + N[3] * node3_data[j];
434  }
435  // pnode->GetValue(IS_VISITED) = 1.0;
436 
437  }
438  //projecting a scalar 2Dversion
439  void Interpolate(
440  Element::Pointer el_it,
441  const array_1d<double,3>& N,
442  Node::Pointer pnode,
443  Variable<double>& rOriginVariable,
444  Variable<double>& rDestinationVariable) {
445 
446  //Geometry element of the rOrigin_ModelPart
447  Geometry< Node >& geom = el_it->GetGeometry();
448 
449 
450  //getting the data of the solution step
451  double& step_data = (pnode)->FastGetSolutionStepValue(rDestinationVariable, 0);
452  const double node0_data = geom[0].FastGetSolutionStepValue(rOriginVariable, 0);
453  const double node1_data = geom[1].FastGetSolutionStepValue(rOriginVariable, 0);
454  const double node2_data = geom[2].FastGetSolutionStepValue(rOriginVariable, 0);
455 
456  //copying this data in the position of the vector we are interested in
457 
458  step_data = N[0] * node0_data + N[1] * node1_data + N[2] * node2_data;
459 
460 
461  // // pnode->GetValue(IS_VISITED) = 1.0;
462 
463  }
464  //projecting a scalar 3Dversion
465  void Interpolate(
466  Element::Pointer el_it,
467  const array_1d<double,4>& N,
468  Node::Pointer pnode,
469  Variable<double>& rOriginVariable,
470  Variable<double>& rDestinationVariable)
471 {
472  //Geometry element of the rOrigin_ModelPart
473  Geometry< Node >& geom = el_it->GetGeometry();
474 
475 
476 
477  //getting the data of the solution step
478  double& step_data = (pnode)->FastGetSolutionStepValue(rDestinationVariable, 0);
479  const double node0_data = geom[0].FastGetSolutionStepValue(rOriginVariable, 0);
480  const double node1_data = geom[1].FastGetSolutionStepValue(rOriginVariable, 0);
481  const double node2_data = geom[2].FastGetSolutionStepValue(rOriginVariable, 0);
482  const double node3_data = geom[3].FastGetSolutionStepValue(rOriginVariable, 0);
483 
484  //copying this data in the position of the vector we are interested in
485 
486  step_data = N[0] * node0_data + N[1] * node1_data + N[2] * node2_data + N[3] * node3_data;
487 
488  // pnode->GetValue(IS_VISITED) = 1.0;
489 
490  }
491  //2Dversion
492  void Transfer(
493  Element::Pointer el_it,
494  const array_1d<double,3>& N,
495  Node::Pointer pnode,
496  Variable<array_1d<double,3> >& rDestinationVariable,
497  Variable<array_1d<double,3> >& rOriginVariable,
498  int n_particles_per_depth_distance) {
499 
500  //Geometry element of the rOrigin_ModelPart
501  Geometry< Node >& geom = el_it->GetGeometry();
502 
503 
504  //getting the data of the solution step
505  const array_1d<double, 3 > & step_data = (pnode)->FastGetSolutionStepValue(rOriginVariable, 0);
506  array_1d<double, 3 > & node0_data = geom[0].FastGetSolutionStepValue(rDestinationVariable, 0);
507  array_1d<double, 3 > & node1_data = geom[1].FastGetSolutionStepValue(rDestinationVariable, 0);
508  array_1d<double, 3 > & node2_data = geom[2].FastGetSolutionStepValue(rDestinationVariable, 0);
509  double origin_data;
510  //copying this data in the position of the vector we are interested in
511  for (unsigned int j = 0; j < TDim; j++) {
512  origin_data = step_data[j];
513  node0_data[j] += -N[0] * origin_data * n_particles_per_depth_distance;
514  node1_data[j] += -N[1] * origin_data * n_particles_per_depth_distance;
515  node2_data[j] += -N[2] * origin_data * n_particles_per_depth_distance;
516  }
517 
518 // pnode->GetValue(IS_VISITED) = 1.0;
519 
520  }
521 
522  //3Dversion
523  void Transfer(
524  Element::Pointer el_it,
525  const array_1d<double,4>& N,
526  Node::Pointer pnode,
527  Variable<array_1d<double,3> >& rOriginVariable,
528  Variable<array_1d<double,3> >& rDestinationVariable)
529 
530  {
531 
532  //Geometry element of the rOrigin_ModelPart
533  Geometry< Node >& geom = el_it->GetGeometry();
534 
535  //unsigned int buffer_size = pnode->GetBufferSize();
536 
537  //for (unsigned int step = 0; step<buffer_size; step++)
538  //{
539  //getting the data of the solution step
540  const array_1d<double,3>& step_data = (pnode)->FastGetSolutionStepValue(rDestinationVariable, 0);
541  array_1d<double,3>& node0_data = geom[0].FastGetSolutionStepValue(rOriginVariable, 0);
542  array_1d<double,3>& node1_data = geom[1].FastGetSolutionStepValue(rOriginVariable, 0);
543  array_1d<double,3>& node2_data = geom[2].FastGetSolutionStepValue(rOriginVariable, 0);
544  array_1d<double,3>& node3_data = geom[3].FastGetSolutionStepValue(rOriginVariable, 0);
545  double origin_data;
546  //copying this data in the position of the vector we are interested in
547  for (unsigned int j= 0; j< TDim; j++)
548  {
549  origin_data = step_data[j];
550  node0_data[j] += -N[0] * origin_data;
551  node1_data[j] += -N[1] * origin_data;
552  node2_data[j] += -N[2] * origin_data;
553  node3_data[j] += -N[3] * origin_data;
554  }
555  //}
556 // pnode->GetValue(IS_VISITED) = 1.0;
557 
558  }
559  inline void Clear(ModelPart::NodesContainerType::iterator node_it, int step_data_size )
560  {
561  unsigned int buffer_size = node_it->GetBufferSize();
562 
563  for (unsigned int step = 0; step<buffer_size; step++)
564  {
565  //getting the data of the solution step
566  double* step_data = (node_it)->SolutionStepData().Data(step);
567 
568  //copying this data in the position of the vector we are interested in
569  for (int j= 0; j< step_data_size; j++)
570  {
571  step_data[j] = 0.0;
572  }
573  }
574 
575  }
576 
577  inline void ClearVariables(ModelPart::NodesContainerType::iterator node_it , Variable<array_1d<double,3> >& rVariable)
578  {
579  array_1d<double, 3>& Aux_var = node_it->FastGetSolutionStepValue(rVariable, 0);
580 
581  noalias(Aux_var) = ZeroVector(3);
582 
583  }
584 
585 
586  inline void ClearVariables(ModelPart::NodesContainerType::iterator node_it, Variable<double>& rVariable)
587  {
588  double& Aux_var = node_it->FastGetSolutionStepValue(rVariable, 0);
589 
590  Aux_var = 0.0;
591 
592  }
593 
594 
595 
596 
600 
604 
605 
609 
610 
614 
615 
619 
621  VolumeAveragingTool& operator=(VolumeAveragingTool const& rOther);
622 
623 
625 
626 }; // Class VolumeAveragingTool
627 
629 
632 
633 
637 
638 
639 
640 
642 template<std::size_t TDim>
643 inline std::ostream& operator << (std::ostream& rOStream,
644  const VolumeAveragingTool<TDim>& rThis)
645 {
646  rThis.PrintInfo(rOStream);
647  rOStream << std::endl;
648  rThis.PrintData(rOStream);
649 
650  return rOStream;
651 }
653 
654 
655 } // namespace Kratos.
656 
657 #endif // KRATOS_VOLUME_AVERAGING_TOOL defined
658 
659 
binbased_fast_point_locator.h
binbased_nodes_in_element_locator.h
Kratos::BinBasedFastPointLocator
This class is designed to allow the fast location of MANY points on the top of a 3D mesh.
Definition: binbased_fast_point_locator.h:68
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::VolumeAveragingTool
This class is designed to manage the volume averaging of DEM properties and their projection onto a m...
Definition: volume_averaging_tool.h:65
Kratos::VolumeAveragingTool::PrintData
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: volume_averaging_tool.h:132
Kratos::VolumeAveragingTool::CalculateSimpleCellBasedPorosity
void CalculateSimpleCellBasedPorosity(ModelPart &rFluid_ModelPart, ModelPart &rDEM_ModelPart, BinBasedFastPointLocator< TDim > &bin_of_objects_fluid)
Definition: volume_averaging_tool.h:88
Kratos::VolumeAveragingTool::VolumeAveragingTool
VolumeAveragingTool()
Default constructor.
Definition: volume_averaging_tool.h:78
Kratos::VolumeAveragingTool::Info
virtual std::string Info() const
Turn back information as a stemplate<class T, std::size_t dim> tring.
Definition: volume_averaging_tool.h:123
Kratos::VolumeAveragingTool::PrintInfo
virtual void PrintInfo(std::ostream &rOStream) const
Print information about this object.
Definition: volume_averaging_tool.h:129
Kratos::VolumeAveragingTool::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(VolumeAveragingTool< TDim >)
Pointer definition of VolumeAveragingTool.
Kratos::VolumeAveragingTool::~VolumeAveragingTool
virtual ~VolumeAveragingTool()
Destructor.
Definition: volume_averaging_tool.h:81
Kratos::array_1d< double, 3 >
model_part.h
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::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::operator<<
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
face_heat.step
int step
Definition: face_heat.py:88
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.R
R
Definition: isotropic_damage_automatic_differentiation.py:172
quadrature.j
int j
Definition: quadrature.py:648
rotating_cone.xc
float xc
Definition: rotating_cone.py:77
rotating_cone.yc
float yc
Definition: rotating_cone.py:78
sensitivityMatrix.N
N
Definition: sensitivityMatrix.py:29
spatial_containers.h
triangle_2d_3.h
Generated on Sat Feb 3 2024 02:26:33 for KratosMultiphysics by doxygen 1.9.1