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.
residual_based_bossak_displacement_rotation_scheme.hpp
Go to the documentation of this file.
1 //
2 // Project Name: KratosSolidMechanicsApplication $
3 // Created by: $Author: JMCarbonell $
4 // Last modified by: $Co-Author: $
5 // Date: $Date: October 2016 $
6 // Revision: $Revision: 0.0 $
7 //
8 //
9 
10 #if !defined(KRATOS_RESIDUAL_BASED_BOSSAK_DISPLACEMENT_ROTATION_SCHEME)
11 #define KRATOS_RESIDUAL_BASED_BOSSAK_DISPLACEMENT_ROTATION_SCHEME
12 
13 /* System includes */
14 
15 /* External includes */
16 #include "boost/smart_ptr.hpp"
17 
18 /* Project includes */
19 #include "includes/define.h"
20 #include "includes/model_part.h"
21 #include "solving_strategies/schemes/scheme.h"
22 #include "includes/variables.h"
23 #include "containers/array_1d.h"
24 #include "includes/element.h"
25 #include "utilities/beam_math_utilities.hpp"
26 
27 #include "solid_mechanics_application_variables.h"
28 
29 namespace Kratos
30 {
47 // Covariant implicit time stepping algorithm: the classical Newmark algorithm of nonlinear elastodynamics and a canonical extension of the Newmark formulas to the orthogonal group SO(3) for the rotational part.
48 
49 template<class TSparseSpace, class TDenseSpace >
50 class ResidualBasedBossakDisplacementRotationScheme: public Scheme<TSparseSpace,TDenseSpace>
51 {
52 protected:
53 
54 
55  struct GeneralAlphaMethod
56  {
57  double f; //alpha Hilbert
58  double m; //alpha Bosssak
59 
60  };
61 
62  struct NewmarkMethod
63  {
64 
65  double beta;
66  double gamma;
67  double deltatime;
68 
69  //system constants
70  double c0;
71  double c1;
72  double c2;
73  double c3;
74  double c4;
75  double c5;
76  double c6;
77  double c7;
78 
79  //static-dynamic parameter
80  double static_dynamic;
81 
82  };
83 
84 
85  struct GeneralMatrices
86  {
87  std::vector< Matrix > M; //first derivative matrix (usually mass matrix)
88  std::vector< Matrix > D; //second derivative matrix (usually damping matrix)
89 
90  };
91 
92  struct GeneralVectors
93  {
94  std::vector< Vector > fm; //fist derivatives vector
95  std::vector< Vector > fd; //second derivative vector
96  };
97 
98 
99 
100 public:
101 
102 
106  KRATOS_CLASS_POINTER_DEFINITION( ResidualBasedBossakDisplacementRotationScheme );
107 
108  typedef Scheme<TSparseSpace,TDenseSpace> BaseType;
109 
110  typedef typename BaseType::TDataType TDataType;
111 
112  typedef typename BaseType::DofsArrayType DofsArrayType;
113 
114  typedef typename Element::DofsVectorType DofsVectorType;
115 
116  typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
117 
118  typedef typename BaseType::TSystemVectorType TSystemVectorType;
119 
120  typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
121 
122  typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
123 
124  typedef ModelPart::ElementsContainerType ElementsArrayType;
125 
126  typedef ModelPart::ConditionsContainerType ConditionsArrayType;
127 
128  typedef BeamMathUtils<double> BeamMathUtilsType;
129 
130  typedef Quaternion<double> QuaternionType;
131 
132 
139  ResidualBasedBossakDisplacementRotationScheme(double rDynamic = 1, double mAlphaM = 0.0)
140  :Scheme<TSparseSpace,TDenseSpace>()
141  {
142  //For pure stable Newmark Scheme
143  // mNewmark.beta= 0.25;
144  // mNewmark.gamma= 0.5;
145 
146  //For Bossak Scheme
147  mAlpha.f= 0;
148  mAlpha.m= mAlphaM; //0.0 to -0.3
149 
150  mNewmark.beta= (1.0+mAlpha.f-mAlpha.m)*(1.0+mAlpha.f-mAlpha.m)*0.25;
151  mNewmark.gamma= 0.5+mAlpha.f-mAlpha.m;
152 
153 
154  mNewmark.static_dynamic= rDynamic;
155 
156 
157  //std::cout << " MECHANICAL SCHEME: The Newmark Simo Time Integration Scheme [ beta= "<<mNewmark.beta<<" gamma= "<<mNewmark.gamma<<"]"<<std::endl;
158 
159 
160  //Allocate auxiliary memory
161  int NumThreads = OpenMPUtils::GetNumThreads();
162 
163  mMatrix.M.resize(NumThreads);
164  mMatrix.D.resize(NumThreads);
165 
166  mVector.fm.resize(NumThreads);
167  mVector.fd.resize(NumThreads);
168 
169  }
170 
173  virtual ~ResidualBasedBossakDisplacementRotationScheme
174  () {}
175 
188  //***************************************************************************
197  void Update(
198  ModelPart& r_model_part,
199  DofsArrayType& rDofSet,
200  TSystemMatrixType& A,
201  TSystemVectorType& Dx,
202  TSystemVectorType& b )
203  {
204  KRATOS_TRY
205 
206 
207  //store previous iteration rotation in previous step rotation
208  for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
209  i != r_model_part.NodesEnd(); ++i)
210  {
211  if( i->IsNot(SLAVE) && i->IsNot(RIGID) ){
212 
213 
214  //DELTA_ROTATION variable used to store previous iteration rotation
215  array_1d<double, 3 > & CurrentRotation = (i)->FastGetSolutionStepValue(ROTATION);
216  array_1d<double, 3 > & ReferenceRotation = (i)->FastGetSolutionStepValue(DELTA_ROTATION, 1);
217 
218  // Rotation at iteration i
219  ReferenceRotation = CurrentRotation;
220 
221  //STEP_DISPLACEMENT variable used to store previous iteration displacement
222  array_1d<double, 3 > & CurrentDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT);
223  array_1d<double, 3 > & ReferenceDisplacement= (i)->FastGetSolutionStepValue(STEP_DISPLACEMENT, 1);
224 
225  // Displacement at iteration i
226  ReferenceDisplacement = CurrentDisplacement;
227 
228  }
229  }
230 
231 
232  //std::cout << " Update " << std::endl;
233  //update of displacement (by DOF)
234  for (typename DofsArrayType::iterator i_dof = rDofSet.begin(); i_dof != rDofSet.end(); ++i_dof)
235  {
236  if (i_dof->IsFree() )
237  {
238  i_dof->GetSolutionStepValue() += Dx[i_dof->EquationId()];
239  }
240  }
241 
242  //LINEAR VELOCITIES AND ACCELERATIONS
243 
244  //updating time derivatives (nodally for efficiency)
245  array_1d<double, 3 > DeltaDisplacement;
246 
247  for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
248  i != r_model_part.NodesEnd(); ++i)
249  {
250 
251  if( i->IsNot(SLAVE) && i->IsNot(RIGID) ){
252 
253  array_1d<double, 3 > & CurrentDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT);
254  array_1d<double, 3 > & ReferenceDisplacement = (i)->FastGetSolutionStepValue(STEP_DISPLACEMENT, 1); //the CurrentDisplacement of the last iteration was stored here
255 
256  noalias(DeltaDisplacement) = CurrentDisplacement-ReferenceDisplacement;
257 
258  array_1d<double, 3 > & CurrentStepDisplacement = (i)->FastGetSolutionStepValue(STEP_DISPLACEMENT);
259 
260  noalias(CurrentStepDisplacement) = CurrentDisplacement - (i)->FastGetSolutionStepValue(DISPLACEMENT,1);
261 
262  //CurrentStepDisplacement += DeltaDisplacement;
263 
264  array_1d<double, 3 > & CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY);
265  array_1d<double, 3 > & CurrentAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION);
266 
267  array_1d<double, 3 > & PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY,1);
268  array_1d<double, 3 > & PreviousAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION,1);
269 
270  UpdateAcceleration ((*i), CurrentAcceleration, DeltaDisplacement, CurrentStepDisplacement, PreviousAcceleration, PreviousVelocity);
271  UpdateVelocity ((*i), CurrentVelocity, DeltaDisplacement, CurrentAcceleration, PreviousAcceleration, PreviousVelocity);
272 
273  }
274  }
275 
276  //updating time derivatives (nodally for efficiency)
277 
278  //UPDATE COMPOUND AND DELTA ROTATIONS BY QUATERNION COMPOSITION and ANGULAR VELOCITIES AND ACCELERATIONS
279  Vector TotalRotationVector = ZeroVector(3);
280  Vector StepRotationVector = ZeroVector(3);
281  Vector DeltaRotationVector = ZeroVector(3);
282 
283  Vector LinearDeltaRotationVector = ZeroVector(3);
284 
285  QuaternionType DeltaRotationQuaternion;
286  QuaternionType CurrentRotationQuaternion;
287  QuaternionType ReferenceRotationQuaternion;
288 
289 
290  for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
291  i != r_model_part.NodesEnd(); ++i)
292  {
293  if( i->IsNot(SLAVE) && i->IsNot(RIGID) ){
294 
295  // Rotation at iteration i+1
296  array_1d<double, 3 > & CurrentRotation = (i)->FastGetSolutionStepValue(ROTATION);
297  // Rotation at iteration i
298  array_1d<double, 3 > & ReferenceRotation = (i)->FastGetSolutionStepValue(DELTA_ROTATION, 1);
299 
300  //StepRotation
301  array_1d<double, 3 > & CurrentStepRotation = (i)->FastGetSolutionStepValue(STEP_ROTATION);
302 
303  //DeltaRotation
304  array_1d<double, 3 > & CurrentDeltaRotation= (i)->FastGetSolutionStepValue(DELTA_ROTATION);
305 
306  CurrentDeltaRotation = CurrentRotation-ReferenceRotation;
307 
308  // if( (i)->Id()==26 )
309  // {
310  // std::cout<<" CurrentDeltaRotation "<<CurrentDeltaRotation<<std::endl;
311  // std::cout<<" ReferenceDeltaRotation "<<ReferenceRotation<<" CurrentRotation "<<CurrentRotation<<std::endl;
312  // }
313 
314 
315  for( unsigned int j=0; j<3; j++)
316  {
317  DeltaRotationVector[j] = CurrentDeltaRotation[j]; //delta rotation computed in the iteration
318 
319  TotalRotationVector[j] = ReferenceRotation[j]; //previous iteration total rotation
320 
321  StepRotationVector[j] = CurrentStepRotation[j]; //previous iteration step rotation
322  }
323 
324 
325  DeltaRotationQuaternion = QuaternionType::FromRotationVector( DeltaRotationVector );
326 
327  // updated step rotations:
328  ReferenceRotationQuaternion = QuaternionType::FromRotationVector( StepRotationVector );
329 
330  CurrentRotationQuaternion = DeltaRotationQuaternion * ReferenceRotationQuaternion;
331 
332  CurrentRotationQuaternion.ToRotationVector( StepRotationVector );
333 
334  // ---------------
335 
336  // std::cout<<" Step Rotation B: "<<StepRotationVector<<std::endl;
337 
338  // updated compound rotations:
339  ReferenceRotationQuaternion = QuaternionType::FromRotationVector( TotalRotationVector );
340 
341  CurrentRotationQuaternion = DeltaRotationQuaternion * ReferenceRotationQuaternion;
342 
343  CurrentRotationQuaternion.ToRotationVector( TotalRotationVector );
344 
345  // ---------------
346 
347  for( unsigned int j=0; j<3; j++)
348  {
349  LinearDeltaRotationVector[j] = StepRotationVector[j] - CurrentStepRotation[j];
350  }
351 
352  // update delta rotation composed:
353  Matrix RotationMatrix;
354  (ReferenceRotationQuaternion.conjugate()).ToRotationMatrix(RotationMatrix);
355  LinearDeltaRotationVector = prod( RotationMatrix, LinearDeltaRotationVector );
356 
357  //(ReferenceRotationQuaternion.conjugate()).RotateVector3(LinearDeltaRotationVector); //exp[-X_n]*(Rotation_i+1-Rotation_i)
358  //(CurrentRotationQuaternion.conjugate()).RotateVector3(LinearDeltaRotationVector); //exp[-X_n+1]*(Rotation_i+1-Rotation_i)
359 
360  array_1d<double, 3 > LinearDeltaRotation;
361 
362  for( unsigned int j=0; j<3; j++)
363  {
364  //update rotation and step_rotation
365  CurrentRotation[j] = TotalRotationVector[j];
366  CurrentStepRotation[j] = StepRotationVector[j];
367 
368  //array_1d to pass in the acceleration and velocity updates
369  LinearDeltaRotation[j] = LinearDeltaRotationVector[j];
370  }
371 
372 
373  array_1d<double, 3 > & CurrentVelocity = (i)->FastGetSolutionStepValue(ANGULAR_VELOCITY);
374  array_1d<double, 3 > & CurrentAcceleration = (i)->FastGetSolutionStepValue(ANGULAR_ACCELERATION);
375 
376 
377  array_1d<double, 3 > & PreviousAngularAcceleration = (i)->FastGetSolutionStepValue(ANGULAR_ACCELERATION,1);
378  array_1d<double, 3 > & PreviousAngularVelocity = (i)->FastGetSolutionStepValue(ANGULAR_VELOCITY,1);
379 
380 
381  UpdateAngularAcceleration ((*i), CurrentAcceleration, LinearDeltaRotation, CurrentStepRotation, PreviousAngularAcceleration, PreviousAngularVelocity);
382  UpdateAngularVelocity ((*i), CurrentVelocity, LinearDeltaRotation, CurrentAcceleration, PreviousAngularAcceleration, PreviousAngularVelocity);
383 
384 
385  }
386 
387  }
388 
389 
390  //UPDATE SLAVE NODES
391  SlaveNodesUpdate(r_model_part);
392 
393 
394  KRATOS_CATCH( "" )
395  }
396 
397 
398 
399 
400  //***************************************************************************
401  //***************************************************************************
402 
403  void SlaveNodesUpdate(ModelPart& r_model_part)
404  {
405  KRATOS_TRY
406 
407  Matrix SkewSymVariable = ZeroMatrix(3,3);
408  Vector RadiusVector = ZeroVector(3);
409  Vector Variable = ZeroVector(3);
410  Vector AngularVariable = ZeroVector(3);
411  array_1d<double,3> VariableArray;
412 
413  array_1d<double,3> Radius;
414 
415  for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
416  i != r_model_part.NodesEnd(); ++i)
417  {
418  if( (i)->Is(SLAVE) && i->IsNot(RIGID) ){
419 
420  Element& MasterElement = *(i)->GetValue(MASTER_ELEMENTS).back();
421 
422  Node& rCenterOfGravity = MasterElement.GetGeometry()[0];
423 
424  array_1d<double, 3 >& Center = rCenterOfGravity.GetInitialPosition();
425  array_1d<double, 3 >& Displacement = rCenterOfGravity.FastGetSolutionStepValue(DISPLACEMENT);
426  array_1d<double, 3 >& Rotation = rCenterOfGravity.FastGetSolutionStepValue(ROTATION);
427  array_1d<double, 3 >& StepRotation = rCenterOfGravity.FastGetSolutionStepValue(STEP_ROTATION);
428  array_1d<double, 3 >& DeltaRotation = rCenterOfGravity.FastGetSolutionStepValue(DELTA_ROTATION);
429 
430  array_1d<double, 3 >& Velocity = rCenterOfGravity.FastGetSolutionStepValue(VELOCITY);
431  array_1d<double, 3 >& Acceleration = rCenterOfGravity.FastGetSolutionStepValue(ACCELERATION);
432  array_1d<double, 3 >& AngularVelocity = rCenterOfGravity.FastGetSolutionStepValue(ANGULAR_VELOCITY);
433  array_1d<double, 3 >& AngularAcceleration = rCenterOfGravity.FastGetSolutionStepValue(ANGULAR_ACCELERATION);
434 
435  //std::cout<<" [ MasterElement "<<MasterElement.Id() ];
436  //std::cout<<" [ Rotation:"<<Rotation<<",StepRotation:"<<StepRotation<<",DeltaRotation:"<<DeltaRotation<<"]"<<std::endl;
437  //std::cout<<" [ Velocity:"<<Velocity<<",Acceleration:"<<Acceleration<<",Displacement:"<<Displacement<<",DeltaDisplacement"<<Displacement-rCenterOfGravity.FastGetSolutionStepValue(DISPLACEMENT,1)<<"]"<<std::endl;
438 
439  //Get rotation matrix
440  QuaternionType TotalQuaternion = QuaternionType::FromRotationVector<array_1d<double,3> >(Rotation);
441 
442  Radius = (i)->GetInitialPosition() - Center;
443 
444  Matrix RotationMatrix;
445  TotalQuaternion.ToRotationMatrix(RotationMatrix);
446 
447  for(int j=0; j<3; j++)
448  RadiusVector[j] = Radius[j];
449 
450  RadiusVector = prod( RotationMatrix, RadiusVector );
451 
452  for(int j=0; j<3; j++)
453  Radius[j] = RadiusVector[j];
454 
455  //TotalQuaternion.RotateVector3<array_1d<double,3> >(Radius);
456 
457  array_1d<double, 3 >& NodeDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT);
458  array_1d<double, 3 >& NodeRotation = (i)->FastGetSolutionStepValue(ROTATION);
459  array_1d<double, 3 >& NodeStepRotation = (i)->FastGetSolutionStepValue(STEP_ROTATION);
460  array_1d<double, 3 >& NodeDeltaRotation = (i)->FastGetSolutionStepValue(DELTA_ROTATION);
461 
462  noalias(NodeDisplacement) = ( (Center + Displacement) + Radius ) - (i)->GetInitialPosition();
463  noalias(NodeRotation) = Rotation;
464  noalias(NodeStepRotation) = StepRotation;
465  noalias(NodeDeltaRotation) = DeltaRotation;
466 
467  for(int j=0; j<3; j++)
468  RadiusVector[j] = Radius[j];
469 
470  //********************
471 
472  for(int j=0; j<3; j++)
473  Variable[j] = AngularVelocity[j];
474 
475  //compute the skewsymmmetric tensor of the angular velocity
476  BeamMathUtilsType::VectorToSkewSymmetricTensor(Variable, SkewSymVariable);
477 
478  //compute the contribution of the angular velocity to the velocity v = Wxr
479  Variable = prod(SkewSymVariable,RadiusVector);
480 
481  for(int j=0; j<3; j++)
482  VariableArray[j] = Variable[j];
483 
484  (i)->FastGetSolutionStepValue(VELOCITY) = Velocity + VariableArray;
485 
486  //********************
487 
488  //centripetal acceleration:
489  for(int j=0; j<3; j++)
490  AngularVariable[j] = AngularVelocity[j];
491 
492  //compute the skewsymmmetric tensor of the angular velocity
493  BeamMathUtilsType::VectorToSkewSymmetricTensor(AngularVariable, SkewSymVariable);
494 
495  AngularVariable = prod(SkewSymVariable,Variable); //ac = Wx(Wxr)
496 
497  for(int j=0; j<3; j++)
498  Variable[j] = AngularAcceleration[j];
499 
500  //compute the skewsymmmetric tensor of the angular acceleration
501  BeamMathUtilsType::VectorToSkewSymmetricTensor(Variable, SkewSymVariable);
502 
503  //compute the contribution of the angular velocity to the velocity a = Axr
504  Variable = prod(SkewSymVariable,RadiusVector);
505 
506  for(int j=0; j<3; j++)
507  VariableArray[j] = Variable[j] + AngularVariable[j];
508 
509  (i)->FastGetSolutionStepValue(ACCELERATION) = Acceleration + VariableArray;
510 
511 
512  //********************
513  (i)->FastGetSolutionStepValue(ANGULAR_VELOCITY) = AngularVelocity;
514  (i)->FastGetSolutionStepValue(ANGULAR_ACCELERATION) = AngularAcceleration;
515 
516  // std::cout<<" [ Finalize Rigid Body Link Point : [Id:"<<(i)->Id()<<"] "<<std::endl;
517  // std::cout<<" [ Displacement:"<<NodeDisplacement<<" / StepRotation"<<NodeStepRotation<<" ] "<<std::endl;
518  // std::cout<<" [ Rotation:"<<NodeRotation<<" / Angular Acceleration"<<AngularAcceleration<<" ] "<<std::endl;
519 
520 
521  }
522 
523  }
524 
525  KRATOS_CATCH( "" )
526  }
527 
528  //***************************************************************************
529  //***************************************************************************
530 
531  //predicts the solution for the current step:
532  // x = xold + vold * Dt
533 
534  void Predict(
535  ModelPart& r_model_part,
536  DofsArrayType& rDofSet,
537  TSystemMatrixType& A,
538  TSystemVectorType& Dx,
539  TSystemVectorType& b
540  )
541  {
542 
543  KRATOS_TRY
544 
545  //std::cout << " Prediction " << std::endl;
546 
547  //double DeltaTime = r_model_part.GetProcessInfo()[DELTA_TIME];
548 
549  //DISPLACEMENTS
550 
551  for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
552  i != r_model_part.NodesEnd(); ++i)
553  {
554 
555  if( i->IsNot(SLAVE) && i->IsNot(RIGID) ){
556  //predicting displacement = PreviousDisplacement + PreviousVelocity * DeltaTime;
557  //ATTENTION::: the prediction is performed only on free nodes
558 
559 
560  array_1d<double, 3 > & CurrentDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT);
561  array_1d<double, 3 > & CurrentVelocity = (i)->FastGetSolutionStepValue(VELOCITY);
562  array_1d<double, 3 > & CurrentAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION);
563 
564  array_1d<double, 3 > & CurrentStepDisplacement = (i)->FastGetSolutionStepValue(STEP_DISPLACEMENT);
565 
566 
567  //DISPLACEMENT AND STEP DISPLACEMENT
568  unsigned int previous = 1;
569  array_1d<double, 3 > & ReferenceDisplacement = (i)->FastGetSolutionStepValue(DISPLACEMENT,previous);
570  array_1d<double, 3 > & PreviousVelocity = (i)->FastGetSolutionStepValue(VELOCITY,previous);
571  array_1d<double, 3 > & PreviousAcceleration = (i)->FastGetSolutionStepValue(ACCELERATION,previous);
572 
573  noalias(CurrentStepDisplacement) = CurrentDisplacement - ReferenceDisplacement;
574 
575  //only used in a special prediction of the displacement
576  array_1d<double, 3 > & CurrentAngularVelocity = (i)->FastGetSolutionStepValue(ANGULAR_VELOCITY);
577 
578  PredictStepDisplacement( (*i),CurrentStepDisplacement, CurrentVelocity, PreviousAcceleration, CurrentAcceleration, CurrentAngularVelocity);
579 
580  if (i->HasDofFor(LAGRANGE_MULTIPLIER_NORMAL))
581  {
582  double& PreviousLM = (i)->FastGetSolutionStepValue(LAGRANGE_MULTIPLIER_NORMAL, 1);
583  double& CurrentLM = (i)->FastGetSolutionStepValue(LAGRANGE_MULTIPLIER_NORMAL);
584 
585  CurrentLM = PreviousLM;
586  }
587 
588  //DISPLACEMENT, LINEAR VELOCITIES AND ACCELERATIONS
589 
590  PredictDisplacement(CurrentDisplacement, CurrentStepDisplacement);
591 
592  //updating time derivatives ::: please note that displacements and its time derivatives
593  //can not be consistently fixed separately
594 
595  PredictAcceleration ( (*i), CurrentAcceleration, CurrentStepDisplacement, PreviousAcceleration, PreviousVelocity);
596  PredictVelocity ( (*i), CurrentVelocity, CurrentStepDisplacement, CurrentAcceleration, PreviousAcceleration, PreviousVelocity);
597 
598  //std::cout<<" Prediction: Velocity "<<CurrentVelocity<<" Accel "<<CurrentAcceleration<<" Disp "<<CurrentDisplacement<<" StepDisp "<<CurrentStepDisplacement<<std::endl;
599  }
600 
601  }
602 
603 
604  //ROTATIONS
605 
606  for (ModelPart::NodeIterator i = r_model_part.NodesBegin();
607  i != r_model_part.NodesEnd(); ++i)
608  {
609 
610  if( i->IsNot(SLAVE) && i->IsNot(RIGID) ){
611 
612  //ATTENTION::: the prediction is performed only on free nodes
613  array_1d<double, 3 >& CurrentRotation = (i)->FastGetSolutionStepValue(ROTATION);
614  array_1d<double, 3 >& CurrentStepRotation = (i)->FastGetSolutionStepValue(STEP_ROTATION);
615  array_1d<double, 3 >& CurrentAngularVelocity = (i)->FastGetSolutionStepValue(ANGULAR_VELOCITY);
616  array_1d<double, 3 >& CurrentAngularAcceleration = (i)->FastGetSolutionStepValue(ANGULAR_ACCELERATION);
617 
618  //array_1d<double, 3 >& ImposedRotation = (i)->FastGetSolutionStepValue(IMPOSED_ROTATION);
619 
620  //ROTATIONS AND STEP ROTATIONS
621  unsigned int previous = 1;
622  array_1d<double, 3 > PreviousAngularVelocity = (i)->FastGetSolutionStepValue(ANGULAR_VELOCITY,previous);
623  array_1d<double, 3 > PreviousAngularAcceleration = (i)->FastGetSolutionStepValue(ANGULAR_ACCELERATION,previous);
624 
625  PredictStepRotation( (*i), CurrentStepRotation, CurrentAngularVelocity, PreviousAngularAcceleration, CurrentAngularAcceleration);
626 
627  //ROTATIONS, ANGULAR VELOCITIES AND ACCELERATIONS
628 
629  //updating time derivatives ::: please note that rotations and its time derivatives
630  //can not be consistently fixed separately
631 
632  PredictRotation(CurrentRotation, CurrentStepRotation);
633 
634 
635  PredictAngularAcceleration ( (*i), CurrentAngularAcceleration, CurrentStepRotation, PreviousAngularAcceleration, PreviousAngularVelocity);
636  PredictAngularVelocity ( (*i), CurrentAngularVelocity, CurrentStepRotation, CurrentAngularAcceleration, PreviousAngularAcceleration, PreviousAngularVelocity);
637 
638  }
639  }
640 
641  //UPDATE SLAVE NODES
642  SlaveNodesUpdate(r_model_part);
643 
644 
645  KRATOS_CATCH( "" )
646  }
647 
648 
649  //***************************************************************************
650  //***************************************************************************
655  void InitializeElements(ModelPart& rModelPart)
656  {
657  KRATOS_TRY
658 
659  int NumThreads = OpenMPUtils::GetNumThreads();
660  OpenMPUtils::PartitionVector ElementPartition;
661  OpenMPUtils::DivideInPartitions(rModelPart.Elements().size(), NumThreads, ElementPartition);
662 
663  #pragma omp parallel
664  {
665  int k = OpenMPUtils::ThisThread();
666  ElementsArrayType::iterator ElemBegin = rModelPart.Elements().begin() + ElementPartition[k];
667  ElementsArrayType::iterator ElemEnd = rModelPart.Elements().begin() + ElementPartition[k + 1];
668 
669  for (ElementsArrayType::iterator itElem = ElemBegin; itElem != ElemEnd; itElem++)
670  {
671  itElem->Initialize(); //function to initialize the element
672  }
673 
674  }
675 
676  this->mElementsAreInitialized = true;
677  //std::cout<<" mechanical elements are initialized "<<std::endl;
678 
679  KRATOS_CATCH( "" )
680  }
681 
682 
687  void InitializeConditions(ModelPart& rModelPart)
688  {
689  KRATOS_TRY
690 
691  if(this->mElementsAreInitialized==false)
692  KRATOS_THROW_ERROR( std::logic_error, "Before initilizing Conditions, initialize Elements FIRST", "" )
693 
694  int NumThreads = OpenMPUtils::GetNumThreads();
695 
696  OpenMPUtils::PartitionVector ConditionPartition;
697  OpenMPUtils::DivideInPartitions(rModelPart.Conditions().size(), NumThreads, ConditionPartition);
698 
699  #pragma omp parallel
700  {
701  int k = OpenMPUtils::ThisThread();
702  ConditionsArrayType::iterator CondBegin = rModelPart.Conditions().begin() + ConditionPartition[k];
703  ConditionsArrayType::iterator CondEnd = rModelPart.Conditions().begin() + ConditionPartition[k + 1];
704 
705  for (ConditionsArrayType::iterator itCond = CondBegin; itCond != CondEnd; itCond++)
706  {
707  itCond->Initialize(); //function to initialize the condition
708  }
709 
710  }
711 
712  this->mConditionsAreInitialized = true;
713  KRATOS_CATCH( "" )
714  }
715 
724  void InitializeSolutionStep(
725  ModelPart& r_model_part,
726  TSystemMatrixType& A,
727  TSystemVectorType& Dx,
728  TSystemVectorType& b)
729  {
730  KRATOS_TRY
731 
732  ProcessInfo& CurrentProcessInfo= r_model_part.GetProcessInfo();
733 
734 
735  Scheme<TSparseSpace,TDenseSpace>::InitializeSolutionStep(r_model_part,A,Dx,b);
736 
737 
738  double DeltaTime = CurrentProcessInfo[DELTA_TIME];
739 
740  if (DeltaTime == 0)
741  KRATOS_THROW_ERROR(std::logic_error, "detected delta_time = 0 in the Solution Scheme ... check if the time step is created correctly for the current model part", "" )
742 
743  //Set Newmark coefficients
744 
745  if( mNewmark.static_dynamic != 0 ){
746  CurrentProcessInfo[NEWMARK_BETA] = mNewmark.beta;
747  CurrentProcessInfo[NEWMARK_GAMMA] = mNewmark.gamma;
748  CurrentProcessInfo[BOSSAK_ALPHA] = mAlpha.m;
749  CurrentProcessInfo[COMPUTE_DYNAMIC_TANGENT] = true;
750  }
751 
752  //initializing Newmark constants
753  mNewmark.deltatime = DeltaTime;
754 
755  mNewmark.c0 = ( mNewmark.gamma / ( DeltaTime * mNewmark.beta ) );
756  mNewmark.c1 = ( 1.0 / ( DeltaTime * DeltaTime * mNewmark.beta ) );
757 
758  mNewmark.c2 = ( DeltaTime * ( 1.0 - mNewmark.gamma ) );
759  mNewmark.c3 = ( DeltaTime * mNewmark.gamma );
760  mNewmark.c4 = ( DeltaTime / mNewmark.beta );
761  mNewmark.c5 = ( DeltaTime * DeltaTime * ( 0.5 - mNewmark.beta ) / mNewmark.beta );
762 
763  //extra
764  mNewmark.c6 = ( 1.0 / (mNewmark.beta * DeltaTime) );
765  mNewmark.c7 = ( 0.5 / (mNewmark.beta) - 1.0 );
766 
767  //std::cout<<" Newmark Variables "<<mNewmark.c0<<" "<<mNewmark.c1<<" "<<mNewmark.c2<<" "<<mNewmark.c3<<" "<<mNewmark.c4<<std::endl;
768 
769 
770  KRATOS_CATCH( "" )
771  }
772  //***************************************************************************
773  //***************************************************************************
774 
779  void FinalizeSolutionStep(
780  ModelPart& rModelPart,
781  TSystemMatrixType& A,
782  TSystemVectorType& Dx,
783  TSystemVectorType& b)
784  {
785  KRATOS_TRY
786  //finalizes solution step for all of the elements
787  ElementsArrayType& rElements = rModelPart.Elements();
788  ProcessInfo& CurrentProcessInfo = rModelPart.GetProcessInfo();
789 
790  int NumThreads = OpenMPUtils::GetNumThreads();
791  OpenMPUtils::PartitionVector ElementPartition;
792  OpenMPUtils::DivideInPartitions(rElements.size(), NumThreads, ElementPartition);
793 
794  #pragma omp parallel
795  {
796  int k = OpenMPUtils::ThisThread();
797 
798  ElementsArrayType::iterator ElementsBegin = rElements.begin() + ElementPartition[k];
799  ElementsArrayType::iterator ElementsEnd = rElements.begin() + ElementPartition[k + 1];
800 
801  for (ElementsArrayType::iterator itElem = ElementsBegin; itElem != ElementsEnd; itElem++)
802  {
803  itElem->FinalizeSolutionStep(CurrentProcessInfo);
804  }
805  }
806 
807  ConditionsArrayType& rConditions = rModelPart.Conditions();
808 
809  OpenMPUtils::PartitionVector ConditionPartition;
810  OpenMPUtils::DivideInPartitions(rConditions.size(), NumThreads, ConditionPartition);
811 
812  #pragma omp parallel
813  {
814  int k = OpenMPUtils::ThisThread();
815 
816  ConditionsArrayType::iterator ConditionsBegin = rConditions.begin() + ConditionPartition[k];
817  ConditionsArrayType::iterator ConditionsEnd = rConditions.begin() + ConditionPartition[k + 1];
818 
819  for (ConditionsArrayType::iterator itCond = ConditionsBegin; itCond != ConditionsEnd; itCond++)
820  {
821  itCond->FinalizeSolutionStep(CurrentProcessInfo);
822  }
823  }
824  KRATOS_CATCH( "" )
825  }
826 
827  //***************************************************************************
828  //***************************************************************************
829 
830  void InitializeNonLinIteration(ModelPart& r_model_part,
831  TSystemMatrixType& A,
832  TSystemVectorType& Dx,
833  TSystemVectorType& b)
834  {
835  KRATOS_TRY
836  ElementsArrayType& pElements = r_model_part.Elements();
837  ProcessInfo& CurrentProcessInfo = r_model_part.GetProcessInfo();
838 
839  for (ElementsArrayType::iterator it = pElements.begin(); it != pElements.end(); ++it)
840  {
841  (it) -> InitializeNonLinearIteration(CurrentProcessInfo);
842  }
843 
844  ConditionsArrayType& pConditions = r_model_part.Conditions();
845  for (ConditionsArrayType::iterator it = pConditions.begin(); it != pConditions.end(); ++it)
846  {
847  (it) -> InitializeNonLinearIteration(CurrentProcessInfo);
848  }
849  KRATOS_CATCH( "" )
850  }
851 
852  //***************************************************************************
853  //***************************************************************************
854 
855  void InitializeNonLinearIteration(Condition::Pointer rCurrentCondition,
856  ProcessInfo& CurrentProcessInfo)
857  {
858  (rCurrentCondition) -> InitializeNonLinearIteration(CurrentProcessInfo);
859  }
860 
861 
862  //***************************************************************************
863  //***************************************************************************
864 
865  void InitializeNonLinearIteration(Element::Pointer rCurrentElement,
866  ProcessInfo& CurrentProcessInfo)
867  {
868  (rCurrentElement) -> InitializeNonLinearIteration(CurrentProcessInfo);
869  }
870 
871 
872 
873  //***************************************************************************
874  //***************************************************************************
875 
876 
877 
878 
879  //***************************************************************************
880  //***************************************************************************
881 
884  void CalculateSystemContributions(
885  Element::Pointer rCurrentElement,
886  LocalSystemMatrixType& LHS_Contribution,
887  LocalSystemVectorType& RHS_Contribution,
888  Element::EquationIdVectorType& EquationId,
889  ProcessInfo& CurrentProcessInfo)
890  {
891  KRATOS_TRY
892 
893  int thread = OpenMPUtils::ThisThread();
894 
895  //(rCurrentElement) -> InitializeNonLinearIteration(CurrentProcessInfo);
896 
897  (rCurrentElement) -> CalculateLocalSystem(LHS_Contribution,RHS_Contribution,CurrentProcessInfo);
898 
899 
900  if(mNewmark.static_dynamic !=0)
901  {
902 
903  (rCurrentElement) -> CalculateSecondDerivativesContributions(mMatrix.M[thread],mVector.fm[thread],CurrentProcessInfo);
904 
905  (rCurrentElement) -> CalculateFirstDerivativesContributions(mMatrix.D[thread],mVector.fd[thread],CurrentProcessInfo);
906 
907  }
908 
909 
910  (rCurrentElement) -> EquationIdVector(EquationId,CurrentProcessInfo);
911 
912 
913  if(mNewmark.static_dynamic !=0)
914  {
915 
916  AddDynamicsToLHS(LHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], CurrentProcessInfo);
917 
918  AddDynamicsToRHS(RHS_Contribution, mVector.fd[thread], mVector.fm[thread], CurrentProcessInfo);
919 
920  }
921 
922  //AssembleTimeSpaceLHS(rCurrentElement, LHS_Contribution, DampMatrix, MassMatrix,CurrentProcessInfo);
923 
924  KRATOS_CATCH( "" )
925  }
926 
927  void Calculate_RHS_Contribution(
928  Element::Pointer rCurrentElement,
929  LocalSystemVectorType& RHS_Contribution,
930  Element::EquationIdVectorType& EquationId,
931  ProcessInfo& CurrentProcessInfo)
932  {
933 
934  KRATOS_TRY
935 
936  int thread = OpenMPUtils::ThisThread();
937 
938  //Initializing the non linear iteration for the current element
939  //(rCurrentElement) -> InitializeNonLinearIteration(CurrentProcessInfo);
940 
941  //basic operations for the element considered
942  (rCurrentElement) -> CalculateRightHandSide(RHS_Contribution,CurrentProcessInfo);
943 
944  if(mNewmark.static_dynamic !=0)
945  {
946 
947  (rCurrentElement) -> CalculateSecondDerivativesRHS(mVector.fm[thread],CurrentProcessInfo);
948 
949  (rCurrentElement) -> CalculateFirstDerivativesRHS(mVector.fd[thread],CurrentProcessInfo);
950 
951  }
952 
953  (rCurrentElement) -> EquationIdVector(EquationId,CurrentProcessInfo);
954 
955  if(mNewmark.static_dynamic !=0)
956  {
957 
958  AddDynamicsToRHS(RHS_Contribution, mVector.fd[thread], mVector.fm[thread], CurrentProcessInfo);
959 
960  }
961 
962  KRATOS_CATCH( "" )
963 
964  }
965 
966  //***************************************************************************
967  //***************************************************************************
968 
972  void Condition_CalculateSystemContributions(
973  Condition::Pointer rCurrentCondition,
974  LocalSystemMatrixType& LHS_Contribution,
975  LocalSystemVectorType& RHS_Contribution,
976  Element::EquationIdVectorType& EquationId,
977  ProcessInfo& CurrentProcessInfo)
978  {
979 
980 
981  KRATOS_TRY
982 
983  int thread = OpenMPUtils::ThisThread();
984 
985  //Initializing the non linear iteration for the current element
986  //(rCurrentCondition) -> InitializeNonLinearIteration(CurrentProcessInfo);
987 
988  //basic operations for the element considered
989  (rCurrentCondition) -> CalculateLocalSystem(LHS_Contribution,RHS_Contribution,CurrentProcessInfo);
990 
991  if(mNewmark.static_dynamic !=0)
992  {
993 
994  (rCurrentCondition) -> CalculateSecondDerivativesContributions(mMatrix.M[thread],mVector.fm[thread],CurrentProcessInfo);
995 
996  (rCurrentCondition) -> CalculateFirstDerivativesContributions(mMatrix.D[thread],mVector.fd[thread],CurrentProcessInfo);
997 
998  }
999 
1000  (rCurrentCondition) -> EquationIdVector(EquationId,CurrentProcessInfo);
1001 
1002  if(mNewmark.static_dynamic !=0)
1003  {
1004 
1005  AddDynamicsToLHS(LHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], CurrentProcessInfo);
1006 
1007  AddDynamicsToRHS(RHS_Contribution, mVector.fd[thread], mVector.fm[thread], CurrentProcessInfo);
1008 
1009  }
1010 
1011  //AssembleTimeSpaceLHS_Condition(rCurrentCondition, LHS_Contribution,DampMatrix, MassMatrix,CurrentProcessInfo);
1012 
1013  KRATOS_CATCH( "" )
1014  }
1015 
1016 
1017  //***************************************************************************
1018  //***************************************************************************
1019 
1020  void Condition_Calculate_RHS_Contribution(
1021  Condition::Pointer rCurrentCondition,
1022  LocalSystemVectorType& RHS_Contribution,
1023  Element::EquationIdVectorType& EquationId,
1024  ProcessInfo& CurrentProcessInfo)
1025  {
1026  KRATOS_TRY
1027 
1028  int thread = OpenMPUtils::ThisThread();
1029 
1030  //Initializing the non linear iteration for the current condition
1031  //(rCurrentCondition) -> InitializeNonLinearIteration(CurrentProcessInfo);
1032 
1033  //basic operations for the element considered
1034  (rCurrentCondition) -> CalculateRightHandSide(RHS_Contribution, CurrentProcessInfo);
1035 
1036  if(mNewmark.static_dynamic !=0)
1037  {
1038 
1039  (rCurrentCondition) -> CalculateSecondDerivativesRHS(mVector.fm[thread],CurrentProcessInfo);
1040 
1041  (rCurrentCondition) -> CalculateFirstDerivativesRHS(mVector.fd[thread],CurrentProcessInfo);
1042 
1043  }
1044 
1045  (rCurrentCondition) -> EquationIdVector(EquationId, CurrentProcessInfo);
1046 
1047  //adding the dynamic contributions (static is already included)
1048 
1049  if(mNewmark.static_dynamic !=0)
1050  {
1051 
1052  AddDynamicsToRHS (RHS_Contribution, mVector.fd[thread], mVector.fm[thread], CurrentProcessInfo);
1053 
1054  }
1055 
1056  KRATOS_CATCH( "" )
1057  }
1058 
1059 
1060  //***************************************************************************
1061  //***************************************************************************
1062 
1066  void GetElementalDofList(
1067  Element::Pointer rCurrentElement,
1068  Element::DofsVectorType& ElementalDofList,
1069  ProcessInfo& CurrentProcessInfo)
1070  {
1071  rCurrentElement->GetDofList(ElementalDofList, CurrentProcessInfo);
1072  }
1073 
1074  //***************************************************************************
1075  //***************************************************************************
1076 
1080  void GetConditionDofList(
1081  Condition::Pointer rCurrentCondition,
1082  Element::DofsVectorType& ConditionDofList,
1083  ProcessInfo& CurrentProcessInfo)
1084  {
1085  rCurrentCondition->GetDofList(ConditionDofList, CurrentProcessInfo);
1086  }
1087 
1088  //***************************************************************************
1089  //***************************************************************************
1090 
1098  virtual int Check(ModelPart& r_model_part)
1099  {
1100  KRATOS_TRY
1101 
1102  int err = Scheme<TSparseSpace, TDenseSpace>::Check(r_model_part);
1103  if(err!=0) return err;
1104 
1105  //check for variables keys
1106  //verify that the variables are correctly initialized
1107  if(DISPLACEMENT.Key() == 0)
1108  KRATOS_THROW_ERROR( std::invalid_argument,"DISPLACEMENT has Key zero! (check if the application is correctly registered", "" )
1109  if(VELOCITY.Key() == 0)
1110  KRATOS_THROW_ERROR( std::invalid_argument,"VELOCITY has Key zero! (check if the application is correctly registered", "" )
1111  if(ACCELERATION.Key() == 0)
1112  KRATOS_THROW_ERROR( std::invalid_argument,"ACCELERATION has Key zero! (check if the application is correctly registered", "" )
1113 
1114  //check that variables are correctly allocated
1115  for(ModelPart::NodesContainerType::iterator it=r_model_part.NodesBegin();
1116  it!=r_model_part.NodesEnd(); it++)
1117  {
1118  if (it->SolutionStepsDataHas(DISPLACEMENT) == false)
1119  KRATOS_THROW_ERROR( std::logic_error, "DISPLACEMENT variable is not allocated for node ", it->Id() )
1120  if (it->SolutionStepsDataHas(VELOCITY) == false)
1121  KRATOS_THROW_ERROR( std::logic_error, "DISPLACEMENT variable is not allocated for node ", it->Id() )
1122  if (it->SolutionStepsDataHas(ACCELERATION) == false)
1123  KRATOS_THROW_ERROR( std::logic_error, "DISPLACEMENT variable is not allocated for node ", it->Id() )
1124  }
1125 
1126  //check that dofs exist
1127  for(ModelPart::NodesContainerType::iterator it=r_model_part.NodesBegin();
1128  it!=r_model_part.NodesEnd(); it++)
1129  {
1130  if(it->HasDofFor(DISPLACEMENT_X) == false)
1131  KRATOS_THROW_ERROR( std::invalid_argument,"missing DISPLACEMENT_X dof on node ",it->Id() )
1132  if(it->HasDofFor(DISPLACEMENT_Y) == false)
1133  KRATOS_THROW_ERROR( std::invalid_argument,"missing DISPLACEMENT_Y dof on node ",it->Id() )
1134  if(it->HasDofFor(DISPLACEMENT_Z) == false)
1135  KRATOS_THROW_ERROR( std::invalid_argument,"missing DISPLACEMENT_Z dof on node ",it->Id() )
1136  }
1137 
1138 
1139  //check for minimum value of the buffer index
1140  //verify buffer size
1141  if (r_model_part.GetBufferSize() < 2)
1142  KRATOS_THROW_ERROR( std::logic_error, "insufficient buffer size. Buffer size should be greater than 2. Current size is", r_model_part.GetBufferSize() )
1143 
1144 
1145  return 0;
1146  KRATOS_CATCH( "" )
1147  }
1148 
1149 
1169 protected:
1170 
1171  GeneralAlphaMethod mAlpha;
1172  NewmarkMethod mNewmark;
1173 
1174  GeneralMatrices mMatrix;
1175  GeneralVectors mVector;
1176 
1177 
1182  //*********************************************************************************
1183  //Predicting Step Displacement variable
1184  //*********************************************************************************
1185 
1186  inline void PredictStepDisplacement(ModelPart::NodeType& Node,
1187  array_1d<double, 3 > & CurrentStepDisplacement,
1188  const array_1d<double, 3 > & PreviousVelocity,
1189  const array_1d<double, 3 > & PreviousAcceleration,
1190  const array_1d<double, 3 > & CurrentAcceleration,
1191  const array_1d<double, 3 > & CurrentAngularVelocity)
1192  {
1193 
1194  //CurrentStepDisplacement.clear();
1195 
1196  if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1197  {
1198  CurrentStepDisplacement[0] = 0;
1199  }
1200 
1201  if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1202  {
1203  CurrentStepDisplacement[1] = 0;
1204  }
1205 
1206  if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1207  {
1208  CurrentStepDisplacement[2] = 0;
1209  }
1210 
1211  //noalias(CurrentStepDisplacement) = ( mNewmark.c4 * PreviousVelocity + mNewmark.c5 * PreviousAcceleration + CurrentAcceleration) * mNewmark.beta * mNewmark.static_dynamic;
1212 
1213  //noalias(CurrentStepDisplacement) = ( mNewmark.deltatime * PreviousVelocity ) * mNewmark.static_dynamic;
1214 
1215  // // Angular velocity correction:
1216 
1217  // QuaternionType AngularVelocityQuaternion;
1218 
1219  // Vector AngularVelocityVector = ZeroVector(3);
1220 
1221  // for( unsigned int j=0; j<3; j++)
1222  // {
1223  // AngularVelocityVector[j] = mNewmark.c4 * mNewmark.beta * CurrentAngularVelocity[j];
1224  // }
1225 
1226  // //updated total rotations:
1227  // AngularVelocityQuaternion = QuaternionType::FromRotationVector( CurrentAngularVelocity );
1228 
1229  // Vector CurrentStepDisplacementVector = ZeroVector(3);
1230  // for( unsigned int j=0; j<3; j++)
1231  // {
1232  // CurrentStepDisplacementVector[j] = CurrentStepDisplacement[j];
1233  // }
1234 
1235  // AngularVelocityQuaternion.RotateVector3(CurrentStepDisplacementVector);
1236 
1237  // for( unsigned int j=0; j<3; j++)
1238  // {
1239  // CurrentStepDisplacement[j] = CurrentStepDisplacementVector[j];
1240  // }
1241  }
1242 
1243  //*********************************************************************************
1244  //Predicting displacement variable
1245  //*********************************************************************************
1246 
1247  inline void PredictDisplacement(array_1d<double, 3 > & CurrentDisplacement,
1248  const array_1d<double, 3 > & StepDisplacement)
1249  {
1250  noalias(CurrentDisplacement) += StepDisplacement; //needed:: to capture the correct displacements
1251  }
1252 
1253  //*********************************************************************************
1254  //Predicting first time Derivative
1255  //*********************************************************************************
1256 
1257 
1258  inline void PredictVelocity(ModelPart::NodeType& Node,
1259  array_1d<double, 3 > & CurrentVelocity,
1260  const array_1d<double, 3 > & StepDisplacement,
1261  const array_1d<double, 3 > & CurrentAcceleration,
1262  const array_1d<double, 3 > & PreviousAcceleration,
1263  const array_1d<double, 3 > & PreviousVelocity)
1264 
1265  {
1266 
1267  // BELYTSCHKO:
1268  // if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1269  // {
1270 
1271  // CurrentVelocity[0] = ( mNewmark.c0 * StepDisplacement[0]
1272  // - ( (mNewmark.gamma / mNewmark.beta) - 1.0 ) * PreviousVelocity[0]
1273  // - ( mNewmark.deltatime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) ) * PreviousAcceleration[0]) * mNewmark.static_dynamic;
1274  // }
1275 
1276  // if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1277  // {
1278  // CurrentVelocity[1] = ( mNewmark.c0 * StepDisplacement[1]
1279  // - ( (mNewmark.gamma / mNewmark.beta) - 1.0 ) * PreviousVelocity[1]
1280  // - ( mNewmark.deltatime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) ) * PreviousAcceleration[1]) * mNewmark.static_dynamic;
1281  // }
1282 
1283  // if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1284  // {
1285  // CurrentVelocity[2] = ( mNewmark.c0 * StepDisplacement[2]
1286  // - ( (mNewmark.gamma / mNewmark.beta) - 1.0 ) * PreviousVelocity[2]
1287  // - ( mNewmark.deltatime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) ) * PreviousAcceleration[2]) * mNewmark.static_dynamic;
1288  // }
1289 
1290 
1291  //SIMO:
1292  if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1293  {
1294  CurrentVelocity[0] = PreviousVelocity[0] + (mNewmark.c2 * PreviousAcceleration[0]
1295  + mNewmark.c3 * CurrentAcceleration[0] ) * mNewmark.static_dynamic;
1296  }
1297  if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1298  {
1299  CurrentVelocity[1] = PreviousVelocity[1] + (mNewmark.c2 * PreviousAcceleration[1]
1300  + mNewmark.c3 * CurrentAcceleration[1] ) * mNewmark.static_dynamic;
1301  }
1302  if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1303  {
1304  CurrentVelocity[2] = PreviousVelocity[2] + (mNewmark.c2 * PreviousAcceleration[2]
1305  + mNewmark.c3 * CurrentAcceleration[2] ) * mNewmark.static_dynamic;
1306  }
1307 
1308  }
1309  //*********************************************************************************
1310  //Predicting second time Derivative
1311  //*********************************************************************************
1312 
1313  inline void PredictAcceleration(ModelPart::NodeType& Node,
1314  array_1d<double, 3 > & CurrentAcceleration,
1315  const array_1d<double, 3 > & StepDisplacement,
1316  const array_1d<double, 3 > & PreviousAcceleration,
1317  const array_1d<double, 3 > & PreviousVelocity)
1318  {
1319  //if(Node.Is(SLAVE)) return;
1320 
1321  // BELYTSCHKO::
1322  if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1323  {
1324  CurrentAcceleration[0] = ( mNewmark.c1 * StepDisplacement[0] - mNewmark.c6 * PreviousVelocity[0] - mNewmark.c7 * PreviousAcceleration[0]) * mNewmark.static_dynamic;
1325  }
1326 
1327  if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1328  {
1329  CurrentAcceleration[1] = ( mNewmark.c1 * StepDisplacement[1] - mNewmark.c6 * PreviousVelocity[1] - mNewmark.c7 * PreviousAcceleration[1]) * mNewmark.static_dynamic;
1330  }
1331 
1332  if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1333  {
1334  CurrentAcceleration[2] = ( mNewmark.c1 * StepDisplacement[2] - mNewmark.c6 * PreviousVelocity[2] - mNewmark.c7 * PreviousAcceleration[2]) * mNewmark.static_dynamic;
1335  }
1336 
1337  // SIMO:
1338  // if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1339  // {
1340  // //CurrentAcceleration[0] = (- mNewmark.c4 * PreviousVelocity[0] + mNewmark.c5 * PreviousAcceleration[0]) * mNewmark.static_dynamic;
1341 
1342  // //CurrentAcceleration[0] = PreviousAcceleration[0]; //not stable
1343  // CurrentAcceleration[0] = 0; //ok
1344  // }
1345 
1346  // if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1347  // {
1348  // //CurrentAcceleration[1] = (- mNewmark.c4 * PreviousVelocity[1] + mNewmark.c5 * PreviousAcceleration[1]) * mNewmark.static_dynamic;
1349 
1350  // //CurrentAcceleration[1] = PreviousAcceleration[1]; //not stable
1351  // CurrentAcceleration[1] = 0; //ok
1352  // }
1353 
1354  // if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1355  // {
1356  // //CurrentAcceleration[2] = (- mNewmark.c4 * PreviousVelocity[2] + mNewmark.c5 * PreviousAcceleration[2]) * mNewmark.static_dynamic;
1357 
1358  // //CurrentAcceleration[2] = PreviousAcceleration[2]; //not stable
1359  // CurrentAcceleration[2] = 0; //ok
1360  // }
1361 
1362  }
1363 
1364  //*********************************************************************************
1365  //Predicting step rotation variable
1366  //*********************************************************************************
1367 
1368  inline void PredictStepRotation(ModelPart::NodeType& Node,
1369  array_1d<double, 3 > & CurrentStepRotation,
1370  const array_1d<double, 3 > & PreviousVelocity,
1371  const array_1d<double, 3 > & PreviousAcceleration,
1372  const array_1d<double, 3 > & CurrentAcceleration)
1373  {
1374 
1375  //CurrentStepRotation.clear(); //not needed :rigid bodies remain without rotation
1376  if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1377  {
1378  CurrentStepRotation[0] = 0;
1379  }
1380 
1381  if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1382  {
1383  CurrentStepRotation[1] = 0;
1384  }
1385 
1386  if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1387  {
1388  CurrentStepRotation[2] = 0;
1389  }
1390 
1391  //noalias(CurrentStepRotation) = ( mNewmark.c4 * PreviousVelocity + mNewmark.c5 * PreviousAcceleration + CurrentAcceleration) * mNewmark.beta * mNewmark.static_dynamic;
1392 
1393  //noalias(CurrentStepRotation) = ( mNewmark.deltatime * PreviousVelocity ) * mNewmark.static_dynamic;
1394 
1395 
1396  }
1397 
1398 
1399 
1400  //*********************************************************************************
1401  //Predicting rotation variable
1402  //*********************************************************************************
1403 
1404  inline void PredictRotation(array_1d<double, 3 > & CurrentRotation,
1405  const array_1d<double, 3 > & StepRotation)
1406  {
1407 
1408  noalias(CurrentRotation) += StepRotation; //needed for imposed rotations
1409 
1410  // QuaternionType StepRotationQuaternion;
1411  // QuaternionType CurrentRotationQuaternion;
1412  // QuaternionType ReferenceRotationQuaternion;
1413 
1414  // Vector PreviousRotationVector = ZeroVector(3);
1415  // Vector StepRotationVector = ZeroVector(3);
1416 
1417  // for( unsigned int j=0; j<3; j++)
1418  // {
1419  // PreviousRotationVector[j] = PreviousRotation[j];
1420  // StepRotationVector[j] = StepRotation[j];
1421  // }
1422 
1423  // // updated total rotations:
1424  // StepRotationQuaternion = QuaternionType::FromRotationVector( StepRotationVector );
1425  // ReferenceRotationQuaternion = QuaternionType::FromRotationVector( PreviousRotationVector );
1426 
1427  // CurrentRotationQuaternion = StepRotationQuaternion * ReferenceRotationQuaternion;
1428  // CurrentRotationQuaternion.ToRotationVector( PreviousRotationVector );
1429 
1430  // for( unsigned int j=0; j<3; j++)
1431  // {
1432  // CurrentRotation[j] = PreviousRotationVector[j];
1433  // }
1434 
1435  }
1436 
1437  //*********************************************************************************
1438  //Predicting first time Derivative
1439  //*********************************************************************************
1440 
1441  inline void PredictAngularVelocity(ModelPart::NodeType& Node,
1442  array_1d<double, 3 > & CurrentAngularVelocity,
1443  const array_1d<double, 3 > & StepRotation,
1444  const array_1d<double, 3 > & CurrentAngularAcceleration,
1445  const array_1d<double, 3 > & PreviousAngularAcceleration,
1446  const array_1d<double, 3 > & PreviousAngularVelocity)
1447  {
1448 
1449  if(Node.Is(SLAVE)) return;
1450 
1451  // BELYTSCHKO:
1452  // if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1453  // {
1454 
1455  // CurrentAngularVelocity[0] = ( mNewmark.c0 * StepRotation[0]
1456  // - ( (mNewmark.gamma / mNewmark.beta) - 1.0 ) * PreviousAngularVelocity[0]
1457  // - ( mNewmark.deltatime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) ) * PreviousAngularAcceleration[0]) * mNewmark.static_dynamic;
1458  // }
1459 
1460  // if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1461  // {
1462  // CurrentAngularVelocity[1] = ( mNewmark.c0 * StepRotation[1]
1463  // - ( (mNewmark.gamma / mNewmark.beta) - 1.0 ) * PreviousAngularVelocity[1]
1464  // - ( mNewmark.deltatime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) ) * PreviousAngularAcceleration[1]) * mNewmark.static_dynamic;
1465  // }
1466 
1467  // if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1468  // {
1469  // CurrentAngularVelocity[2] = ( mNewmark.c0 * StepRotation[2]
1470  // - ( (mNewmark.gamma / mNewmark.beta) - 1.0 ) * PreviousAngularVelocity[2]
1471  // - ( mNewmark.deltatime * 0.5 * ( ( mNewmark.gamma / mNewmark.beta ) - 2 ) ) * PreviousAngularAcceleration[2]) * mNewmark.static_dynamic;
1472  // }
1473 
1474 
1475  // SIMO:
1476  if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1477  {
1478  CurrentAngularVelocity[0] = PreviousAngularVelocity[0] + (mNewmark.c2 * PreviousAngularAcceleration[0]
1479  + mNewmark.c3 * CurrentAngularAcceleration[0] ) * mNewmark.static_dynamic;
1480 
1481  }
1482  if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1483  {
1484  CurrentAngularVelocity[1] = PreviousAngularVelocity[1] + (mNewmark.c2 * PreviousAngularAcceleration[1]
1485  + mNewmark.c3 * CurrentAngularAcceleration[1] ) * mNewmark.static_dynamic;
1486 
1487  }
1488  if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1489  {
1490  CurrentAngularVelocity[2] = PreviousAngularVelocity[2] + (mNewmark.c2 * PreviousAngularAcceleration[2]
1491  + mNewmark.c3 * CurrentAngularAcceleration[2] ) * mNewmark.static_dynamic;
1492 
1493  }
1494  }
1495 
1496  //*********************************************************************************
1497  //Predicting second time Derivative
1498  //*********************************************************************************
1499 
1500  inline void PredictAngularAcceleration(ModelPart::NodeType& Node,
1501  array_1d<double, 3 > & CurrentAngularAcceleration,
1502  const array_1d<double, 3 > & StepRotation,
1503  const array_1d<double, 3 > & PreviousAngularAcceleration,
1504  const array_1d<double, 3 > & PreviousAngularVelocity)
1505  {
1506  if(Node.Is(SLAVE)) return;
1507 
1508  // BELYTSCHKO:
1509  if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1510  {
1511  CurrentAngularAcceleration[0] = ( mNewmark.c1 * StepRotation[0] - mNewmark.c6 * PreviousAngularVelocity[0] - mNewmark.c7 * PreviousAngularAcceleration[0]) * mNewmark.static_dynamic;
1512  }
1513 
1514  if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1515  {
1516  CurrentAngularAcceleration[1] = ( mNewmark.c1 * StepRotation[1] - mNewmark.c6 * PreviousAngularVelocity[1] - mNewmark.c7 * PreviousAngularAcceleration[1]) * mNewmark.static_dynamic;
1517 
1518  }
1519 
1520  if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1521  {
1522  CurrentAngularAcceleration[2] = ( mNewmark.c1 * StepRotation[2] - mNewmark.c6 * PreviousAngularVelocity[2] - mNewmark.c7 * PreviousAngularAcceleration[2]) * mNewmark.static_dynamic;
1523  }
1524 
1525 
1526  //SIMO
1527  // if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1528  // {
1529  // //CurrentAngularAcceleration[0] = (- mNewmark.c4 * PreviousAngularVelocity[0] + mNewmark.c5 * PreviousAngularAcceleration[0]) * mNewmark.static_dynamic;
1530 
1531  // //CurrentAngularAcceleration[0] = PreviousAngularAcceleration[0]; //not stable
1532  // CurrentAngularAcceleration[0] = 0; //ok
1533  // }
1534 
1535  // if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1536  // {
1537  // //CurrentAngularAcceleration[1] = (- mNewmark.c4 * PreviousAngularVelocity[1] + mNewmark.c5 * PreviousAngularAcceleration[1]) * mNewmark.static_dynamic;
1538 
1539  // //CurrentAngularAcceleration[1] = PreviousAngularAcceleration[1]; //not stable
1540  // CurrentAngularAcceleration[1] = 0; //ok
1541  // }
1542 
1543  // if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1544  // {
1545  // //CurrentAngularAcceleration[2] = (- mNewmark.c4 * PreviousAngularVelocity[2] + mNewmark.c5 * PreviousAngularAcceleration[2]) * mNewmark.static_dynamic;
1546 
1547  // //CurrentAngularAcceleration[2] = PreviousAngularAcceleration[2]; //not stable
1548  // CurrentAngularAcceleration[2] = 0; //ok
1549  // }
1550 
1551 
1552  }
1553 
1554 
1555  //*********************************************************************************
1556  //Updating first time Derivative
1557  //*********************************************************************************
1558 
1559  inline void UpdateVelocity(ModelPart::NodeType& Node,
1560  array_1d<double, 3 > & CurrentVelocity,
1561  const array_1d<double, 3 > & DeltaDisplacement,
1562  const array_1d<double, 3 > & CurrentAcceleration,
1563  const array_1d<double, 3 > & PreviousAcceleration,
1564  const array_1d<double, 3 > & PreviousVelocity)
1565  {
1566 
1567  if(Node.Is(SLAVE)) return;
1568 
1569  // BELYTSCHKO:
1570  if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1571  {
1572  CurrentVelocity[0] = PreviousVelocity[0] + mNewmark.c2 * PreviousAcceleration[0]
1573  + mNewmark.c3 * CurrentAcceleration[0];
1574  }
1575 
1576  if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1577  {
1578  CurrentVelocity[1] = PreviousVelocity[1] + mNewmark.c2 * PreviousAcceleration[1]
1579  + mNewmark.c3 * CurrentAcceleration[1];
1580  }
1581 
1582  if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1583  {
1584  CurrentVelocity[2] = PreviousVelocity[2] + mNewmark.c2 * PreviousAcceleration[2]
1585  + mNewmark.c3 * CurrentAcceleration[2];
1586  }
1587 
1588  // SIMO:
1589  // if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1590  // {
1591  // CurrentVelocity[0] += mNewmark.c0 * DeltaDisplacement[0];
1592  // }
1593 
1594  // if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1595  // {
1596  // CurrentVelocity[1] += mNewmark.c0 * DeltaDisplacement[1];
1597  // }
1598 
1599  // if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1600  // {
1601  // CurrentVelocity[2] += mNewmark.c0 * DeltaDisplacement[2];
1602  // }
1603 
1604 
1605  }
1606 
1607  //*********************************************************************************
1608  //Updating second time Derivative
1609  //*********************************************************************************
1610 
1611 
1612  inline void UpdateAcceleration(ModelPart::NodeType& Node,
1613  array_1d<double, 3 > & CurrentAcceleration,
1614  const array_1d<double, 3 > & DeltaDisplacement,
1615  const array_1d<double, 3 > & StepDisplacement,
1616  const array_1d<double, 3 > & PreviousAcceleration,
1617  const array_1d<double, 3 > & PreviousVelocity)
1618  {
1619 
1620  if(Node.Is(SLAVE)) return;
1621 
1622  // BELYTSCHKO:
1623  if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1624  {
1625  CurrentAcceleration[0] = ( mNewmark.c1 * StepDisplacement[0]- mNewmark.c6 * PreviousVelocity[0] - mNewmark.c7 * PreviousAcceleration[0]) * mNewmark.static_dynamic;
1626  }
1627 
1628  if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1629  {
1630  CurrentAcceleration[1] = ( mNewmark.c1 * StepDisplacement[1]- mNewmark.c6 * PreviousVelocity[1] - mNewmark.c7 * PreviousAcceleration[1]) * mNewmark.static_dynamic;
1631  }
1632 
1633  if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1634  {
1635  CurrentAcceleration[2] = ( mNewmark.c1 * StepDisplacement[2]- mNewmark.c6 * PreviousVelocity[2] - mNewmark.c7 * PreviousAcceleration[2]) * mNewmark.static_dynamic;
1636  }
1637 
1638  // SIMO:
1639  // if ((Node.pGetDof(DISPLACEMENT_X))->IsFixed() == false)
1640  // {
1641  // CurrentAcceleration[0] += mNewmark.c1 * DeltaDisplacement[0];
1642  // }
1643 
1644  // if ((Node.pGetDof(DISPLACEMENT_Y))->IsFixed() == false)
1645  // {
1646  // CurrentAcceleration[1] += mNewmark.c1 * DeltaDisplacement[1];
1647  // }
1648 
1649  // if ((Node.pGetDof(DISPLACEMENT_Z))->IsFixed() == false)
1650  // {
1651  // CurrentAcceleration[2] += mNewmark.c1 * DeltaDisplacement[2];
1652  // }
1653 
1654  }
1655 
1656 
1657 
1658 
1659  //*********************************************************************************
1660  //Updating first time Derivative
1661  //*********************************************************************************
1662 
1663 
1664  inline void UpdateAngularVelocity(ModelPart::NodeType& Node,
1665  array_1d<double, 3 > & CurrentAngularVelocity,
1666  const array_1d<double, 3 > & DeltaRotation,
1667  const array_1d<double, 3 > & CurrentAngularAcceleration,
1668  const array_1d<double, 3 > & PreviousAngularAcceleration,
1669  const array_1d<double, 3 > & PreviousAngularVelocity)
1670  {
1671 
1672  if(Node.Is(SLAVE)) return;
1673 
1674  // BELYTSCHKO:
1675  if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1676  {
1677  CurrentAngularVelocity[0] = ( PreviousAngularVelocity[0] + (1-mNewmark.gamma) * mNewmark.deltatime * PreviousAngularAcceleration[0]
1678  + mNewmark.gamma * mNewmark.deltatime * CurrentAngularAcceleration[0]) * mNewmark.static_dynamic;
1679  }
1680 
1681  if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1682  {
1683  CurrentAngularVelocity[1] = ( PreviousAngularVelocity[1] + (1-mNewmark.gamma) * mNewmark.deltatime * PreviousAngularAcceleration[1]
1684  + mNewmark.gamma * mNewmark.deltatime * CurrentAngularAcceleration[1]) * mNewmark.static_dynamic;
1685  }
1686 
1687  if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1688  {
1689  CurrentAngularVelocity[2] = ( PreviousAngularVelocity[2] + (1-mNewmark.gamma) * mNewmark.deltatime * PreviousAngularAcceleration[2]
1690  + mNewmark.gamma * mNewmark.deltatime * CurrentAngularAcceleration[2]) * mNewmark.static_dynamic;
1691  }
1692 
1693  // SIMO:
1694  // if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1695  // {
1696  // CurrentAngularVelocity[0] += mNewmark.c0 * DeltaRotation[0];
1697  // }
1698 
1699  // if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1700  // {
1701  // CurrentAngularVelocity[1] += mNewmark.c0 * DeltaRotation[1];
1702  // }
1703 
1704  // if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1705  // {
1706  // CurrentAngularVelocity[2] += mNewmark.c0 * DeltaRotation[2];
1707  // }
1708  }
1709 
1710 
1711  //*********************************************************************************
1712  //Updating second time Derivative
1713  //*********************************************************************************
1714 
1715 
1716  inline void UpdateAngularAcceleration(ModelPart::NodeType& Node,
1717  array_1d<double, 3 > & CurrentAngularAcceleration,
1718  const array_1d<double, 3 > & DeltaRotation,
1719  const array_1d<double, 3 > & StepRotation,
1720  const array_1d<double, 3 > & PreviousAngularAcceleration,
1721  const array_1d<double, 3 > & PreviousAngularVelocity)
1722  {
1723  if(Node.Is(SLAVE)) return;
1724 
1725  //BELYTSCHKO:
1726  if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1727  {
1728  CurrentAngularAcceleration[0] = ( mNewmark.c1 * StepRotation[0]- mNewmark.c6 * PreviousAngularVelocity[0] - mNewmark.c7 * PreviousAngularAcceleration[0]) * mNewmark.static_dynamic;
1729  }
1730 
1731  if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1732  {
1733  CurrentAngularAcceleration[1] = ( mNewmark.c1 * StepRotation[1]- mNewmark.c6 * PreviousAngularVelocity[1] - mNewmark.c7 * PreviousAngularAcceleration[1]) * mNewmark.static_dynamic;
1734  }
1735 
1736  if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1737  {
1738  CurrentAngularAcceleration[2] = ( mNewmark.c1 * StepRotation[2]- mNewmark.c6 * PreviousAngularVelocity[2] - mNewmark.c7 * PreviousAngularAcceleration[2]) * mNewmark.static_dynamic;
1739  }
1740 
1741  // SIMO:
1742  // if ((Node.pGetDof(ROTATION_X))->IsFixed() == false)
1743  // {
1744  // CurrentAngularAcceleration[0] += mNewmark.c1 * DeltaRotation[0];
1745  // }
1746 
1747  // if ((Node.pGetDof(ROTATION_Y))->IsFixed() == false)
1748  // {
1749  // CurrentAngularAcceleration[1] += mNewmark.c1 * DeltaRotation[1];
1750  // }
1751 
1752  // if ((Node.pGetDof(ROTATION_Z))->IsFixed() == false)
1753  // {
1754  // CurrentAngularAcceleration[2] += mNewmark.c1 * DeltaRotation[2];
1755  // }
1756  }
1757 
1758 
1759  //Elements and Conditions:
1760  //****************************************************************************
1761 
1766  void AddDynamicsToLHS(
1767  LocalSystemMatrixType& LHS_Contribution,
1768  LocalSystemMatrixType& D,
1769  LocalSystemMatrixType& M,
1770  ProcessInfo& CurrentProcessInfo)
1771  {
1772 
1773  // adding mass contribution to the dynamic stiffness
1774  if (M.size1() != 0) // if M matrix declared
1775  {
1776  noalias(LHS_Contribution) += M * mNewmark.static_dynamic;
1777 
1778  //std::cout<<" Mass Matrix "<<M<<" coeficient "<<mNewmark.c0<<std::endl;
1779  }
1780 
1781  //adding damping contribution
1782  if (D.size1() != 0) // if M matrix declared
1783  {
1784  noalias(LHS_Contribution) += D * mNewmark.static_dynamic;
1785 
1786  }
1787 
1788  }
1789 
1794  void AddDynamicsToRHS(
1795  LocalSystemVectorType& RHS_Contribution,
1796  LocalSystemVectorType& fd,
1797  LocalSystemVectorType& fm,
1798  ProcessInfo& CurrentProcessInfo)
1799  {
1800 
1801  //adding inertia contribution
1802  if (fm.size() != 0)
1803  {
1804  noalias(RHS_Contribution) -= mNewmark.static_dynamic * fm;
1805  }
1806 
1807  //adding damping contribution
1808  if (fd.size() != 0)
1809  {
1810  noalias(RHS_Contribution) -= mNewmark.static_dynamic * fd;
1811  }
1812 
1813 
1814  }
1815 
1816 
1817  //Elements:
1818  //****************************************************************************
1819 
1824  void AddDynamicsToRHS(
1825  Element::Pointer rCurrentElement,
1826  LocalSystemVectorType& RHS_Contribution,
1827  LocalSystemMatrixType& D,
1828  LocalSystemMatrixType& M,
1829  ProcessInfo& CurrentProcessInfo)
1830  {
1831  int thread = OpenMPUtils::ThisThread();
1832 
1833  //adding inertia contribution
1834  if (M.size1() != 0)
1835  {
1836  rCurrentElement->GetSecondDerivativesVector(mVector.a[thread], 0);
1837 
1838  (mVector.a[thread]) *= mNewmark.static_dynamic ;
1839 
1840  rCurrentElement->GetSecondDerivativesVector(mVector.ap[thread], 1);
1841 
1842  noalias(mVector.a[thread]) += mVector.ap[thread] * mNewmark.static_dynamic;
1843 
1844  noalias(RHS_Contribution) -= prod(M, mVector.a[thread]);
1845  //KRATOS_WATCH( prod(M, macc[thread] ) )
1846 
1847  }
1848 
1849  //adding damping contribution
1850  if (D.size1() != 0)
1851  {
1852  rCurrentElement->GetFirstDerivativesVector(mVector.v[thread], 0);
1853 
1854  (mVector.v[thread]) *= mNewmark.static_dynamic ;
1855 
1856  noalias(RHS_Contribution) -= prod(D, mVector.v[thread]);
1857  }
1858 
1859 
1860  }
1861 
1862 
1863 
1864  //Conditions:
1865  //****************************************************************************
1866 
1867 
1868  void AddDynamicsToRHS(
1869  Condition::Pointer rCurrentCondition,
1870  LocalSystemVectorType& RHS_Contribution,
1871  LocalSystemMatrixType& D,
1872  LocalSystemMatrixType& M,
1873  ProcessInfo& CurrentProcessInfo)
1874  {
1875  int thread = OpenMPUtils::ThisThread();
1876 
1877  //adding inertia contribution
1878  if (M.size1() != 0)
1879  {
1880  rCurrentCondition->GetSecondDerivativesVector(mVector.a[thread], 0);
1881 
1882  (mVector.a[thread]) *= mNewmark.static_dynamic;
1883 
1884  rCurrentCondition->GetSecondDerivativesVector(mVector.ap[thread], 1);
1885 
1886  noalias(mVector.a[thread]) += mVector.ap[thread] * mNewmark.static_dynamic;
1887 
1888  noalias(RHS_Contribution) -= prod(M, mVector.a[thread]);
1889  }
1890 
1891  //adding damping contribution
1892  //damping contribution
1893  if (D.size1() != 0)
1894  {
1895  rCurrentCondition->GetFirstDerivativesVector(mVector.v[thread], 0);
1896 
1897  (mVector.v[thread]) *= mNewmark.static_dynamic ;
1898 
1899  noalias(RHS_Contribution) -= prod(D, mVector.v [thread]);
1900  }
1901 
1902  }
1903 
1904 
1905 
1924 private:
1930  //DofsVectorType mElementalDofList;
1933 
1946 }; /* Class Scheme */
1947 } /* namespace Kratos.*/
1948 
1949 #endif /* KRATOS_RESIDUAL_BASED_BOSSAK_DISPLACEMENT_ROTATION_SCHEME defined */
array_1d.h
beam_math_utilities.hpp
Kratos::BeamMathUtils
Definition: beam_math_utilities.hpp:31
Kratos::BeamMathUtils::VectorToSkewSymmetricTensor
static void VectorToSkewSymmetricTensor(const TVector3 &rVector, TMatrix3 &rSkewSymmetricTensor)
Definition: beam_math_utilities.hpp:231
Kratos::Element
Base class for all Elements.
Definition: element.h:60
Kratos::Element::DofsVectorType
std::vector< DofType::Pointer > DofsVectorType
Definition: element.h:100
Kratos::Element::EquationIdVectorType
std::vector< std::size_t > EquationIdVectorType
Definition: element.h:98
Kratos::Flags::Is
bool Is(Flags const &rOther) const
Definition: flags.h:274
Kratos::GeometricalObject::GetGeometry
GeometryType & GetGeometry()
Returns the reference of the geometry.
Definition: geometrical_object.h:158
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::ConditionsContainerType
MeshType::ConditionsContainerType ConditionsContainerType
Condintions container. A vector set of Conditions with their Id's as key.
Definition: model_part.h:183
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::GetBufferSize
IndexType GetBufferSize() const
This method gets the suffer size of the model part database.
Definition: model_part.h:1876
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::NodeIterator
MeshType::NodeIterator NodeIterator
Definition: model_part.h:134
Kratos::ModelPart::Elements
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
Kratos::ModelPart::NodesEnd
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::Node::FastGetSolutionStepValue
TVariableType::Type & FastGetSolutionStepValue(const TVariableType &rThisVariable)
Definition: node.h:435
Kratos::Node::pGetDof
DofType::Pointer pGetDof(TVariableType const &rDofVariable) const
Get DoF counted pointer for a given variable.
Definition: node.h:711
Kratos::Node::GetInitialPosition
const PointType & GetInitialPosition() const
Definition: node.h:559
Kratos::OpenMPUtils::DivideInPartitions
static void DivideInPartitions(const int NumTerms, const int NumThreads, PartitionVector &Partitions)
Divide an array of length NumTerms between NumThreads threads.
Definition: openmp_utils.h:158
Kratos::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::OpenMPUtils::PartitionVector
std::vector< int > PartitionVector
Vector type for the output of DivideInPartitions method.
Definition: openmp_utils.h:53
Kratos::PointerVectorSet
A sorted associative container similar to an STL set, but uses a vector to store pointers to its data...
Definition: pointer_vector_set.h:72
Kratos::PointerVectorSet::iterator
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: pointer_vector_set.h:95
Kratos::PointerVectorSet::begin
iterator begin()
Returns an iterator pointing to the beginning of the container.
Definition: pointer_vector_set.h:278
Kratos::PointerVectorSet::end
iterator end()
Returns an iterator pointing to the end of the container.
Definition: pointer_vector_set.h:314
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::Quaternion::ToRotationMatrix
void ToRotationMatrix(TMatrix3x3 &R) const
Definition: quaternion.h:213
Kratos::Quaternion::ToRotationVector
void ToRotationVector(T &rx, T &ry, T &rz) const
Definition: quaternion.h:262
Kratos::Quaternion< double >::FromRotationVector
static Quaternion FromRotationVector(double rx, double ry, double rz)
Definition: quaternion.h:427
Kratos::ResidualBasedBossakDisplacementRotationScheme
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:51
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictAngularVelocity
void PredictAngularVelocity(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentAngularVelocity, const array_1d< double, 3 > &StepRotation, const array_1d< double, 3 > &CurrentAngularAcceleration, const array_1d< double, 3 > &PreviousAngularAcceleration, const array_1d< double, 3 > &PreviousAngularVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1441
Kratos::ResidualBasedBossakDisplacementRotationScheme::mVector
GeneralVectors mVector
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1175
Kratos::ResidualBasedBossakDisplacementRotationScheme::BaseType
Scheme< TSparseSpace, TDenseSpace > BaseType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:108
Kratos::ResidualBasedBossakDisplacementRotationScheme::CalculateSystemContributions
void CalculateSystemContributions(Element::Pointer rCurrentElement, LocalSystemMatrixType &LHS_Contribution, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:884
Kratos::ResidualBasedBossakDisplacementRotationScheme::InitializeElements
void InitializeElements(ModelPart &rModelPart)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:655
Kratos::ResidualBasedBossakDisplacementRotationScheme::Condition_CalculateSystemContributions
void Condition_CalculateSystemContributions(Condition::Pointer rCurrentCondition, LocalSystemMatrixType &LHS_Contribution, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:972
Kratos::ResidualBasedBossakDisplacementRotationScheme::InitializeNonLinearIteration
void InitializeNonLinearIteration(Element::Pointer rCurrentElement, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:865
Kratos::ResidualBasedBossakDisplacementRotationScheme::QuaternionType
Quaternion< double > QuaternionType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:130
Kratos::ResidualBasedBossakDisplacementRotationScheme::AddDynamicsToRHS
void AddDynamicsToRHS(LocalSystemVectorType &RHS_Contribution, LocalSystemVectorType &fd, LocalSystemVectorType &fm, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1794
Kratos::ResidualBasedBossakDisplacementRotationScheme::ConditionsArrayType
ModelPart::ConditionsContainerType ConditionsArrayType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:126
Kratos::ResidualBasedBossakDisplacementRotationScheme::GetElementalDofList
void GetElementalDofList(Element::Pointer rCurrentElement, Element::DofsVectorType &ElementalDofList, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1066
Kratos::ResidualBasedBossakDisplacementRotationScheme::SlaveNodesUpdate
void SlaveNodesUpdate(ModelPart &r_model_part)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:403
Kratos::ResidualBasedBossakDisplacementRotationScheme::Condition_Calculate_RHS_Contribution
void Condition_Calculate_RHS_Contribution(Condition::Pointer rCurrentCondition, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1020
Kratos::ResidualBasedBossakDisplacementRotationScheme::DofsArrayType
BaseType::DofsArrayType DofsArrayType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:112
Kratos::ResidualBasedBossakDisplacementRotationScheme::UpdateAcceleration
void UpdateAcceleration(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentAcceleration, const array_1d< double, 3 > &DeltaDisplacement, const array_1d< double, 3 > &StepDisplacement, const array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &PreviousVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1612
Kratos::ResidualBasedBossakDisplacementRotationScheme::BeamMathUtilsType
BeamMathUtils< double > BeamMathUtilsType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:128
Kratos::ResidualBasedBossakDisplacementRotationScheme::DofsVectorType
Element::DofsVectorType DofsVectorType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:114
Kratos::ResidualBasedBossakDisplacementRotationScheme::Predict
void Predict(ModelPart &r_model_part, DofsArrayType &rDofSet, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b)
Performing the prediction of the solution.
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:534
Kratos::ResidualBasedBossakDisplacementRotationScheme::~ResidualBasedBossakDisplacementRotationScheme
virtual ~ResidualBasedBossakDisplacementRotationScheme()
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:174
Kratos::ResidualBasedBossakDisplacementRotationScheme::InitializeConditions
void InitializeConditions(ModelPart &rModelPart)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:687
Kratos::ResidualBasedBossakDisplacementRotationScheme::TSystemMatrixType
BaseType::TSystemMatrixType TSystemMatrixType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:116
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictAcceleration
void PredictAcceleration(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentAcceleration, const array_1d< double, 3 > &StepDisplacement, const array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &PreviousVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1313
Kratos::ResidualBasedBossakDisplacementRotationScheme::TDataType
BaseType::TDataType TDataType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:110
Kratos::ResidualBasedBossakDisplacementRotationScheme::AddDynamicsToLHS
void AddDynamicsToLHS(LocalSystemMatrixType &LHS_Contribution, LocalSystemMatrixType &D, LocalSystemMatrixType &M, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1766
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictVelocity
void PredictVelocity(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentVelocity, const array_1d< double, 3 > &StepDisplacement, const array_1d< double, 3 > &CurrentAcceleration, const array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &PreviousVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1258
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictStepDisplacement
void PredictStepDisplacement(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentStepDisplacement, const array_1d< double, 3 > &PreviousVelocity, const array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &CurrentAcceleration, const array_1d< double, 3 > &CurrentAngularVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1186
Kratos::ResidualBasedBossakDisplacementRotationScheme::InitializeNonLinIteration
void InitializeNonLinIteration(ModelPart &r_model_part, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b)
unction to be called when it is needed to initialize an iteration. It is designed to be called at the...
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:830
Kratos::ResidualBasedBossakDisplacementRotationScheme::ElementsArrayType
ModelPart::ElementsContainerType ElementsArrayType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:124
Kratos::ResidualBasedBossakDisplacementRotationScheme::InitializeNonLinearIteration
void InitializeNonLinearIteration(Condition::Pointer rCurrentCondition, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:855
Kratos::ResidualBasedBossakDisplacementRotationScheme::GetConditionDofList
void GetConditionDofList(Condition::Pointer rCurrentCondition, Element::DofsVectorType &ConditionDofList, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1080
Kratos::ResidualBasedBossakDisplacementRotationScheme::ResidualBasedBossakDisplacementRotationScheme
ResidualBasedBossakDisplacementRotationScheme(double rDynamic=1, double mAlphaM=0.0)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:139
Kratos::ResidualBasedBossakDisplacementRotationScheme::LocalSystemMatrixType
BaseType::LocalSystemMatrixType LocalSystemMatrixType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:122
Kratos::ResidualBasedBossakDisplacementRotationScheme::Check
virtual int Check(ModelPart &r_model_part)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1098
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictAngularAcceleration
void PredictAngularAcceleration(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentAngularAcceleration, const array_1d< double, 3 > &StepRotation, const array_1d< double, 3 > &PreviousAngularAcceleration, const array_1d< double, 3 > &PreviousAngularVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1500
Kratos::ResidualBasedBossakDisplacementRotationScheme::mMatrix
GeneralMatrices mMatrix
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1174
Kratos::ResidualBasedBossakDisplacementRotationScheme::FinalizeSolutionStep
void FinalizeSolutionStep(ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:779
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictStepRotation
void PredictStepRotation(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentStepRotation, const array_1d< double, 3 > &PreviousVelocity, const array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &CurrentAcceleration)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1368
Kratos::ResidualBasedBossakDisplacementRotationScheme::UpdateAngularAcceleration
void UpdateAngularAcceleration(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentAngularAcceleration, const array_1d< double, 3 > &DeltaRotation, const array_1d< double, 3 > &StepRotation, const array_1d< double, 3 > &PreviousAngularAcceleration, const array_1d< double, 3 > &PreviousAngularVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1716
Kratos::ResidualBasedBossakDisplacementRotationScheme::InitializeSolutionStep
void InitializeSolutionStep(ModelPart &r_model_part, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:724
Kratos::ResidualBasedBossakDisplacementRotationScheme::Update
void Update(ModelPart &r_model_part, DofsArrayType &rDofSet, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:197
Kratos::ResidualBasedBossakDisplacementRotationScheme::AddDynamicsToRHS
void AddDynamicsToRHS(Element::Pointer rCurrentElement, LocalSystemVectorType &RHS_Contribution, LocalSystemMatrixType &D, LocalSystemMatrixType &M, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1824
Kratos::ResidualBasedBossakDisplacementRotationScheme::KRATOS_CLASS_POINTER_DEFINITION
KRATOS_CLASS_POINTER_DEFINITION(ResidualBasedBossakDisplacementRotationScheme)
Kratos::ResidualBasedBossakDisplacementRotationScheme::mNewmark
NewmarkMethod mNewmark
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1172
Kratos::ResidualBasedBossakDisplacementRotationScheme::mAlpha
GeneralAlphaMethod mAlpha
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1171
Kratos::ResidualBasedBossakDisplacementRotationScheme::LocalSystemVectorType
BaseType::LocalSystemVectorType LocalSystemVectorType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:120
Kratos::ResidualBasedBossakDisplacementRotationScheme::UpdateAngularVelocity
void UpdateAngularVelocity(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentAngularVelocity, const array_1d< double, 3 > &DeltaRotation, const array_1d< double, 3 > &CurrentAngularAcceleration, const array_1d< double, 3 > &PreviousAngularAcceleration, const array_1d< double, 3 > &PreviousAngularVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1664
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictRotation
void PredictRotation(array_1d< double, 3 > &CurrentRotation, const array_1d< double, 3 > &StepRotation)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1404
Kratos::ResidualBasedBossakDisplacementRotationScheme::UpdateVelocity
void UpdateVelocity(ModelPart::NodeType &Node, array_1d< double, 3 > &CurrentVelocity, const array_1d< double, 3 > &DeltaDisplacement, const array_1d< double, 3 > &CurrentAcceleration, const array_1d< double, 3 > &PreviousAcceleration, const array_1d< double, 3 > &PreviousVelocity)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1559
Kratos::ResidualBasedBossakDisplacementRotationScheme::AddDynamicsToRHS
void AddDynamicsToRHS(Condition::Pointer rCurrentCondition, LocalSystemVectorType &RHS_Contribution, LocalSystemMatrixType &D, LocalSystemMatrixType &M, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1868
Kratos::ResidualBasedBossakDisplacementRotationScheme::Calculate_RHS_Contribution
void Calculate_RHS_Contribution(Element::Pointer rCurrentElement, LocalSystemVectorType &RHS_Contribution, Element::EquationIdVectorType &EquationId, ProcessInfo &CurrentProcessInfo)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:927
Kratos::ResidualBasedBossakDisplacementRotationScheme::PredictDisplacement
void PredictDisplacement(array_1d< double, 3 > &CurrentDisplacement, const array_1d< double, 3 > &StepDisplacement)
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:1247
Kratos::ResidualBasedBossakDisplacementRotationScheme::TSystemVectorType
BaseType::TSystemVectorType TSystemVectorType
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:118
Kratos::Scheme
This class provides the implementation of the basic tasks that are needed by the solution strategy.
Definition: scheme.h:56
Kratos::Scheme::TSystemMatrixType
typename TSparseSpace::MatrixType TSystemMatrixType
Matrix type definition.
Definition: scheme.h:71
Kratos::Scheme::mConditionsAreInitialized
bool mConditionsAreInitialized
Flag taking in account if the elements were initialized correctly or not.
Definition: scheme.h:757
Kratos::Scheme::ElementsArrayType
ModelPart::ElementsContainerType ElementsArrayType
Elements containers definition.
Definition: scheme.h:89
Kratos::Scheme::EquationId
virtual void EquationId(const Element &rElement, Element::EquationIdVectorType &rEquationId, const ProcessInfo &rCurrentProcessInfo)
This method gets the eqaution id corresponding to the current element.
Definition: scheme.h:636
Kratos::Scheme::TSystemVectorType
typename TSparseSpace::VectorType TSystemVectorType
Vector type definition.
Definition: scheme.h:74
Kratos::Scheme::LocalSystemVectorType
typename TDenseSpace::VectorType LocalSystemVectorType
Local system vector type definition.
Definition: scheme.h:80
Kratos::Scheme::InitializeSolutionStep
virtual void InitializeSolutionStep(ModelPart &rModelPart, TSystemMatrixType &A, TSystemVectorType &Dx, TSystemVectorType &b)
Function called once at the beginning of each solution step.
Definition: scheme.h:272
Kratos::Scheme::TDataType
typename TSparseSpace::DataType TDataType
Data type definition.
Definition: scheme.h:68
Kratos::Scheme::LocalSystemMatrixType
typename TDenseSpace::MatrixType LocalSystemMatrixType
Local system matrix type definition.
Definition: scheme.h:77
Kratos::Scheme::Check
virtual int Check(const ModelPart &rModelPart) const
This function is designed to be called once to perform all the checks needed on the input provided....
Definition: scheme.h:508
Kratos::Scheme::mElementsAreInitialized
bool mElementsAreInitialized
Flag to be used in controlling if the Scheme has been initialized or not.
Definition: scheme.h:756
Kratos::Scheme::ConditionsArrayType
ModelPart::ConditionsContainerType ConditionsArrayType
Conditions containers definition.
Definition: scheme.h:92
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 >
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_TRY
#define KRATOS_TRY
Definition: define.h:109
model_part.h
Kratos::Python::GetValue
Parameters GetValue(Parameters &rParameters, const std::string &rEntry)
Definition: add_kratos_parameters_to_python.cpp:53
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::prod
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
generate_total_lagrangian_mixed_volumetric_strain_element.b
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
sensitivityMatrix.A
A
Definition: sensitivityMatrix.py:70
tensors::i
integer i
Definition: TensorModule.f:17
solid_mechanics_application_variables.h
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralAlphaMethod
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:56
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralAlphaMethod::m
double m
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:58
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralAlphaMethod::f
double f
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:57
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralMatrices
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:86
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralMatrices::D
std::vector< Matrix > D
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:88
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralMatrices::M
std::vector< Matrix > M
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:87
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralVectors
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:93
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralVectors::fd
std::vector< Vector > fd
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:95
Kratos::ResidualBasedBossakDisplacementRotationScheme::GeneralVectors::fm
std::vector< Vector > fm
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:94
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:63
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c5
double c5
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:75
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::deltatime
double deltatime
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:67
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c1
double c1
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:71
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::beta
double beta
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:65
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c0
double c0
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:70
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::static_dynamic
double static_dynamic
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:80
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c3
double c3
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:73
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::gamma
double gamma
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:66
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c2
double c2
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:72
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c6
double c6
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:76
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c4
double c4
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:74
Kratos::ResidualBasedBossakDisplacementRotationScheme::NewmarkMethod::c7
double c7
Definition: residual_based_bossak_displacement_rotation_scheme.hpp:77
Generated on Sat Feb 3 2024 01:58:18 for KratosMultiphysics by doxygen 1.9.1