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.
laplacian_smoothing.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosDelaunayMeshingApplication $
3 // Created by: $Author: JMCarbonell $
4 // Last modified by: $Co-Author: $
5 // Date: $Date: April 2018 $
6 // Revision: $Revision: 0.0 $
7 //
8 //
9 
10 #if !defined(KRATOS_LAPLACIAN_SMOOTHING_H_INCLUDED)
11 #define KRATOS_LAPLACIAN_SMOOTHING_H_INCLUDED
12 
13 // System includes
14 
15 // Project includes
16 #include "includes/variables.h"
17 #include "spatial_containers/spatial_containers.h"
18 #include "custom_utilities/mesher_utilities.hpp"
19 
20 #include "delaunay_meshing_application_variables.h"
21 
22 #ifdef SINGLE
23 #define REAL float
24 #else // not SINGLE
25 #define REAL double
26 #endif // not SINGLE
27 
28 #include "triangle.h"
29 
30 // External includes
31 #if !defined(KRATOS_TETGEN_EXTERNAL_H_INCLUDED)
32 #define KRATOS_TETGEN_EXTERNAL_H_INCLUDED
33 #include "tetgen.h"
34 #endif
35 
36 
37 //VARIABLES used:
38 //Data:
39 //StepData: DISPLACEMENT, CONTACT_FORCE, NORMAL, OFFSET
40 //Flags: (checked) BOUNDARY, TO_ERASE, INSIDE
41 // (set)
42 // (modified) INSIDE
43 // (reset)
44 
45 namespace Kratos
46 {
49 
53 
57 
61 
65 
67 
70  class LaplacianSmoothing
71  {
72  public:
73 
76 
78  KRATOS_CLASS_POINTER_DEFINITION( LaplacianSmoothing );
79 
80  typedef ModelPart::PropertiesType PropertiesType;
81  typedef ModelPart::MeshType MeshType;
82  typedef ModelPart::ElementsContainerType ElementsContainerType;
83  typedef ModelPart::NodesContainerType NodesContainerType;
84  typedef ModelPart::MeshType::GeometryType::PointsArrayType PointsArrayType;
85 
89 
91  LaplacianSmoothing(ModelPart& rModelPart)
92  {
93  for (unsigned int i = 0; i < rModelPart.NumberOfMeshes(); ++i)
94  mPreviousMeshes.push_back(Kratos::make_shared<MeshType>((rModelPart.GetMesh(i)).Clone())); //in fact do not clones
95 
96  } //
97 
99  virtual ~LaplacianSmoothing() {}
100 
101 
105 
106 
110 
111 
112  //*******************************************************************************************
113  //*******************************************************************************************
114 
115  void ApplyMeshSmoothing(ModelPart& rModelPart,
116  std::vector<int> & PreservedElements,
117  const int* pElementsList,
118  const int& NumberOfPoints)
119  {
120 
121  KRATOS_TRY
122 
123  NodesContainerType& rNodes = rModelPart.Nodes();
124 
125  ModelPart::ElementsContainerType::iterator element_begin = rModelPart.ElementsBegin();
126 
127  const unsigned int nds = element_begin->GetGeometry().size();
128 
129  //*******************************************************************
130  //NEIGHBOUR NODES:
131 
132  std::vector<int> EmptyVector(0);
133  std::vector<std::vector<int> > NeigbourNodesList(rNodes.size());
134  std::fill( NeigbourNodesList.begin(), NeigbourNodesList.end(), EmptyVector );
135 
136  this->GetNeigbourNodes(NeigbourNodesList, PreservedElements, pElementsList, NumberOfPoints, nds);
137 
138 
139  //*******************************************************************
140  //MOVE NODES: LAPLACIAN SMOOTHING:
141 
142 
143  double convergence_tol = 0.001;
144  double smoothing_factor = 0.1;
145  double smoothing_iters = 3;//4
146  double iters = 0;
147 
148  bool simple = true; //weight = 1;
149 
150  bool converged = false;
151  double MaxLength = 0;
152  double NewMaxLength = 0;
153 
154  bool contact_active = false;
155  double boundary_weight = 0.9; //(0,1]
156  double contact_weight = 0.8; //(0,1]
157 
158  unsigned int number_of_nodes = 0;
159 
160  ModelPart::NodesContainerType::iterator nodes_begin = rModelPart.NodesBegin();
161 
162  while ( iters<smoothing_iters && converged==false ){
163 
164  //std::cout<<" Iter "<< iters <<std::endl;
165 
166  array_1d<double,3> P;
167  array_1d<double,3> Q;//neighbour position
168  array_1d<double,3> D;
169 
170 
171  double TotalWeight = 0;
172  double Weight = 0;
173  array_1d<double,3> TotalDistance;
174 
175 
176  //convergence variables
177  double Length = 0;
178  MaxLength = NewMaxLength;
179  NewMaxLength = 0;
180 
181 
182  number_of_nodes = 0;
183  for(unsigned int in = 0; in<rNodes.size(); in++)
184  {
185 
186  unsigned int NumberOfNeighbours = NeigbourNodesList[in+1].size();
187 
188  if(rNodes[in+1].IsNot(BOUNDARY) && rNodes[in+1].IsNot(RIGID) && rNodes[in+1].IsNot(TO_ERASE) && NumberOfNeighbours>1)
189  {
190 
191  TotalDistance.clear();
192  TotalWeight = 0;
193  Weight = 0;
194 
195  //point position
196  P[0] = (nodes_begin+in)->X();
197  P[1] = (nodes_begin+in)->Y();
198  P[2] = (nodes_begin+in)->Z();
199 
200  //std::cout<<" Initial Position: "<<P<<std::endl;
201  Length = 0;
202 
203  for(unsigned int i = 0; i < NumberOfNeighbours; ++i)
204  {
205  //neighbour position
206  Q[0] = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->X();
207  Q[1] = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->Y();
208  Q[2] = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->Z();
209 
210  D = P-Q;
211 
212  Length =sqrt(D[0]*D[0]+D[1]*D[1]+D[2]*D[2]);
213 
214 
215  if( simple ){
216 
217  Weight = 1;
218 
219  }
220  else{
221 
222  if(Length !=0)
223  Weight = ( 1.0/Length );
224  else
225  Weight = 0;
226  }
227 
228  if( (nodes_begin+(NeigbourNodesList[in+1][i]-1))->Is(BOUNDARY) ){
229 
230  contact_active = false;
231  if( (nodes_begin+(NeigbourNodesList[in+1][i]-1))->SolutionStepsDataHas(CONTACT_FORCE) ){
232  array_1d<double, 3 > & ContactForce = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->FastGetSolutionStepValue(CONTACT_FORCE);
233  if( norm_2(ContactForce) !=0 )
234  contact_active = true;
235  }
236 
237  if( contact_active )
238  Weight *= contact_weight;
239  else
240  Weight *= boundary_weight;
241 
242  }
243 
244  if(NewMaxLength<Length)
245  NewMaxLength = Length;
246 
247  TotalDistance += (Weight*(Q-P)) ;
248  TotalWeight += Weight ;
249 
250  }
251 
252 
253  if(TotalWeight!=0)
254  D = ( smoothing_factor / TotalWeight ) * TotalDistance;
255  else
256  D.clear();
257 
258 
259  P += D;
260 
261  (nodes_begin+in)->X() = P[0];
262  (nodes_begin+in)->Y() = P[1];
263  (nodes_begin+in)->Z() = P[2];
264 
265  number_of_nodes +=1;
266  }
267 
268  }
269 
270 
271  if( (NewMaxLength-MaxLength)/NewMaxLength < convergence_tol ){
272  converged = true;
273  if( GetEchoLevel() > 0 )
274  std::cout<<" Laplacian smoothing convergence achieved "<<std::endl;
275  }
276 
277 
278  iters++;
279 
280  }
281 
282  if(iters==smoothing_iters && !converged){
283  if( GetEchoLevel() > 0 )
284  std::cout<<" WARNING: Laplacian smoothing convergence NOT achieved (iters:"<<iters<<")"<<std::endl;
285  }
286 
287  //*******************************************************************
288  //MOVE NODES: BOUNDARY SMOOTHING
289  SetBoundarySmoothing (rModelPart, PreservedElements, pElementsList, NumberOfPoints);
290 
291 
292  //*******************************************************************
293  //MOVE NODES: BOUNDARY PROJECTION
294  //SetInsideProjection (rModelPart, out, NeigbourNodesList);
295 
296  //*******************************************************************
297  //TRANSFER VARIABLES TO NEW NODES POSITION:
298  ProjectVariablesToNewPositions(rModelPart);
299 
300  KRATOS_CATCH( "" )
301 
302  }
303 
304 
305  //*******************************************************************************************
306  //*******************************************************************************************
307 
308  void ProjectVariablesToNewPositions(ModelPart& rModelPart)
309  {
310 
311  KRATOS_TRY
312 
313  bool transfer=true; //transfer active or inactive
314 
315  if(transfer==true){
316 
317  //create the list of the nodes to be check during the search (new positions after the laplacian smoothing)
318 
319  std::vector<Node::Pointer> list_of_nodes;
320 
321  for(NodesContainerType::iterator i_node = rModelPart.NodesBegin() ; i_node != rModelPart.NodesEnd() ; i_node++)
322  {
323 
324  //std::cout<<" ID: "<<i_node->Id()<<" NODAL_H "<<i_node->FastGetSolutionStepValue(NODAL_H)<<std::endl;
325 
326  //if(rNodes[in+1].IsNot(BOUNDARY) && rNodes[in+1].IsNot(TO_ERASE) && NumberOfNeighbours>1)
327  if( i_node->IsNot(TO_ERASE) ){
328  (list_of_nodes).push_back(*(i_node.base()));
329  }
330  // else {
331  // std::cout <<" LLM:PSEUDOERROR node to erase : " << i_node->Id() << std::endl;
332  // }
333  }
334 
335 
336  //defintions for spatial search
337  typedef Node PointType;
338  typedef Node::Pointer PointPointerType;
339  typedef std::vector<PointPointerType> PointVector;
340  typedef PointVector::iterator PointIterator;
341  //typedef std::vector<double> DistanceVector;
342  typedef std::vector<double>::iterator DistanceIterator;
343 
344  typedef Bucket<3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType;
345  typedef Tree< KDTreePartition<BucketType> > KdtreeType; //Kdtree
346  //defintions for spatial search
347 
348  unsigned int bucket_size = 20;
349  KdtreeType NodesTree(list_of_nodes.begin(),list_of_nodes.end(),bucket_size);
350 
351  //Find out where the new nodes belong to:
352  unsigned int number_of_nodes = list_of_nodes.size()+1;
353 
354  ModelPart::ElementsContainerType::iterator element_begin = rModelPart.ElementsBegin();
355  const unsigned int nds = element_begin->GetGeometry().size();
356 
357  std::vector<double> ShapeFunctionsN(nds);
358 
359 
360  if( number_of_nodes < rModelPart.Nodes().size()+1 ){
361  number_of_nodes = rModelPart.Nodes().size()+1;
362  std::cout<<" WARNING: ISOLATED NODES DELETED "<<std::endl;
363  }
364 
365 
366  std::vector<VariablesListDataValueContainer> VariablesListVector(number_of_nodes);
367  std::vector<int> UniquePosition (number_of_nodes);
368  std::fill( UniquePosition.begin(), UniquePosition.end(), 0 );
369 
370  //unsigned int step_data_size = rModelPart.GetNodalSolutionStepDataSize();
371  VariablesList& rVariablesList = rModelPart.GetNodalSolutionStepVariablesList();
372 
373  //find the center and "radius" of the element
374  double radius = 0;
375  Node center(0,0.0,0.0,0.0);
376 
377  unsigned int MaximumNumberOfPointsInRadius = list_of_nodes.size();
378  std::vector<Node::Pointer> PointsInRadius (MaximumNumberOfPointsInRadius);
379  std::vector<double> PointsInRadiusDistances (MaximumNumberOfPointsInRadius);
380 
381  //geometry
382  std::vector<std::vector<double> > ElementPointCoordinates(nds);
383  std::vector<double> PointCoordinates(3); //dimension=3
384  std::fill( PointCoordinates.begin(), PointCoordinates.end(), 0.0 );
385  std::fill( ElementPointCoordinates.begin(), ElementPointCoordinates.end(), PointCoordinates );
386 
387 
388  for(ModelPart::ElementsContainerType::const_iterator ie = rModelPart.ElementsBegin();
389  ie != rModelPart.ElementsEnd(); ie++)
390  {
391 
392  for(unsigned int i=0; i<nds; ++i)
393  {
394  //ID[cn] = rElementsList[el][cn].Id();
395  const array_1d<double,3>& Displacement = ie->GetGeometry()[i].FastGetSolutionStepValue(DISPLACEMENT);
396  PointCoordinates[0] = ie->GetGeometry()[i].X0() + Displacement[0];
397  PointCoordinates[1] = ie->GetGeometry()[i].Y0() + Displacement[1];
398  PointCoordinates[2] = ie->GetGeometry()[i].Z0() + Displacement[2];
399 
400  ElementPointCoordinates[i] = PointCoordinates;
401  }
402 
403  std::fill( PointCoordinates.begin(), PointCoordinates.end(), 0.0 );
404  MeshDataTransferUtilities DataTransferUtilities;
405  DataTransferUtilities.CalculateCenterAndSearchRadius( ElementPointCoordinates, PointCoordinates, radius );
406 
407  //find all of the new nodes within the radius
408  center.X() = PointCoordinates[0];
409  center.Y() = PointCoordinates[1];
410  center.Z() = PointCoordinates[2];
411 
412  double Radius = radius * 1.15;
413  int NumberOfPointsInRadius = NodesTree.SearchInRadius (center, Radius, PointsInRadius.begin(), PointsInRadiusDistances.begin(), MaximumNumberOfPointsInRadius);
414 
415  //std::cout<<"[ID:"<<ie->Id()<<"]: NumberOfPointsInRadius "<<NumberOfPointsInRadius<<" Radius "<<Radius<<std::endl;
416 
417  //check if inside and eventually interpolate
418  for(std::vector<Node::Pointer>::iterator it_found = PointsInRadius.begin(); it_found != (PointsInRadius.begin() + NumberOfPointsInRadius) ; ++it_found)
419  {
420 
421  if((*it_found)->IsNot(TO_ERASE)){
422 
423  //std::cout<<" Found ID "<<(*it_found)->Id()<<std::endl;
424 
425  PointCoordinates[0] = (*it_found)->X();
426  PointCoordinates[1] = (*it_found)->Y();
427  PointCoordinates[2] = (*it_found)->Z();
428 
429  bool is_inside = MesherUtilities::CalculatePosition( ElementPointCoordinates, PointCoordinates, ShapeFunctionsN );
430 
431  if(is_inside == true)
432  {
433  //std::cout<<" Node interpolation: "<<(*it_found)->Id()<<" VariablesList size "<<VariablesListVector.size()<<std::endl;
434  //std::cout<<" Node interpolation: "<<(*it_found)->Id()<<" N ["<<ie->GetGeometry()[0].Id()<<"] "<<ShapeFunctionsN[0]<<" ["<<ie->GetGeometry()[1].Id()<<"] "<<ShapeFunctionsN[1]<<" ["<<ie->GetGeometry()[2].Id()<<"] "<<ShapeFunctionsN[2]<<std::endl;
435  if(UniquePosition [(*it_found)->Id()] == 0){
436 
437  UniquePosition [(*it_found)->Id()] = 1;
438 
439  double alpha = 0.25; //[0,1] //smoothing level of the interpolation
440 
441  MeshDataTransferUtilities DataTransferUtilities;
442  VariablesListVector[(*it_found)->Id()] = DataTransferUtilities.InterpolateVariables( ie->GetGeometry(), ShapeFunctionsN, rVariablesList, (*it_found), alpha );
443 
444  }
445  else{
446  UniquePosition [(*it_found)->Id()] += 1;
447  //std::cout<<" Node "<<(*it_found)->Id()<<" is relocated again in a element:: "<<ie->Id()<<" num locations "<<UniquePosition [(*it_found)->Id()]<<std::endl;
448  //std::cout<<" ShapeFunctionsN "<<ShapeFunctionsN.size()<<std::endl;
449  }
450 
451  }
452  }
453  }
454 
455  }
456 
457 
458 
459  //the search moves the nodes order using its PreIds
460  rModelPart.Nodes().Sort();
461 
462 
463  //*******************************************************************
464  //CREATE NEW NODE INFORMATION:
465 
466  int id=0;
467  for(NodesContainerType::iterator i_node = rModelPart.NodesBegin() ; i_node != rModelPart.NodesEnd() ; ++i_node)
468  {
469 
470  if( UniquePosition[i_node->Id()] && i_node->IsNot(TO_ERASE) ){
471 
472  if ( i_node->SolutionStepsDataHas(DISPLACEMENT) == false)
473  {
474  std::cout << " WEIRD " << std::endl;
475  std::cout << " Laplacian. ThisNode Does not have displacemenet " << i_node->Id() << std::endl;
476  std::cout << " X: " << i_node->X() << " Y: " << i_node->Y() << " Z: " << i_node->Z() << std::endl;
477  }
478 
479  if( i_node->IsNot(BOUNDARY) ){
480 
481  //recover the original position of the node
482  id = i_node->Id();
483 
484  i_node->SolutionStepData() = VariablesListVector[id];
485 
486  const array_1d<double,3>& disp = i_node->FastGetSolutionStepValue(DISPLACEMENT);
487 
488  i_node->X0() = i_node->X() - disp[0];
489  i_node->Y0() = i_node->Y() - disp[1];
490  i_node->Z0() = i_node->Z() - disp[2];
491 
492  }
493  else if ( i_node->Is(BOUNDARY) && i_node->Is(INSIDE) ){ // Set the position of boundary laplacian
494 
495  i_node->Set(INSIDE, false);
496 
497  //recover the original position of the node
498  id = i_node->Id();
499 
500  i_node->SolutionStepData() = VariablesListVector[id];
501 
502  const array_1d<double,3>& disp = i_node->FastGetSolutionStepValue(DISPLACEMENT);
503 
504  bool MoveFixedNodes = false;
505  if (MoveFixedNodes)
506  {
507  i_node->X0() = i_node->X() - disp[0];
508  i_node->Y0() = i_node->Y() - disp[1];
509  i_node->Z0() = i_node->Z() - disp[2];
510  }
511  else {
512  if ( i_node->pGetDof(DISPLACEMENT_X)->IsFixed() == false) {
513  i_node->X0() = i_node->X() - disp[0];
514  }
515  if ( i_node->pGetDof(DISPLACEMENT_Y)->IsFixed() == false) {
516  i_node->Y0() = i_node->Y() - disp[1];
517  }
518 
519  if ( i_node->pGetDof(DISPLACEMENT_Z)->IsFixed() == false) {
520  i_node->Z0() = i_node->Z() - disp[2];
521  }
522  }
523 
524  }
525 
526 
527  }
528  else{
529 
530 
531  if ( i_node->SolutionStepsDataHas(DISPLACEMENT) == false)
532  {
533  std::cout << " WEIRD " << std::endl;
534  std::cout << " Laplacian. ThisNode Does not have displacemenet " << i_node->Id() << std::endl;
535  std::cout << " X: " << i_node->X() << " Y: " << i_node->Y() << " Z: " << i_node->Z() << std::endl;
536  }
537 
538  //recover the original position of the node
539  id = i_node->Id();
540 
541  const array_1d<double,3>& disp = i_node->FastGetSolutionStepValue(DISPLACEMENT);
542 
543  i_node->X0() = i_node->X() - disp[0];
544  i_node->Y0() = i_node->Y() - disp[1];
545  i_node->Z0() = i_node->Z() - disp[2];
546 
547  if( GetEchoLevel() > 1 ){
548  std::cout << " OUT::PSEUDOERROR: IN this line, there is a node that does not have displacement" << std::endl;
549  std::cout << " Laplacian. ThisNode new information does not have displacement " << i_node->Id() << std::endl;
550  std::cout << " THIS IS BECAUSE THE NODE is out of the DOMAIN and the interpolation is wrong" << std::endl;
551  std::cout << " X: " << i_node->X() << " Y: " << i_node->Y() << " Z: " << i_node->Z() << std::endl;
552  }
553 
554  }
555 
556  //std::cout<<"(B) ID: "<<i_node->Id()<<" NODAL_H "<<i_node->FastGetSolutionStepValue(NODAL_H)<<std::endl;
557  }
558 
559 
560  }
561  else{
562 
563  for(NodesContainerType::iterator i_node = rModelPart.NodesBegin() ; i_node != rModelPart.NodesEnd() ; ++i_node)
564  {
565 
566  //recover the original position of the node
567  array_1d<double,3>& disp = i_node->FastGetSolutionStepValue(DISPLACEMENT);
568  if(norm_2(disp)>0){
569  i_node->X0() = i_node->X() - disp[0];
570  i_node->Y0() = i_node->Y() - disp[1];
571  i_node->Z0() = i_node->Z() - disp[2];
572  }
573 
574  //Set the position of boundary laplacian (Reset the flag)
575  if ( i_node->Is(BOUNDARY) && i_node->IsNot(TO_ERASE) && i_node->Is(INSIDE) )
576  {
577  i_node->Set(INSIDE, false); //LMV.
578  }
579 
580  }
581 
582  }
583 
584 
585  KRATOS_CATCH( "" )
586  }
587 
588 
589 
590  //*******************************************************************************************
591  //*******************************************************************************************
592 
593  void SetBoundarySmoothing(ModelPart& rModelPart,
594  std::vector<int> & PreservedElements,
595  const int* pElementsList,
596  const int& NumberOfPoints)
597  {
598  KRATOS_TRY
599 
600  NodesContainerType& rNodes = rModelPart.Nodes();
601 
602  ModelPart::ElementsContainerType::iterator element_begin = rModelPart.ElementsBegin();
603 
604  const unsigned int nds = element_begin->GetGeometry().size();
605 
606 
607  //*******************************************************************
608  //NEIGHBOR NODES:
609 
610  std::vector<int> EmptyVector(0);
611  std::vector<std::vector<int> > NeigbourNodesList(rNodes.size());
612  std::fill( NeigbourNodesList.begin(), NeigbourNodesList.end(), EmptyVector );
613 
614  this->GetBoundaryNeigbourNodes(rNodes, NeigbourNodesList, PreservedElements, pElementsList, NumberOfPoints, nds);
615 
616  //*******************************************************************
617  //MOVE BOUNDARY NODES: LAPLACIAN SMOOTHING:
618 
619  double convergence_tol =0.001;
620  double smoothing_factor=0.1; //0.1
621  double smoothing_iters =4; //3,4
622  double iters=0;
623 
624  bool simple = true; //weight = 1;
625  bool converged = false;
626 
627  double MaxLength=0;
628  double NewMaxLength=0;
629 
630  unsigned int number_of_nodes = 0;
631 
632  ModelPart::NodesContainerType::iterator nodes_begin = rModelPart.NodesBegin();
633 
634  while ( iters<smoothing_iters && converged==false ){
635 
636  //std::cout<<" Iter "<< iters <<std::endl;
637 
638  array_1d<double,3> P;
639  array_1d<double,3> Q;//neighbour position
640  array_1d<double,3> D;
641 
642 
643  double TotalWeight = 0;
644  double Weight = 0;
645  array_1d<double,3> TotalDistance;
646 
647 
648  //convergence variables
649  double Length = 0;
650  MaxLength = NewMaxLength;
651  NewMaxLength = 0;
652 
653  for(unsigned int in = 0; in<rNodes.size(); ++in)
654  {
655  unsigned int NumberOfNeighbours = NeigbourNodesList[in+1].size();
656 
657  if(rNodes[in+1].Is(BOUNDARY) && rNodes[in+1].IsNot(TO_ERASE) && rNodes[in+1].IsNot(BOUNDARY) &&
658  rNodes[in+1].Is(INSIDE) && NumberOfNeighbours>1 )
659  {
660  TotalDistance.clear();
661  TotalWeight = 0;
662  Weight = 0;
663 
664  //point position
665  P[0] = (nodes_begin+in)->X();
666  P[1] = (nodes_begin+in)->Y();
667  P[2] = (nodes_begin+in)->Z();
668 
669 
670  //std::cout<<" Initial Position: "<<P<<std::endl;
671  Length = 0;
672 
673  for(unsigned int i = 0; i < NumberOfNeighbours; ++i)
674  {
675  //neighbour position
676  Q[0] = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->X();
677  Q[1] = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->Y();
678  Q[2] = (nodes_begin+(NeigbourNodesList[in+1][i]-1))->Z();
679 
680 
681  D = P-Q;
682 
683  Length =sqrt(D[0]*D[0]+D[1]*D[1]+D[2]*D[2]);
684 
685 
686  if( simple ){
687 
688  Weight = 1;
689 
690  }
691  else{
692 
693  if(Length !=0)
694  Weight = ( 1.0/Length );
695  else
696  Weight = 0;
697  }
698 
699  if(NewMaxLength<Length)
700  NewMaxLength = Length;
701 
702  TotalDistance += (Weight*(Q-P)) ;
703  TotalWeight += Weight ;
704 
705  }
706 
707 
708  if(TotalWeight!=0)
709  D = ( smoothing_factor / TotalWeight ) * TotalDistance;
710  else
711  D.clear();
712 
713 
714  P += D;
715 
716 
717  (nodes_begin+in)->X() = P[0];
718  (nodes_begin+in)->Y() = P[1];
719  (nodes_begin+in)->Z() = P[2];
720 
721 
722  number_of_nodes +=1;
723 
724  }
725 
726  //rNodes[in+1].Set(INSIDE,false); //LMV: Reset the flag after interpolation. Indeed, if the flag is set, only one iteration takes place
727  }
728 
729 
730  if( (NewMaxLength-MaxLength)/NewMaxLength < convergence_tol ){
731  converged = true;
732  if( GetEchoLevel() > 0 )
733  std::cout<<" Laplacian smoothing convergence achieved "<<std::endl;
734  }
735 
736 
737  iters++;
738 
739  }
740 
741  if(iters==smoothing_iters && !converged){
742  if( GetEchoLevel() > 0 )
743  std::cout<<" WARNING: Boundary Laplacian smoothing convergence NOT achieved (iters:"<<iters<<")"<<std::endl;
744  }
745 
746  KRATOS_CATCH( "" )
747  }
748 
749  //*******************************************************************************************
750  //*******************************************************************************************
751 
755  virtual void SetEchoLevel(int Level)
756  {
757  mEchoLevel = Level;
758  }
759 
760  int GetEchoLevel()
761  {
762  return mEchoLevel;
763  }
764 
768 
769 
773 
774 
778 
780  virtual std::string Info() const
781  {
782  return "";
783  }
784 
786  virtual void PrintInfo(std::ostream& rOStream) const {}
787 
789  virtual void PrintData(std::ostream& rOStream) const {}
790 
791 
795 
796 
798 
799  protected:
802 
803 
807 
808 
812 
813 
817 
818 
822 
823 
827 
828 
832 
833 
835 
836  private:
839 
840 
844 
845  ModelPart::MeshesContainerType mPreviousMeshes;
846 
847  int mEchoLevel;
848 
852 
854  LaplacianSmoothing& operator=(LaplacianSmoothing const& rOther);
855 
859 
860 
864 
865 
869 
870 
874 
875  //*******************************************************************************************
876  //*******************************************************************************************
877 
878  void GetNeigbourNodes (std::vector<std::vector<int> >& list_of_neighbor_nodes, std::vector<int> & PreservedElements,const int* ElementList, const int& NumberOfPoints, const unsigned int& nds)
879  {
880 
881  KRATOS_TRY
882 
883  if( (int)list_of_neighbor_nodes.size() != NumberOfPoints+1 ){
884  list_of_neighbor_nodes.resize(NumberOfPoints+1);
885  }
886  std::vector<int> empty_vector(0);
887  std::fill( list_of_neighbor_nodes.begin(), list_of_neighbor_nodes.end(), empty_vector );
888 
889  bool neighb_set = false;
890  int neighb_size = 0;
891 
892  for(unsigned int el = 0; el<PreservedElements.size(); ++el)
893  {
894  if(PreservedElements[el])
895  {
896  //a) Create list of node neighbors (list_of_neighbor_nodes)
897  for(unsigned int ipn=0; ipn<nds; ++ipn)
898  {
899 
900  for(unsigned int jpn=0; jpn<nds; ++jpn)
901  {
902  if(ipn!=jpn){
903  //add unique node neighbor
904  neighb_size = list_of_neighbor_nodes[ElementList[el*nds+ipn]].size();
905  neighb_set = false;
906  for(int npn=0; npn<neighb_size; ++npn)
907  {
908  if( list_of_neighbor_nodes[ElementList[el*nds+ipn]][npn]==(ElementList[el*nds+jpn]) ){
909  neighb_set=true;
910  }
911  }
912  if(neighb_set==false){
913  list_of_neighbor_nodes[ElementList[el*nds+ipn]].push_back(ElementList[el*nds+jpn]);
914  }
915  }
916  }
917  }
918  }
919  }
920 
921 
922  KRATOS_CATCH( "" )
923 
924  }
925 
926  //*******************************************************************************************
927  //*******************************************************************************************
928 
929  void GetBoundaryNeigbourNodes (NodesContainerType& rNodes, std::vector<std::vector<int> >& list_of_neighbor_nodes, std::vector<int> & PreservedElements,const int* ElementList, const int& NumberOfPoints, const unsigned int& nds)
930  {
931 
932  KRATOS_TRY
933 
934  if( (int)list_of_neighbor_nodes.size() != NumberOfPoints+1 ){
935  list_of_neighbor_nodes.resize(NumberOfPoints+1);
936  std::vector<int> empty_vector(0);
937  std::fill( list_of_neighbor_nodes.begin(), list_of_neighbor_nodes.end(), empty_vector );
938  }
939 
940  bool neighb_set = false;
941  int neighb_size = 0;
942 
943  for(unsigned int el = 0; el<PreservedElements.size(); ++el)
944  {
945  if(PreservedElements[el])
946  {
947  //a) Create list of node neighbors (list_of_neighbor_nodes)
948  for(unsigned int ipn=0; ipn<nds; ++ipn)
949  {
950  if(rNodes[ElementList[el*nds+ipn]].Is(BOUNDARY)){
951  for(unsigned int jpn=0; jpn<nds; ++jpn)
952  {
953  if(ipn!=jpn && rNodes[ElementList[el*nds+jpn]].Is(BOUNDARY)){
954  //add unique node neighbor
955  neighb_size = list_of_neighbor_nodes[ElementList[el*nds+ipn]].size();
956  neighb_set = false;
957  for(int npn=0; npn<neighb_size; ++npn)
958  {
959  if( list_of_neighbor_nodes[ElementList[el*nds+ipn]][npn]==(ElementList[el*nds+jpn]) ){
960  neighb_set=true;
961  }
962  }
963  if(neighb_set==false){
964  list_of_neighbor_nodes[ElementList[el*nds+ipn]].push_back(ElementList[el*nds+jpn]);
965  }
966  }
967  }
968  }
969  }
970  }
971  }
972 
973 
974  KRATOS_CATCH( "" )
975 
976  }
977 
978 
979  //*******************************************************************************************
980  //*******************************************************************************************
981 
982  std::vector<double> SetRanks (ModelPart& rModelPart,
983  struct triangulateio &out,
984  std::vector<std::vector<int> > & list_of_neighbor_nodes)
985  {
986 
987  KRATOS_TRY
988 
989  //set ranks
990 
991  std::vector<double> nodes_ranks;
992  nodes_ranks.resize(MesherUtilities::GetMaxNodeId(rModelPart)+1);
993  //nodes_ranks.resize(rModelPart.NumberOfNodes()+1); //mesh 0
994  std::fill( nodes_ranks.begin(), nodes_ranks.end(), 0 );
995 
996 
997  //initial assignation
998  for(int in = 0; in<out.numberofpoints; ++in)
999  {
1000  bool contact_active = false;
1001 
1002  if( rModelPart.Nodes()[in+1].SolutionStepsDataHas(CONTACT_FORCE) ){
1003  array_1d<double, 3 > & ContactForceNormal = rModelPart.Nodes()[in+1].FastGetSolutionStepValue(CONTACT_FORCE);
1004  if(norm_2(ContactForceNormal) !=0 )
1005  contact_active = true;
1006  }
1007 
1008  if( contact_active && rModelPart.Nodes()[in+1].IsNot(BOUNDARY)){
1009  nodes_ranks[in+1]=5;
1010  }
1011  // else if( norm_2(ContactForceNormal)==0 && rModelPart.Nodes()[in+1].Is(BOUNDARY)){
1012  // nodes_ranks[in+1]=1;
1013  // }
1014 
1015  }
1016 
1017 
1018  //RANK assignation:
1019  double rang_assign = 0;
1020  double rang_top = 5; //3;
1021 
1022  while (rang_assign<rang_top){
1023 
1024  for(int in = 0; in<out.numberofpoints; ++in)
1025  {
1026  if(nodes_ranks[in+1]==rang_assign){
1027 
1028  //Rank 0
1029  unsigned int shared_node=1;
1030  unsigned int NumberOfNeighbours = list_of_neighbor_nodes[in+1].size();
1031  for(unsigned int i = 0; i < NumberOfNeighbours; ++i)
1032  {
1033  shared_node = list_of_neighbor_nodes[in+1][i];
1034 
1035  if(nodes_ranks[shared_node]>rang_assign)
1036  nodes_ranks[shared_node]=rang_assign+1;
1037  }
1038  }
1039 
1040  }
1041 
1042  rang_assign++;
1043  }
1044 
1045  return nodes_ranks;
1046 
1047  KRATOS_CATCH( "" )
1048 
1049  }
1050 
1051  //*******************************************************************************************
1052  //*******************************************************************************************
1053 
1054  std::vector<double> SetFirstLayer (ModelPart& rModelPart,
1055  struct triangulateio &out,
1056  std::vector<std::vector<int> > & list_of_neighbor_nodes)
1057  {
1058  KRATOS_TRY
1059 
1060  //set ranks
1061 
1062  std::vector<double> nodes_layer;
1063  nodes_layer.resize(MesherUtilities::GetMaxNodeId(rModelPart)+1);
1064  //nodes_layer.resize(rModelPart.NumberOfNodes()+1); //mesh 0
1065  std::fill( nodes_layer.begin(), nodes_layer.end(), 0 );
1066 
1067 
1068  //initial assignation
1069  for(int in = 0; in<out.numberofpoints; ++in)
1070  {
1071  if(rModelPart.Nodes()[in+1].IsNot(BOUNDARY)){
1072  nodes_layer[in+1]=2;
1073  }
1074 
1075  }
1076 
1077 
1078  //LAYER assignation:
1079  double layer_assign = 0;
1080  double layer_top = 1;
1081 
1082  while (layer_assign<layer_top){
1083 
1084  for(int in = 0; in<out.numberofpoints; ++in)
1085  {
1086  if(nodes_layer[in+1]==layer_assign){
1087 
1088  //Rank 0
1089  unsigned int shared_node=1;
1090  unsigned int NumberOfNeighbours = list_of_neighbor_nodes[in+1].size();
1091  for(unsigned int i = 0; i < NumberOfNeighbours; ++i)
1092  {
1093  shared_node = list_of_neighbor_nodes[in+1][i];
1094 
1095  if(nodes_layer[shared_node]>layer_assign)
1096  nodes_layer[shared_node]=layer_assign+1;
1097  }
1098  }
1099 
1100  }
1101 
1102  layer_assign++;
1103  }
1104 
1105  return nodes_layer;
1106 
1107  KRATOS_CATCH( "" )
1108 
1109  }
1110 
1111 
1112  //*******************************************************************************************
1113  //*******************************************************************************************
1114 
1115  //GENERAL :: TODO
1116  //Note : to increase de robustness I propose to detect the layer nodes, via setting ranks to nodes
1117  // then ensure that the layer nodes have a certain distance to the boundary. Bigger than the offset applied
1118 
1119 
1120  void SetInsideProjection (ModelPart& rModelPart,
1121  struct triangulateio &out,
1122  std::vector<std::vector<int> > & list_of_neighbor_nodes)
1123  {
1124 
1125  KRATOS_TRY
1126 
1127  std::vector<double> nodes_ranks = SetRanks(rModelPart,out,list_of_neighbor_nodes);
1128 
1129  std::vector<double> nodes_layer = SetFirstLayer(rModelPart,out,list_of_neighbor_nodes);
1130 
1131  double movement_factor = 1.2;
1132  double contact_factor = 2.0;
1133  const array_1d<double,3> ZeroPoint(3,0.0);
1134 
1135  std::vector<array_1d<double,3> > initial_nodes_distances (MesherUtilities::GetMaxNodeId(rModelPart)+1);
1136  //std::vector<array_1d<double,3> > initial_nodes_distances (rModelPart.NumberOfNodes()+1);
1137  std::fill( initial_nodes_distances.begin(), initial_nodes_distances.end(), ZeroPoint );
1138 
1139  for(int in = 0; in<out.numberofpoints; ++in)
1140  {
1141 
1142  //if(nodes_ranks[in+1]<=1)
1143  if(nodes_ranks[in+1]<1)
1144  {
1145  array_1d<double, 3 > DeltaDisplacement = rModelPart.Nodes()[in+1].FastGetSolutionStepValue(DISPLACEMENT) - rModelPart.Nodes()[in+1].FastGetSolutionStepValue(DISPLACEMENT,1);
1146 
1147  array_1d<double, 3>& Normal= rModelPart.Nodes()[in+1].FastGetSolutionStepValue(NORMAL); //BOUNDARY_NORMAL must be set as nodal variable
1148 
1149  double projection=inner_prod(DeltaDisplacement,Normal);
1150 
1151  bool contact_active = false;
1152  if( rModelPart.Nodes()[in+1].SolutionStepsDataHas(CONTACT_FORCE) ){
1153  array_1d<double, 3 > & ContactForce = rModelPart.Nodes()[in+1].FastGetSolutionStepValue(CONTACT_FORCE);
1154  if(norm_2(ContactForce)!=0)
1155  contact_active = true;
1156  }
1157 
1158  if(contact_active){
1159  initial_nodes_distances[in+1] = (-1)*(movement_factor*contact_factor)*fabs(projection)*Normal;
1160  }
1161  else{
1162  initial_nodes_distances[in+1] = (-1)*(movement_factor)*fabs(projection)*Normal;
1163  }
1164 
1165  }
1166 
1167  //layer modification
1168 
1169  array_1d<double,3> P;
1170  array_1d<double,3> Q;//neighbour position
1171  array_1d<double,3> D;
1172  int dimension = 2;
1173 
1174 
1175  if(nodes_layer[in+1]==1)
1176  {
1177 
1178  unsigned int NumberOfNeighbours = list_of_neighbor_nodes[in+1].size();
1179 
1180  //point position
1181  P[0] = out.pointlist[in*dimension];
1182  P[1] = out.pointlist[in*dimension+1];
1183  P[2] = 0;
1184 
1185  unsigned int shared_node = 1;
1186  for(unsigned int i = 0; i < NumberOfNeighbours; ++i)
1187  {
1188  shared_node = list_of_neighbor_nodes[in+1][i];
1189 
1190  if(rModelPart.Nodes()[shared_node].Is(BOUNDARY)){
1191 
1192  //neighbour position
1193  Q[0] = out.pointlist[(list_of_neighbor_nodes[in+1][i]-1)*dimension];
1194  Q[1] = out.pointlist[(list_of_neighbor_nodes[in+1][i]-1)*dimension+1];
1195  Q[2] = 0;
1196  array_1d<double, 3>& Normal= rModelPart.Nodes()[shared_node].FastGetSolutionStepValue(NORMAL); //BOUNDARY_NORMAL must be set as nodal variable
1197 
1198  D = Q-P;
1199 
1200  double projection=inner_prod(D,Normal);
1201 
1202  array_1d<double, 3>& Offset = rModelPart.Nodes()[shared_node].FastGetSolutionStepValue(OFFSET);
1203  double offset = norm_2(Offset);
1204 
1205  double secure_offset_factor = 1.1;
1206 
1207  if(projection<offset && offset!=0)
1208  initial_nodes_distances[in+1] = (-1)*(secure_offset_factor)*(offset-projection)*Normal;
1209  }
1210 
1211  }
1212 
1213  }
1214  }
1215 
1216 
1217  double smoothing_factor=0.15;
1218  double smoothing_iters =10;
1219 
1220  double iters=0;
1221 
1222  while ( iters<smoothing_iters ){
1223 
1224  //std::cout<<" Iter "<< iters <<std::endl;
1225 
1226  double TotalWeight = 0;
1227  double Weight = 0;
1228  array_1d<double,3> TotalDistance;
1229  array_1d<double,3> Distance;
1230  array_1d<double,3> OffsetDistance;
1231 
1232  //convergence variables
1233  for(int in = 0; in<out.numberofpoints; ++in)
1234  {
1235 
1236  TotalDistance = initial_nodes_distances[in+1];
1237  OffsetDistance.clear();
1238 
1239 
1240  if(nodes_ranks[in+1]>0)
1241  {
1242 
1243  if(nodes_layer[in+1]==1)
1244  OffsetDistance = TotalDistance;
1245 
1246  unsigned int NumberOfNeighbours = list_of_neighbor_nodes[in+1].size();
1247 
1248  unsigned int shared_node=1;
1249 
1250  TotalWeight = 0;
1251  Weight = 0;
1252  Distance.clear();
1253 
1254 
1255  for(unsigned int i = 0; i < NumberOfNeighbours; ++i)
1256  {
1257 
1258  //neighbour
1259  shared_node = list_of_neighbor_nodes[in+1][i];
1260 
1261  Weight = 1.0 / (nodes_ranks[shared_node]+1.0);
1262  TotalWeight += Weight ;
1263  Distance += initial_nodes_distances[shared_node] * Weight;
1264  }
1265 
1266 
1267  if(TotalWeight!=0)
1268  Distance *= ( 1.0 / TotalWeight );
1269  else
1270  Distance = initial_nodes_distances[in+1];
1271 
1272  TotalDistance += smoothing_factor*(Distance-initial_nodes_distances[in+1]);
1273  }
1274 
1275 
1276  if( nodes_layer[in+1]==1 && norm_2(OffsetDistance)>norm_2(TotalDistance)+1e-10 ){
1277  TotalDistance = OffsetDistance;
1278  if( GetEchoLevel() > 0 )
1279  std::cout<<" Layer Correction "<<norm_2(OffsetDistance)<<" > "<<norm_2(TotalDistance)<<std::endl;
1280  }
1281 
1282  initial_nodes_distances[in+1] = TotalDistance;
1283 
1284 
1285  }
1286 
1287 
1288  iters++;
1289 
1290  }
1291 
1292  int dimension = 2;
1293  for(int in = 0; in<out.numberofpoints; ++in)
1294  {
1295 
1296  // array_1d<double, 3>& Projection= rModelPart.Nodes()[in+1].FastGetSolutionStepValue(FORCE_EXTERNAL); //BOUNDARY_NORMAL must be set as nodal variable
1297  // Projection = initial_nodes_distances[in+1];
1298 
1299  if(nodes_ranks[in+1]>0)
1300  {
1301  //std::cout<<" Projection Set "<<initial_nodes_distances[in+1]<<std::endl;
1302  out.pointlist[in*dimension] += initial_nodes_distances[in+1][0];
1303  out.pointlist[in*dimension+1] += initial_nodes_distances[in+1][1];
1304  }
1305 
1306 
1307  }
1308 
1309  KRATOS_CATCH( "" )
1310 
1311 
1312  }
1313 
1314 
1315  //*******************************************************************************************
1316  //*******************************************************************************************
1317 
1318 
1319  inline void CalculateCenterAndSearchRadius(const double x0, const double y0,
1320  const double x1, const double y1,
1321  const double x2, const double y2,
1322  double& xc, double& yc, double& R)
1323 
1324  {
1325  xc = 0.3333333333333333333*(x0+x1+x2);
1326  yc = 0.3333333333333333333*(y0+y1+y2);
1327 
1328  double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0);
1329  double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1);
1330  double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2);
1331 
1332  R = R1;
1333  if(R2 > R) R = R2;
1334  if(R3 > R) R = R3;
1335 
1336  R = sqrt(R);
1337  }
1338 
1339 
1340  //*******************************************************************************************
1341  //*******************************************************************************************
1342 
1343 
1344  inline void Clear(ModelPart::NodesContainerType::iterator node_it, int step_data_size )
1345  {
1346  unsigned int buffer_size = node_it->GetBufferSize();
1347 
1348  for(unsigned int step = 0; step<buffer_size; ++step)
1349  {
1350  //getting the data of the solution step
1351  double* step_data = (node_it)->SolutionStepData().Data(step);
1352 
1353  //copying this data in the position of the vector we are interested in
1354  for(int j= 0; j< step_data_size; ++j)
1355  {
1356  step_data[j] = 0.0;
1357  }
1358  }
1359 
1360  }
1361 
1362  //*******************************************************************************************
1363  //*******************************************************************************************
1364 
1365 
1366  inline void ClearVariables(ModelPart::NodesContainerType::iterator node_it , Variable<array_1d<double,3> >& rVariable)
1367  {
1368  array_1d<double, 3>& Aux_var = node_it->FastGetSolutionStepValue(rVariable, 0);
1369 
1370  noalias(Aux_var) = ZeroVector(3);
1371 
1372  }
1373 
1374  //*******************************************************************************************
1375  //*******************************************************************************************
1376 
1377 
1378  inline void ClearVariables(ModelPart::NodesContainerType::iterator node_it, Variable<double>& rVariable)
1379  {
1380  double& Aux_var = node_it->FastGetSolutionStepValue(rVariable, 0);
1381 
1382  Aux_var = 0.0;
1383 
1384  }
1385 
1386 
1388 
1389  }; // Class LaplacianSmoothing
1390 
1392 
1395 
1396 
1400 
1401 
1403  inline std::istream& operator >> (std::istream& rIStream,
1404  LaplacianSmoothing& rThis);
1405 
1407  inline std::ostream& operator << (std::ostream& rOStream,
1408  const LaplacianSmoothing& rThis)
1409  {
1410  rThis.PrintInfo(rOStream);
1411  rOStream << std::endl;
1412  rThis.PrintData(rOStream);
1413 
1414  return rOStream;
1415  }
1417 
1418 
1419 } // namespace Kratos.
1420 
1421 #endif // KRATOS_LAPLACIAN_SMOOTHING_H_INCLUDED defined
Kratos::Bucket
Short class definition.
Definition: bucket.h:57
Kratos::Geometry::PointsArrayType
PointerVector< TPointType > PointsArrayType
Definition: geometry.h:118
Kratos::LaplacianSmoothing
Short class definition.
Definition: laplacian_smoothing.hpp:71
Kratos::LaplacianSmoothing::SetEchoLevel
virtual void SetEchoLevel(int Level)
Definition: laplacian_smoothing.hpp:755
Kratos::LaplacianSmoothing::~LaplacianSmoothing
virtual ~LaplacianSmoothing()
Destructor.
Definition: laplacian_smoothing.hpp:99
Kratos::LaplacianSmoothing::LaplacianSmoothing
LaplacianSmoothing(ModelPart &rModelPart)
Default constructor.
Definition: laplacian_smoothing.hpp:91
Kratos::LaplacianSmoothing::PrintInfo
virtual void PrintInfo(std::ostream &rOStream) const
Print information about this object.
Definition: laplacian_smoothing.hpp:786
Kratos::LaplacianSmoothing::Info
virtual std::string Info() const
Turn back information as a string.
Definition: laplacian_smoothing.hpp:780
Kratos::LaplacianSmoothing::GetEchoLevel
int GetEchoLevel()
Definition: laplacian_smoothing.hpp:760
Kratos::LaplacianSmoothing::ElementsContainerType
ModelPart::ElementsContainerType ElementsContainerType
Definition: laplacian_smoothing.hpp:82
Kratos::LaplacianSmoothing::PropertiesType
ModelPart::PropertiesType PropertiesType
Definition: laplacian_smoothing.hpp:80
Kratos::LaplacianSmoothing::NodesContainerType
ModelPart::NodesContainerType NodesContainerType
Definition: laplacian_smoothing.hpp:83
Kratos::LaplacianSmoothing::ApplyMeshSmoothing
void ApplyMeshSmoothing(ModelPart &rModelPart, std::vector< int > &PreservedElements, const int *pElementsList, const int &NumberOfPoints)
Definition: laplacian_smoothing.hpp:115
Kratos::LaplacianSmoothing::MeshType
ModelPart::MeshType MeshType
Definition: laplacian_smoothing.hpp:81
Kratos::LaplacianSmoothing::SetBoundarySmoothing
void SetBoundarySmoothing(ModelPart &rModelPart, std::vector< int > &PreservedElements, const int *pElementsList, const int &NumberOfPoints)
Definition: laplacian_smoothing.hpp:593
Kratos::LaplacianSmoothing::PointsArrayType
ModelPart::MeshType::GeometryType::PointsArrayType PointsArrayType
Definition: laplacian_smoothing.hpp:84
Kratos::LaplacianSmoothing::PrintData
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: laplacian_smoothing.hpp:789
Kratos::LaplacianSmoothing::ProjectVariablesToNewPositions
void ProjectVariablesToNewPositions(ModelPart &rModelPart)
Definition: laplacian_smoothing.hpp:308
Kratos::LaplacianSmoothing::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(LaplacianSmoothing)
Pointer definition of data transfer.
Kratos::MeshDataTransferUtilities
Short class definition.
Definition: mesh_data_transfer_utilities.hpp:46
Kratos::MeshDataTransferUtilities::CalculateCenterAndSearchRadius
void CalculateCenterAndSearchRadius(const std::vector< std::vector< double > > &rPointCoordinates, std::vector< double > &rCenter, double &rRadius)
Definition: mesh_data_transfer_utilities.hpp:334
Kratos::MeshDataTransferUtilities::InterpolateVariables
VariablesListDataValueContainer InterpolateVariables(Geometry< Node > &geom, const std::vector< double > &N, VariablesList &rVariablesList, Node::Pointer pnode, double alpha=1.0)
Definition: mesh_data_transfer_utilities.cpp:1374
Kratos::MesherUtilities::GetMaxNodeId
static unsigned int GetMaxNodeId(ModelPart &rModelPart)
Definition: mesher_utilities.hpp:1408
Kratos::MesherUtilities::CalculatePosition
static bool CalculatePosition(const std::vector< std::vector< double >> &rPointCoordinates, const std::vector< double > &rCenter, std::vector< double > &rShapeFunctionsN)
Definition: mesher_utilities.hpp:1272
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::ElementsContainerType
MeshType::ElementsContainerType ElementsContainerType
Element container. A vector set of Elements with their Id's as key.
Definition: model_part.h:168
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::NumberOfMeshes
SizeType NumberOfMeshes()
Definition: model_part.h:1776
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::GetMesh
MeshType & GetMesh(IndexType ThisIndex=0)
Definition: model_part.h:1791
Kratos::ModelPart::Nodes
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
Kratos::ModelPart::GetNodalSolutionStepVariablesList
VariablesList & GetNodalSolutionStepVariablesList()
Definition: model_part.h:549
Kratos::Node
This class defines the node.
Definition: node.h:65
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::push_back
void push_back(const TPointerType &x)
Definition: pointer_vector.h:270
Kratos::Properties
Properties encapsulates data shared by different Elements or Conditions. It can store any type of dat...
Definition: properties.h:69
Kratos::Tree
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
Kratos::VariablesList
Holds a list of variables and their position in VariablesListDataValueContainer.
Definition: variables_list.h:50
Kratos::array_1d< double, 3 >
Kratos::array_1d::clear
BOOST_UBLAS_INLINE void clear()
Definition: array_1d.h:325
Kratos::array_1d::size
BOOST_UBLAS_INLINE size_type size() const
Definition: array_1d.h:370
Kratos::array_1d::begin
BOOST_UBLAS_INLINE const_iterator begin() const
Definition: array_1d.h:606
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
delaunay_meshing_application_variables.h
mesher_utilities.hpp
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::norm_2
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
Kratos::inner_prod
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
Kratos::operator>>
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
Kratos::DistanceIterator
DistanceVector::iterator DistanceIterator
Definition: internal_variables_interpolation_process.h:56
Kratos::PointVector
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::PointIterator
PointVector::iterator PointIterator
Definition: internal_variables_interpolation_process.h:54
Kratos::operator<<
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
Kratos::PointType
Point PointType
Geometric definitions.
Definition: mortar_classes.h:36
face_heat.step
int step
Definition: face_heat.py:88
generate_convection_diffusion_explicit_element.alpha
alpha
Definition: generate_convection_diffusion_explicit_element.py:113
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.out
out
Definition: isotropic_damage_automatic_differentiation.py:200
isotropic_damage_automatic_differentiation.dimension
int dimension
Definition: isotropic_damage_automatic_differentiation.py:123
isotropic_damage_automatic_differentiation.R
R
Definition: isotropic_damage_automatic_differentiation.py:172
isotropic_damage_automatic_differentiation.Q
tuple Q
Definition: isotropic_damage_automatic_differentiation.py:235
mesh_to_mdpa_converter.radius
float radius
Definition: mesh_to_mdpa_converter.py:18
quadrature.j
int j
Definition: quadrature.py:648
read_stl.el
el
Definition: read_stl.py:25
rotating_cone.xc
float xc
Definition: rotating_cone.py:77
rotating_cone.yc
float yc
Definition: rotating_cone.py:78
tensors::i
integer i
Definition: TensorModule.f:17
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
spatial_containers.h
Generated on Sat Feb 3 2024 01:59:57 for KratosMultiphysics by doxygen 1.9.1