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.
monolithic_dem_coupled.h
Go to the documentation of this file.
1 /*
2 ==============================================================================
3 Kratos
4 A General Purpose Software for Multi-Physics Finite Element Analysis
5 Version 1.0 (Released on march 05, 2007).
6 
7 Copyright 2007
8 Pooyan Dadvand, Riccardo Rossi
9 pooyan@cimne.upc.edu
10 rrossi@cimne.upc.edu
11 CIMNE (International Center for Numerical Methods in Engineering),
12 Gran Capita' s/n, 08034 Barcelona, Spain
13 
14 Permission is hereby granted, free of charge, to any person obtaining
15 a copy of this software and associated documentation files (the
16 "Software"), to deal in the Software without restriction, including
17 without limitation the rights to use, copy, modify, merge, publish,
18 distribute, sublicense and/or sell copies of the Software, and to
19 permit persons to whom the Software is furnished to do so, subject to
20 the following condition:
21 
22 Distribution of this code for any commercial purpose is permissible
23 ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNER.
24 
25 The above copyright notice and this permission notice shall be
26 included in all copies or substantial portions of the Software.
27 
28 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
29 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
30 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
31 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
32 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
33 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
34 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
35 
36 ==============================================================================
37  */
38 
39 //
40 // Project Name: Kratos
41 // Last Modified by: $Author: gcasas $
42 // Date: $Date: 2010-10-09 10:34:00 $
43 // Revision: $Revision: 0.1 $
44 //
45 //
46 
47 
48 #ifndef MONOLITHIC_DEM_COUPLED_H
49 #define MONOLITHIC_DEM_COUPLED_H
50 
51 // System includes
52 #include <string>
53 #include <iostream>
54 
55 
56 // External includes
57 
58 
59 // Project includes
60 #include "containers/array_1d.h"
61 #include "includes/define.h"
62 #include "includes/element.h"
63 #include "includes/serializer.h"
64 #include "utilities/geometry_utilities.h"
65 #include "includes/cfd_variables.h"
66 
67 // Application includes
69 
70 namespace Kratos
71 {
72 
75 
78 
82 
86 
90 
94 
96 
134 template< unsigned int TDim,
135  unsigned int TNumNodes = TDim + 1 >
136 class KRATOS_API(SWIMMING_DEM_APPLICATION) MonolithicDEMCoupled : public Element
137 {
138 public:
141 
144 
147 
149  typedef Node NodeType;
150 
156 
159 
162 
164 
166 
167  typedef std::size_t IndexType;
168 
169  typedef std::size_t SizeType;
170 
171  typedef std::vector<std::size_t> EquationIdVectorType;
172 
173  typedef std::vector< Dof<double>::Pointer > DofsVectorType;
174 
176 
177 //G
180 
183 
186 //Z
190 
191  //Constructors.
192 
194 
198  Element(NewId)
199  {}
200 
202 
206  MonolithicDEMCoupled(IndexType NewId, const NodesArrayType& ThisNodes) :
207  Element(NewId, ThisNodes)
208  {}
209 
211 
215  MonolithicDEMCoupled(IndexType NewId, GeometryType::Pointer pGeometry) :
216  Element(NewId, pGeometry)
217  {}
218 
220 
225  MonolithicDEMCoupled(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) :
226  Element(NewId, pGeometry, pProperties)
227  {}
228 
231  {}
232 
233 
237 
238 
242 
244 
251  Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes,
252  PropertiesType::Pointer pProperties) const override
253  {
254 
255  return Element::Pointer(new MonolithicDEMCoupled(NewId, GetGeometry().Create(ThisNodes), pProperties));
256 
257  }
258 
260 
269  virtual void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
270  VectorType& rRightHandSideVector,
271  const ProcessInfo& rCurrentProcessInfo) override
272  {
273  if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1) {
274  const unsigned int LocalSize = (TDim + 1) * TNumNodes;
275 
276  // Check sizes and initialize matrix
277  if (rLeftHandSideMatrix.size1() != LocalSize)
278  rLeftHandSideMatrix.resize(LocalSize, LocalSize, false);
279 
280  noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
281 
282  // Calculate RHS
283  this->CalculateRightHandSide(rRightHandSideVector, rCurrentProcessInfo);
284  }
285 
286  else {
287  const unsigned int LocalSize = TDim * TNumNodes;
288 
289  if (rLeftHandSideMatrix.size1() != LocalSize)
290  rLeftHandSideMatrix.resize(LocalSize, LocalSize, false);
291 
292  noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
293  CalculateLaplacianMassMatrix(rLeftHandSideMatrix, rCurrentProcessInfo);
294  noalias(rRightHandSideVector) = ZeroVector(LocalSize);
295  this->CalculateRightHandSide(rRightHandSideVector, rCurrentProcessInfo);
296  }
297  }
298 
300 
304  virtual void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rCurrentProcessInfo) override
305  {
306  if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1) {
307 
308  const unsigned int LocalSize = (TDim + 1) * TNumNodes;
309 
310  if (rLeftHandSideMatrix.size1() != LocalSize) rLeftHandSideMatrix.resize(LocalSize, LocalSize, false);
311 
312  noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
313  }
314 
315  else {
316  const unsigned int LocalSize = TDim * TNumNodes;
317 
318  if (rLeftHandSideMatrix.size1() != LocalSize) rLeftHandSideMatrix.resize(LocalSize, LocalSize, false);
319 
320  noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
321  CalculateLaplacianMassMatrix(rLeftHandSideMatrix, rCurrentProcessInfo);
322  }
323  }
324 
326 
335  virtual void CalculateRightHandSide(VectorType& rRightHandSideVector,
336  const ProcessInfo& rCurrentProcessInfo) override
337  {
338  // Calculate this element's geometric parameters
339  double Area;
342  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
343 
344  // Calculate this element's fluid properties
345  double Density;
346  this->EvaluateInPoint(Density, DENSITY, N);
347 
348  if (rCurrentProcessInfo[FRACTIONAL_STEP] == 1) {
349 
350  const unsigned int LocalSize = (TDim + 1) * TNumNodes;
351 
352  // Check sizes and initialize
353  if (rRightHandSideVector.size() != LocalSize)
354  rRightHandSideVector.resize(LocalSize, false);
355 
356  noalias(rRightHandSideVector) = ZeroVector(LocalSize);
357 
358  // Calculate Momentum RHS contribution
359  this->AddMomentumRHS(rRightHandSideVector, Density, N, Area);
360  //G
361  const double& DeltaTime = rCurrentProcessInfo[DELTA_TIME];
362  static const double arr[] = {1.0,-1.0};
363  std::vector<double> SchemeWeights (arr, arr + sizeof(arr) / sizeof(arr[0]));
364  this->AddMassRHS(rRightHandSideVector, Density, N, Area, SchemeWeights, DeltaTime);
365  }
366 
367  else {
368  const unsigned int LocalSize = TDim * TNumNodes;
369 
370  // Check sizes and initialize
371  if (rRightHandSideVector.size() != LocalSize)
372  rRightHandSideVector.resize(LocalSize, false);
373 
374  noalias(rRightHandSideVector) = ZeroVector(LocalSize);
375 
376  this->AddRHSLaplacian(rRightHandSideVector, DN_DX, Area);
377  }
378  // Shape functions and integration points
379 
380 /*
381  MatrixType NContainer;
382  ShapeFunctionDerivativesArrayType DN_DXContainer;
383  VectorType GaussWeights;
384  this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
385  const SizeType NumGauss = NContainer.size1();
386 
387  for (SizeType g = 0; g < NumGauss; g++){
388  const double GaussWeight = GaussWeights[g];
389  const ShapeFunctionsType& Ng = row(NContainer, g);
390  this->AddMomentumRHS(rRightHandSideVector, Density, Ng, GaussWeight);
391  this->AddMassRHS(rRightHandSideVector, Density, Ng, GaussWeight, SchemeWeights, DeltaTime);
392  }
393 
394 */
395 
396 //Z
397 
398  // For OSS: Add projection of residuals to RHS
399  if (rCurrentProcessInfo[OSS_SWITCH] == 1)
400  {
401  array_1d<double, 3 > AdvVel;
402  this->GetAdvectiveVel(AdvVel, N);
403 
404  double KinViscosity;
405  this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
406 
407  double Viscosity;
408  this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
409 
410  // Calculate stabilization parameters
411  double TauOne, TauTwo;
412  this->CalculateTau(TauOne, TauTwo, AdvVel, Area,Density, Viscosity, rCurrentProcessInfo);
413 
414  this->AddProjectionToRHS(rRightHandSideVector, AdvVel, Density, TauOne, TauTwo, N, DN_DX, Area,rCurrentProcessInfo[DELTA_TIME]);
415  }
416  /*
417 // else if (this->GetValue(TRACK_SUBSCALES) == 1)// Experimental: Dynamic tracking of ASGS subscales (see Codina 2002 Stabilized finite element ... using orthogonal subscales)
418 // {
419 // * We want to evaluate v * d(u_subscale)/dt. This term is zero in OSS, due the orthogonality of the two terms. in ASGS, we approximate it as
420 // * d(u_s)/dt = u_s(last iteration) - u_s(previous time step) / DeltaTime where u_s(last iteration) = TauOne(MomResidual(last_iteration) + u_s(previous step)/DeltaTime)
421 // *
422 // array_1d<double, 3 > AdvVel;
423 // this->GetAdvectiveVel(AdvVel, N);
424 //
425 // double KinViscosity;
426 // this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
427 //
428 // double Viscosity;
429 // this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
430 //
431 // double StaticTauOne;
432 // this->CalculateStaticTau(StaticTauOne,AdvVel,Area,Viscosity);
433 //
434 // double TauOne,TauTwo;
435 // this->CalculateTau(TauOne,TauTwo,AdvVel,Area,Density,Viscosity,rCurrentProcessInfo);
436 //
437 // array_1d<double,3> ElemMomRes(3,0.0);
438 // this->ASGSMomResidual(AdvVel,Density,ElemMomRes,N,DN_DX,1.0);
439 //
440 // const array_1d<double,3>& rOldSubscale = this->GetValue(SUBSCALE_VELOCITY);
441 // const double DeltaTime = rCurrentProcessInfo[DELTA_TIME];
442 //
443 // // array_1d<double,3> Subscale = TauOne * (ElemMomRes + Density * rOldSubscale / DeltaTime);
444 // const double C1 = 1.0 - TauOne/StaticTauOne;
445 // const double C2 = TauOne / (StaticTauOne*DeltaTime);
446 //
447 // unsigned int LocalIndex = 0;
448 // for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
449 // {
450 // for(unsigned int d = 0; d < TDim; d++)
451 // rRightHandSideVector[LocalIndex++] -= N[iNode] * Area * ( C1 * ElemMomRes[d] - C2 * rOldSubscale[d] );
452 // // rRightHandSideVector[LocalIndex++] -= N[iNode] * Area * ( Subscale[d] - rOldSubscale[d] ) / DeltaTime;
453 // LocalIndex++; // Pressure Dof
454 // }
455 // }
456 */
457  }
458 
460 
467  virtual void CalculateMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo) override
468  {
469 // rMassMatrix.resize(0,0,false);
470  const unsigned int LocalSize = (TDim + 1) * TNumNodes;
471 
472  // Resize and set to zero
473  if (rMassMatrix.size1() != LocalSize)
474  rMassMatrix.resize(LocalSize, LocalSize, false);
475 
476  rMassMatrix = ZeroMatrix(LocalSize, LocalSize);
477 
478  // Get the element's geometric parameters
479  double Area;
482  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
483 
484  // Calculate this element's fluid properties
485  double Density;
486  this->EvaluateInPoint(Density, DENSITY, N);
487 
488  // Add 'classical' mass matrix (lumped)
489  double Coeff = Density * Area / TNumNodes; //Optimize!
490 //G
491  this->CalculateLumpedMassMatrix(rMassMatrix, Coeff);
492 /*
493  MatrixType NContainer;
494  ShapeFunctionDerivativesArrayType DN_DXContainer;
495  VectorType GaussWeights;
496  this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
497  const SizeType NumGauss = NContainer.size1();
498 
499  for (SizeType g = 0; g < NumGauss; g++){
500  const double GaussWeight = GaussWeights[g];
501  const ShapeFunctionsType& Ng = row(NContainer, g);
502  this->AddConsistentMassMatrixContribution(rMassMatrix, Ng, Density, GaussWeight);
503  }
504 */
505 //Z
506  /* For ASGS: add dynamic stabilization terms.
507  These terms are not used in OSS, as they belong to the finite element
508  space and cancel out with their projections.
509  */
510  if (rCurrentProcessInfo[OSS_SWITCH] != 1)
511  {
512  double KinViscosity;
513  this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
514 
515  double Viscosity;
516  this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
517 
518  // Get Advective velocity
519  array_1d<double, 3 > AdvVel;
520  this->GetAdvectiveVel(AdvVel, N);
521 
522  // Calculate stabilization parameters
523  double TauOne, TauTwo;
524  this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
525 
526  // Add dynamic stabilization terms ( all terms involving a delta(u) )
527 /*
528  for (SizeType g = 0; g < NumGauss; g++){
529  const double GaussWeight = GaussWeights[g];
530  const ShapeFunctionsType& Ng = row(NContainer, g);
531  this->GetAdvectiveVel(AdvVel, Ng);
532  this->AddMassStabTerms(rMassMatrix, Density, AdvVel, TauOne, Ng, DN_DX, GaussWeight);
533  }
534 */
535 
536  this->AddMassStabTerms(rMassMatrix, Density, AdvVel, TauOne, N, DN_DX, Area);
537 
538 //Z
539  }
540  }
541 //G
542  virtual void CalculateLaplacianMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo)
543  {
544  const unsigned int LocalSize = TDim * TNumNodes;
545 
546  // Resize and set to zero
547  if (rMassMatrix.size1() != LocalSize)
548  rMassMatrix.resize(LocalSize, LocalSize, false);
549 
550  rMassMatrix = ZeroMatrix(LocalSize, LocalSize);
551 
552  // Get the element's geometric parameters
553  double Area;
556  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
557 
558  // Calculate this element's fluid properties
559 
560  // Add 'classical' mass matrix (lumped)
561  double Coeff = Area / TNumNodes; //Optimize!
562 
563  this->CalculateLaplacianLumpedMassMatrix(rMassMatrix, Coeff);
564  }
565 //Z
566 
568 
576  virtual void CalculateLocalVelocityContribution(MatrixType& rDampingMatrix,
577  VectorType& rRightHandSideVector,
578  const ProcessInfo& rCurrentProcessInfo) override
579  {
580  const unsigned int LocalSize = (TDim + 1) * TNumNodes;
581 
582  // Resize and set to zero the matrix
583  // Note that we don't clean the RHS because it will already contain body force (and stabilization) contributions
584  if (rDampingMatrix.size1() != LocalSize)
585  rDampingMatrix.resize(LocalSize, LocalSize, false);
586 
587  noalias(rDampingMatrix) = ZeroMatrix(LocalSize, LocalSize);
588 
589  // Get this element's geometric properties
590  double Area;
593  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
594 
595  // Calculate this element's fluid properties
596  double Density, KinViscosity;
597  this->EvaluateInPoint(Density, DENSITY, N);
598  this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
599 
600  double Viscosity;
601  this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
602 
603  // Get Advective velocity
604  array_1d<double, 3 > AdvVel;
605  this->GetAdvectiveVel(AdvVel, N);
606 
607  // Calculate stabilization parameters
608  double TauOne, TauTwo;
609  this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
610 
611  /* if(this->GetValue(TRACK_SUBSCALES)==1)
612  {
613  const double DeltaTime = rCurrentProcessInfo[DELTA_TIME];
614  this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo, N, DN_DX, Area,DeltaTime);
615  }
616  else
617  {*/
618 //G
619  this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo, N, DN_DX, Area);
620 /*
621 
622  MatrixType NContainer;
623  ShapeFunctionDerivativesArrayType DN_DXContainer;
624  VectorType GaussWeights;
625  this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
626  const SizeType NumGauss = NContainer.size1();
627 
628  for (SizeType g = 0; g < NumGauss; g++)
629  {
630  const double GaussWeight = GaussWeights[g];
631  const ShapeFunctionsType& Ng = row(NContainer, g);
632  this->GetAdvectiveVel(AdvVel, Ng);
633  this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo, Ng, DN_DX, GaussWeight);
634  }
635 
636 */
637 
638 //Z
639 
640 // this->ModulatedGradientDiffusion(rDampingMatrix,DN_DX,Density*Area);
641 // }
642 
643  // Now calculate an additional contribution to the residual: r -= rDampingMatrix * (u,p)
644  VectorType U = ZeroVector(LocalSize);
645  int LocalIndex = 0;
646 
647  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
648  {
649  array_1d< double, 3 > & rVel = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY);
650  for (unsigned int d = 0; d < TDim; ++d) // Velocity Dofs
651  {
652  U[LocalIndex] = rVel[d];
653  ++LocalIndex;
654  }
655  U[LocalIndex] = this->GetGeometry()[iNode].FastGetSolutionStepValue(PRESSURE); // Pressure Dof
656  ++LocalIndex;
657  }
658 
659  noalias(rRightHandSideVector) -= prod(rDampingMatrix, U);
660  }
661 
662  virtual void FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override
663  {
664 // if(this->GetValue(TRACK_SUBSCALES) == 1)
665 // {
666 // // Get this element's geometric properties
667 // double Area;
668 // array_1d<double, TNumNodes> N;
669 // BoundedMatrix<double, TNumNodes, TDim> DN_DX;
670 // GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
671 //
672 // // Calculate this element's fluid properties
673 // double Density, KinViscosity;
674 // this->EvaluateInPoint(Density, DENSITY, N);
675 // this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
676 //
677 // double Viscosity;
678 // this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
679 //
680 // // Get Advective velocity
681 // array_1d<double, 3 > AdvVel;
682 // this->GetAdvectiveVel(AdvVel, N);
683 //
684 // // Calculate stabilization parameters
685 // double StaticTauOne;
686 // this->CalculateStaticTau(StaticTauOne,AdvVel,Area,Viscosity);
687 //
688 // array_1d<double,3> ElementalMomRes(3,0.0);
689 //
690 // if ( rCurrentProcessInfo[OSS_SWITCH] != 1 ) // ASGS
691 // {
692 // this->ASGSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);
693 // }
694 // else // OSS
695 // {
696 // this->OSSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);;
697 // }
698 //
699 // // Update subscale term
700 // const double DeltaTime = rCurrentProcessInfo.GetValue(DELTA_TIME);
701 // array_1d<double,3>& OldSubscaleVel = this->GetValue(SUBSCALE_VELOCITY);
702 // array_1d<double,3> Tmp = ( 1.0/( 1.0/DeltaTime + 1.0/StaticTauOne)) * (Density*OldSubscaleVel/DeltaTime + ElementalMomRes);
703 // OldSubscaleVel = Tmp;
704 // }
705  }
706 
707  void FinalizeSolutionStep(const ProcessInfo &rCurrentProcessInfo) override {
708  for(unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
709  {
710  this->GetGeometry()[iNode].SetLock();
711  this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION_OLD) = this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION);
712  this->GetGeometry()[iNode].UnSetLock();
713  }
714  }
715 
717 
729  virtual void Calculate(const Variable<double>& rVariable,
730  double& rOutput,
731  const ProcessInfo& rCurrentProcessInfo) override
732  {
733  if (rVariable == ERROR_RATIO)
734  {
735  // Get the element's geometric parameters
736  double Area;
739  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
740 
741  // Calculate this element's fluid properties
742  double Density, KinViscosity;
743  this->EvaluateInPoint(Density, DENSITY, N);
744  this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
745 
746  double Viscosity;
747  this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
748 
749  // Get Advective velocity
750  array_1d<double, 3 > AdvVel;
751  this->GetAdvectiveVel(AdvVel, N);
752 
753  // Output container
754  array_1d< double, 3 > ElementalMomRes(3, 0.0);
755 
756  // Calculate stabilization parameter. Note that to estimate the subscale velocity, the dynamic coefficient in TauOne is assumed zero.
757  double TauOne;
758  this->CalculateStaticTau(TauOne, AdvVel, Area,Density, Viscosity);
759 
760  if ( rCurrentProcessInfo[OSS_SWITCH] != 1 ) // ASGS
761  {
762 //G
763  this->ASGSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);
764 
765 
766 /*
767  MatrixType NContainer;
768  ShapeFunctionDerivativesArrayType DN_DXContainer;
769  VectorType GaussWeights;
770  this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
771  const SizeType NumGauss = NContainer.size1();
772 
773  for (SizeType g = 0; g < NumGauss; g++){
774  const double GaussWeight = GaussWeights[g];
775  const ShapeFunctionsType& Ng = row(NContainer, g);
776  this->GetAdvectiveVel(AdvVel, Ng);
777  this->ASGSMomResidual(AdvVel, Density, ElementalMomRes, Ng, DN_DX, GaussWeight);
778  }
779 
780 */
781 //Z
782  ElementalMomRes *= TauOne;
783  }
784  else // OSS
785  {
786  this->OSSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);;
787  ElementalMomRes *= TauOne;
788  }
789 
790  // Error estimation ( ||U'|| / ||Uh_gauss|| ), taking ||U'|| = TauOne ||MomRes||
791  double ErrorRatio(0.0);//, UNorm(0.0);
792 // array_1d< double, 3 > UGauss(3, 0.0);
793 // this->AddPointContribution(UGauss, VELOCITY, N);
794 
795  for (unsigned int i = 0; i < TDim; ++i)
796  {
797  ErrorRatio += ElementalMomRes[i] * ElementalMomRes[i];
798 // UNorm += UGauss[i] * UGauss[i];
799  }
800  ErrorRatio = sqrt(ErrorRatio); // / UNorm);
801  ErrorRatio /= Density;
802  this->SetValue(ERROR_RATIO, ErrorRatio);
803  rOutput = ErrorRatio;
804  }
805  else if (rVariable == NODAL_AREA)
806  {
807  // Get the element's geometric parameters
808  double Area;
811  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
812 
813  // Carefully write results to nodal variables, to avoid parallelism problems
814  for (unsigned int i = 0; i < TNumNodes; ++i)
815  {
816  this->GetGeometry()[i].SetLock(); // So it is safe to write in the node in OpenMP
817  this->GetGeometry()[i].FastGetSolutionStepValue(NODAL_AREA) += Area * N[i];
818  this->GetGeometry()[i].UnSetLock(); // Free the node for other threads
819  }
820  }
821  }
822 
824 
835  virtual void Calculate(const Variable<array_1d<double, 3 > >& rVariable,
836  array_1d<double, 3 > & rOutput,
837  const ProcessInfo& rCurrentProcessInfo) override
838  {
839  if (rVariable == ADVPROJ) // Compute residual projections for OSS
840  {
841  // Get the element's geometric parameters
842  double Area;
845  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
846 
847  // Calculate this element's fluid properties
848  double Density;
849  this->EvaluateInPoint(Density, DENSITY, N);
850 
851  // Get Advective velocity
852  array_1d<double, 3 > AdvVel;
853  this->GetAdvectiveVel(AdvVel, N);
854 
855  // Output containers
856  array_1d< double, 3 > ElementalMomRes(3, 0.0);
857  double ElementalMassRes(0);
858 
859  this->AddProjectionResidualContribution(AdvVel, Density, ElementalMomRes, ElementalMassRes, N, DN_DX, Area);
860 
861  if (rCurrentProcessInfo[OSS_SWITCH] == 1)
862  {
863  // Carefully write results to nodal variables, to avoid parallelism problems
864  for (unsigned int i = 0; i < TNumNodes; ++i)
865  {
866  this->GetGeometry()[i].SetLock(); // So it is safe to write in the node in OpenMP
867  array_1d< double, 3 > & rAdvProj = this->GetGeometry()[i].FastGetSolutionStepValue(ADVPROJ);
868  for (unsigned int d = 0; d < TDim; ++d)
869  rAdvProj[d] += N[i] * ElementalMomRes[d];
870 
871  this->GetGeometry()[i].FastGetSolutionStepValue(DIVPROJ) += N[i] * ElementalMassRes;
872  this->GetGeometry()[i].FastGetSolutionStepValue(NODAL_AREA) += Area * N[i];
873  this->GetGeometry()[i].UnSetLock(); // Free the node for other threads
874  }
875  }
876 
878  rOutput = ElementalMomRes;
879  }
880  else if (rVariable == SUBSCALE_VELOCITY)
881  {
882  // Get the element's geometric parameters
883  double Area;
886  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
887 
888  // Calculate this element's fluid properties
889  double Density;
890  this->EvaluateInPoint(Density, DENSITY, N);
891 
892  // Get Advective velocity
893  array_1d<double, 3 > AdvVel;
894  this->GetAdvectiveVel(AdvVel, N);
895 
896  // Output containers
897  array_1d< double, 3 > ElementalMomRes(3,0.0);
898  double ElementalMassRes(0.0);
899 
900  this->AddProjectionResidualContribution(AdvVel, Density, ElementalMomRes, ElementalMassRes, N, DN_DX, Area);
901 
902  if (rCurrentProcessInfo[OSS_SWITCH] == 1)
903  {
904  /* Projections of the elemental residual are computed with
905  * Newton-Raphson iterations of type M(lumped) dx = ElemRes - M(consistent) * x
906  */
907  const double Weight = ConsistentMassCoef(Area); // Consistent mass matrix is Weigth * ( Ones(TNumNodes,TNumNodes) + Identity(TNumNodes,TNumNodes) )
908  // Carefully write results to nodal variables, to avoid parallelism problems
909  for (unsigned int i = 0; i < TNumNodes; ++i)
910  {
911  this->GetGeometry()[i].SetLock(); // So it is safe to write in the node in OpenMP
912 
913  // Add elemental residual to RHS
914  array_1d< double, 3 > & rMomRHS = this->GetGeometry()[i].GetValue(ADVPROJ);
915  double& rMassRHS = this->GetGeometry()[i].GetValue(DIVPROJ);
916  for (unsigned int d = 0; d < TDim; ++d)
917  rMomRHS[d] += N[i] * ElementalMomRes[d];
918 
919  rMassRHS += N[i] * ElementalMassRes;
920 
921  // Write nodal area
922  this->GetGeometry()[i].FastGetSolutionStepValue(NODAL_AREA) += Area * N[i];
923 
924  // Substract M(consistent)*x(i-1) from RHS
925  for(unsigned int j = 0; j < TNumNodes; ++j) // RHS -= Weigth * Ones(TNumNodes,TNumNodes) * x(i-1)
926  {
927  for(unsigned int d = 0; d < TDim; ++d)
928  rMomRHS[d] -= Weight * this->GetGeometry()[j].FastGetSolutionStepValue(ADVPROJ)[d];
929  rMassRHS -= Weight * this->GetGeometry()[j].FastGetSolutionStepValue(DIVPROJ);
930  }
931  for(unsigned int d = 0; d < TDim; ++d) // RHS -= Weigth * Identity(TNumNodes,TNumNodes) * x(i-1)
932  rMomRHS[d] -= Weight * this->GetGeometry()[i].FastGetSolutionStepValue(ADVPROJ)[d];
933  rMassRHS -= Weight * this->GetGeometry()[i].FastGetSolutionStepValue(DIVPROJ);
934 
935  this->GetGeometry()[i].UnSetLock(); // Free the node for other threads
936  }
937  }
938 
940  rOutput = ElementalMomRes;
941  }
942  }
943 
944  // The following methods have different implementations depending on TDim
946 
952  virtual void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo) const override;
953 
955 
959  virtual void GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo& rCurrentProcessInfo) const override;
960 
962 
966  virtual void GetFirstDerivativesVector(Vector& Values, int Step = 0) const override;
967 
969 
973  virtual void GetSecondDerivativesVector(Vector& Values, int Step = 0) const override;
974 
976 
986  std::vector<array_1d<double, 3 > >& rOutput,
987  const ProcessInfo& rCurrentProcessInfo) override;
988 
990 
1001  virtual void CalculateOnIntegrationPoints(const Variable<double>& rVariable,
1002  std::vector<double>& rValues,
1003  const ProcessInfo& rCurrentProcessInfo) override
1004  {
1005  if (rVariable == TAUONE || rVariable == TAUTWO || rVariable == MU)
1006  {
1007  double TauOne, TauTwo;
1008  double Area;
1011  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
1012 
1013  array_1d<double, 3 > AdvVel;
1014  this->GetAdvectiveVel(AdvVel, N);
1015 
1016  double Density,KinViscosity;
1017  this->EvaluateInPoint(Density, DENSITY, N);
1018  this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
1019 
1020  double Viscosity;
1021  this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
1022 
1023  this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
1024 
1025  rValues.resize(1, false);
1026  if (rVariable == TAUONE)
1027  {
1028  rValues[0] = TauOne;
1029  }
1030  else if (rVariable == TAUTWO)
1031  {
1032  rValues[0] = TauTwo;
1033  }
1034  else if (rVariable == MU)
1035  {
1036  rValues[0] = Density * Viscosity;
1037  }
1038  }
1039  else if(rVariable == SUBSCALE_PRESSURE)
1040  {
1041  double TauOne, TauTwo;
1042  double Area;
1045  GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
1046 
1047  array_1d<double, 3 > AdvVel;
1048  this->GetAdvectiveVel(AdvVel, N);
1049 
1050  double Density,KinViscosity;
1051  this->EvaluateInPoint(Density, DENSITY, N);
1052  this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
1053 
1054  double Viscosity;
1055  this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
1056 
1057  this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
1058 
1059  double DivU = 0.0;
1060  for(unsigned int i=0; i < TNumNodes; i++)
1061  {
1062  for(unsigned int d = 0; d < TDim; d++)
1063  DivU -= DN_DX(i,d) * this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY)[d];
1064  }
1065 
1066  rValues.resize(1, false);
1067 
1068  rValues[0] = TauTwo * DivU;// *Density?? decide on criteria and use the same for SUBSCALE_VELOCITY
1069 
1070  if(rCurrentProcessInfo[OSS_SWITCH]==1)
1071  {
1072  double Proj = 0.0;
1073  for(unsigned int i=0; i < TNumNodes; i++)
1074  {
1075  Proj += N[i]*this->GetGeometry()[i].FastGetSolutionStepValue(DIVPROJ);
1076  }
1077  rValues[0] -= TauTwo*Proj;
1078  }
1079  }
1080  else if (rVariable == NODAL_AREA && TDim == 3)
1081  {
1082  MatrixType J = ZeroMatrix(3,3);
1083  const array_1d<double,3>& X0 = this->GetGeometry()[0].Coordinates();
1084  const array_1d<double,3>& X1 = this->GetGeometry()[1].Coordinates();
1085  const array_1d<double,3>& X2 = this->GetGeometry()[2].Coordinates();
1086  const array_1d<double,3>& X3 = this->GetGeometry()[3].Coordinates();
1087 
1088  J(0,0) = X1[0]-X0[0];
1089  J(0,1) = X2[0]-X0[0];
1090  J(0,2) = X3[0]-X0[0];
1091  J(1,0) = X1[1]-X0[1];
1092  J(1,1) = X2[1]-X0[1];
1093  J(1,2) = X3[1]-X0[1];
1094  J(2,0) = X1[2]-X0[2];
1095  J(2,1) = X2[2]-X0[2];
1096  J(2,2) = X3[2]-X0[2];
1097 
1098  double DetJ = J(0,0)*( J(1,1)*J(2,2) - J(1,2)*J(2,1) ) + J(0,1)*( J(1,2)*J(2,0) - J(1,0)*J(2,2) ) + J(0,2)*( J(1,0)*J(2,1) - J(1,1)*J(2,0) );
1099  rValues.resize(1, false);
1100  rValues[0] = DetJ;
1101  }
1102  else // Default behaviour (returns elemental data)
1103  {
1104  rValues.resize(1, false);
1105  /*
1106  The cast is done to avoid modification of the element's data. Data modification
1107  would happen if rVariable is not stored now (would initialize a pointer to &rVariable
1108  with associated value of 0.0). This is catastrophic if the variable referenced
1109  goes out of scope.
1110  */
1111  const MonolithicDEMCoupled<TDim, TNumNodes>* const_this = static_cast<const MonolithicDEMCoupled<TDim, TNumNodes>*> (this);
1112  rValues[0] = const_this->GetValue(rVariable);
1113  }
1114  }
1115 
1117  std::vector<array_1d<double, 6 > >& rValues,
1118  const ProcessInfo& rCurrentProcessInfo) override
1119  {}
1120 
1121  virtual void CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
1122  std::vector<Vector>& rValues,
1123  const ProcessInfo& rCurrentProcessInfo) override
1124  {}
1125 
1126  virtual void CalculateOnIntegrationPoints(const Variable<Matrix>& rVariable,
1127  std::vector<Matrix>& rValues,
1128  const ProcessInfo& rCurrentProcessInfo) override
1129  {}
1130 
1134 
1138 
1140 
1148  virtual int Check(const ProcessInfo& rCurrentProcessInfo) const override
1149  {
1150  KRATOS_TRY
1151 
1152  // Perform basic element checks
1153  int ErrorCode = Kratos::Element::Check(rCurrentProcessInfo);
1154  if(ErrorCode != 0) return ErrorCode;
1155 
1156  // Additional variables, only required to print results:
1157  // SUBSCALE_VELOCITY, SUBSCALE_PRESSURE, TAUONE, TAUTWO, MU, VORTICITY.
1158 
1159  // Checks on nodes
1160 
1161  // Check that the element's nodes contain all required SolutionStepData and Degrees of freedom
1162  for(unsigned int i=0; i<this->GetGeometry().size(); ++i) {
1163  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(VELOCITY, this->GetGeometry()[i])
1164  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(PRESSURE, this->GetGeometry()[i])
1165  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(MESH_VELOCITY, this->GetGeometry()[i])
1166  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(ACCELERATION, this->GetGeometry()[i])
1167 
1168  KRATOS_CHECK_DOF_IN_NODE(VELOCITY_X, this->GetGeometry()[i])
1169  KRATOS_CHECK_DOF_IN_NODE(VELOCITY_Y, this->GetGeometry()[i])
1170  KRATOS_CHECK_DOF_IN_NODE(VELOCITY_Z, this->GetGeometry()[i])
1171  KRATOS_CHECK_DOF_IN_NODE(PRESSURE, this->GetGeometry()[i])
1172  }
1173  // Not checking OSS related variables NODAL_AREA, ADVPROJ, DIVPROJ, which are only required as SolutionStepData if OSS_SWITCH == 1
1174 
1175  // If this is a 2D problem, check that nodes are in XY plane
1176  if (this->GetGeometry().WorkingSpaceDimension() == 2) {
1177  for (unsigned int i=0; i<this->GetGeometry().size(); ++i) {
1178  KRATOS_ERROR_IF(this->GetGeometry()[i].Z() != 0.0) << "Node with non-zero Z coordinate found. Id: "<< this->GetGeometry()[i].Id() << std::endl;
1179  }
1180  }
1181 
1182  return 0;
1183 
1184  KRATOS_CATCH("");
1185  }
1186 
1190 
1191 
1195 
1197  virtual std::string Info() const override
1198  {
1199  std::stringstream buffer;
1200  buffer << "MonolithicDEMCoupled #" << Id();
1201  return buffer.str();
1202  }
1203 
1205  virtual void PrintInfo(std::ostream& rOStream) const override
1206  {
1207  rOStream << "MonolithicDEMCoupled" << TDim << "D";
1208  }
1209 
1210 // /// Print object's data.
1211 // virtual void PrintData(std::ostream& rOStream) const;
1212 
1216 
1217 
1219 
1220 protected:
1223 
1224 //G
1226  void CalculateWeights(ShapeFunctionDerivativesArrayType& rDN_DX, Matrix& rNContainer, Vector& rGaussWeights);
1227 //Z
1228 
1232 
1233 
1237 
1238 
1242 
1244 
1256  virtual void CalculateTau(double& TauOne,
1257  double& TauTwo,
1258  const array_1d< double, 3 > & rAdvVel,
1259  const double Area,
1260  const double Density,
1261  const double KinViscosity,
1262  const ProcessInfo& rCurrentProcessInfo)
1263  {
1264  // Compute mean advective velocity norm
1265  double AdvVelNorm = 0.0;
1266  for (unsigned int d = 0; d < TDim; ++d)
1267  AdvVelNorm += rAdvVel[d] * rAdvVel[d];
1268 
1269  AdvVelNorm = sqrt(AdvVelNorm);
1270 
1271  const double Element_Size = this->ElementSize(Area);
1272 //G
1273  // TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 4 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
1274  TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 5.6666666666 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
1275 //Z
1276  TauTwo = Density * (KinViscosity + 0.5 * Element_Size * AdvVelNorm);
1277  //TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 12.0 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
1278  //TauTwo = Density * (KinViscosity + Element_Size * AdvVelNorm / 6.0);
1279 
1280 
1281  }
1282 
1284 
1294  virtual void CalculateStaticTau(double& TauOne,
1295  const array_1d< double, 3 > & rAdvVel,
1296  const double Area,
1297  const double Density,
1298  const double KinViscosity)
1299  {
1300  // Compute mean advective velocity norm
1301  double AdvVelNorm = 0.0;
1302  for (unsigned int d = 0; d < TDim; ++d)
1303  AdvVelNorm += rAdvVel[d] * rAdvVel[d];
1304 
1305  AdvVelNorm = sqrt(AdvVelNorm);
1306 
1307  const double Element_Size = this->ElementSize(Area);
1308 
1309  TauOne = 1.0 / (Density*(4.0 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size));
1310  }
1311 
1313  virtual void AddMomentumRHS(VectorType& F,
1314  const double Density,
1315  const array_1d<double, TNumNodes>& rShapeFunc,
1316  const double Weight)
1317  {
1318  double Coef = Density * Weight;
1319 
1320  array_1d<double, 3 > BodyForce(3, 0.0);
1321  this->EvaluateInPoint(BodyForce, BODY_FORCE, rShapeFunc);
1322 
1323  // Add the results to the velocity components (Local Dofs are vx, vy, [vz,] p for each node)
1324  int LocalIndex = 0;
1325 
1326  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1327  {
1328  for (unsigned int d = 0; d < TDim; ++d)
1329  {
1330  F[LocalIndex++] += Coef * rShapeFunc[iNode] * BodyForce[d];
1331  }
1332  ++LocalIndex; // Skip pressure Dof
1333  }
1334  }
1335 
1337  const BoundedMatrix<double, TNumNodes, TDim>& rShapeDeriv,
1338  const double Weight)
1339  {
1340  double Coef = Weight;
1341  array_1d<double, 3 > Velocity(3, 0.0);
1342 
1343  int LocalIndex = 0;
1344  int LocalNodalIndex = 0;
1345 
1346  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1347  {
1348  Velocity = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY);
1349  for (unsigned int d = 0; d < TDim; ++d)
1350  {
1351  F[LocalIndex++] -= Coef * rShapeDeriv(LocalNodalIndex, d) * Velocity[d] * rShapeDeriv(iNode, d);
1352  }
1353  LocalNodalIndex++;
1354  }
1355  }
1356 //G
1357  virtual void AddMassRHS(VectorType& F,
1358  const double Density,
1359  const array_1d<double, TNumNodes>& rShapeFunc,
1360  const double Weight,
1361  const std::vector<double>& TimeSchemeWeights,
1362  const double& DeltaTime)
1363  {
1364  double FluidFractionRate = 0.0;
1365  this->EvaluateTimeDerivativeInPoint(FluidFractionRate, FLUID_FRACTION_RATE, rShapeFunc, DeltaTime, TimeSchemeWeights);
1366  //this->EvaluateInPoint(FluidFractionRate, FLUID_FRACTION_RATE, rShapeFunc);
1367  // Add the results to the pressure components (Local Dofs are vx, vy, [vz,] p for each node)
1368  int LocalIndex = TDim;
1369 
1370  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
1371  F[LocalIndex] -= Weight * rShapeFunc[iNode] * FluidFractionRate;
1372  LocalIndex += TDim + 1;
1373  }
1374 
1375  }
1376 //Z
1377 
1379  virtual void AddProjectionToRHS(VectorType& RHS,
1380  const array_1d<double, 3 > & rAdvVel,
1381  const double Density,
1382  const double TauOne,
1383  const double TauTwo,
1384  const array_1d<double, TNumNodes>& rShapeFunc,
1385  const BoundedMatrix<double, TNumNodes, TDim>& rShapeDeriv,
1386  const double Weight,
1387  const double DeltaTime = 1.0)
1388  {
1389  const unsigned int BlockSize = TDim + 1;
1390 
1392  this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1393  array_1d<double,3> MomProj(3,0.0);
1394  double DivProj = 0.0;
1395  this->EvaluateInPoint(MomProj,ADVPROJ,rShapeFunc);
1396  this->EvaluateInPoint(DivProj,DIVPROJ,rShapeFunc);
1397 
1398  MomProj *= TauOne;
1399  DivProj *= TauTwo;
1400 
1401  unsigned int FirstRow = 0;
1402 
1403  for (unsigned int i = 0; i < TNumNodes; i++)
1404  {
1405 //G
1406  double FluidFraction = this->GetGeometry()[i].FastGetSolutionStepValue(FLUID_FRACTION);
1407  array_1d<double,3> FluidFractionGradient(3,0.0);
1408  FluidFractionGradient[0] += rShapeDeriv(i, 0) * FluidFraction;
1409  FluidFractionGradient[1] += rShapeDeriv(i, 1) * FluidFraction;
1410  FluidFractionGradient[2] += rShapeDeriv(i, 2) * FluidFraction;
1411 //Z
1412 
1413  for (unsigned int d = 0; d < TDim; d++)
1414  {
1415 //G
1416 // RHS[FirstRow+d] -= Weight * (Density * AGradN[i] * MomProj[d] + rShapeDeriv(i,d) * DivProj); // TauOne * ( a * Grad(v) ) * MomProjection + TauTwo * Div(v) * MassProjection
1417  RHS[FirstRow+d] -= Weight * (Density * AGradN[i] * MomProj[d] + (FluidFraction * rShapeDeriv(i,d) + FluidFractionGradient[d] * rShapeFunc[i]) * DivProj); // TauOne * (a * Grad(v)) * MomProjection + TauTwo * Div(v) * MassProjection
1418 //Z
1419  RHS[FirstRow+TDim] -= Weight * rShapeDeriv(i,d) * MomProj[d]; // TauOne * Grad(q) * MomProjection
1420  }
1421  FirstRow += BlockSize;
1422  }
1423  }
1424 
1426 
1433  const double Mass)
1434  {
1435  unsigned int DofIndex = 0;
1436  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1437  {
1438  for (unsigned int d = 0; d < TDim; ++d)
1439  {
1440  rLHSMatrix(DofIndex, DofIndex) += Mass;
1441  ++DofIndex;
1442  }
1443  ++DofIndex; // Skip pressure Dof
1444  }
1445  }
1446 
1447 //G
1449  const double Mass)
1450  {
1451  unsigned int DofIndex = 0;
1452  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
1453  {
1454  for (unsigned int d = 0; d < TDim; ++d)
1455  {
1456  rLHSMatrix(DofIndex, DofIndex) += Mass;
1457  ++DofIndex;
1458  }
1459  }
1460  }
1461 //Z
1462 
1464  const array_1d<double,TNumNodes>& rShapeFunc,
1465  const double Density,
1466  const double Weight)
1467  {
1468  const unsigned int BlockSize = TDim + 1;
1469 
1470  double Coef = Density * Weight;
1471  unsigned int FirstRow(0), FirstCol(0);
1472  double K; // Temporary results
1473 
1474  // Note: Dof order is (vx,vy,[vz,]p) for each node
1475  for (unsigned int i = 0; i < TNumNodes; ++i)
1476  {
1477  // Loop over columns
1478  for (unsigned int j = 0; j < TNumNodes; ++j)
1479  {
1480  // Delta(u) * TauOne * [ AdvVel * Grad(v) ] in velocity block
1481  K = Coef * rShapeFunc[i] * rShapeFunc[j];
1482 
1483  for (unsigned int d = 0; d < TDim; ++d) // iterate over dimensions for velocity Dofs in this node combination
1484  {
1485  rLHSMatrix(FirstRow + d, FirstCol + d) += K;
1486  }
1487  // Update column index
1488  FirstCol += BlockSize;
1489  }
1490  // Update matrix indices
1491  FirstRow += BlockSize;
1492  FirstCol = 0;
1493  }
1494  }
1495 
1496 
1498 
1510  void AddMassStabTerms(MatrixType& rLHSMatrix,
1511  const double Density,
1512  const array_1d<double, 3 > & rAdvVel,
1513  const double TauOne,
1514  const array_1d<double, TNumNodes>& rShapeFunc,
1515  const BoundedMatrix<double, TNumNodes, TDim>& rShapeDeriv,
1516  const double Weight)
1517  {
1518  const unsigned int BlockSize = TDim + 1;
1519 
1520  double Coef = Weight * TauOne;
1521  unsigned int FirstRow(0), FirstCol(0);
1522  double K; // Temporary results
1523 
1524  // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
1525  array_1d<double, TNumNodes> AGradN, AGradNMod;
1526  this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1527 //G
1528  double AdvVelDiv = 0.0;
1529  this->GetAdvectiveVelDivergence(AdvVelDiv, rShapeDeriv);
1530  //this->GetModifiedConvectionOperator(AGradNMod, rAdvVel, AdvVelDiv, rShapeFunc, rShapeDeriv); // Get a * grad(Ni) + div(a) * Ni
1531  double FluidFraction;
1532  this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
1533 //Z
1534 
1535  // Note: Dof order is (vx,vy,[vz,]p) for each node
1536  for (unsigned int i = 0; i < TNumNodes; ++i)
1537  {
1538  // Loop over columns
1539  for (unsigned int j = 0; j < TNumNodes; ++j)
1540  {
1541  // Delta(u) * TauOne * [ AdvVel * Grad(v) ] in velocity block
1542 //G
1543  K = Coef * Density * AGradN[i] * Density * rShapeFunc[j];
1544  //K = Coef * Density * AGradNMod[i] * Density * rShapeFunc[j];
1545 //Z
1546 
1547  for (unsigned int d = 0; d < TDim; ++d) // iterate over dimensions for velocity Dofs in this node combination
1548  {
1549  rLHSMatrix(FirstRow + d, FirstCol + d) += K;
1550  // Delta(u) * TauOne * Grad(q) in q * Div(u) block
1551 //G
1552  // rLHSMatrix(FirstRow + TDim, FirstCol + d) += Coef * Density * rShapeDeriv(i, d) * rShapeFunc[j];
1553  rLHSMatrix(FirstRow + TDim, FirstCol + d) += Coef * Density * FluidFraction * rShapeDeriv(i, d) * rShapeFunc[j]; // Delta(u) * TauOne * alpha * Grad(q)
1554 //Z
1555  }
1556  // Update column index
1557  FirstCol += BlockSize;
1558  }
1559  // Update matrix indices
1560  FirstRow += BlockSize;
1561  FirstCol = 0;
1562  }
1563  }
1564 
1567  VectorType& rDampRHS,
1568  const double Density,
1569  const double Viscosity,
1570  const array_1d< double, 3 > & rAdvVel,
1571  const double TauOne,
1572  const double TauTwo,
1573  const array_1d< double, TNumNodes >& rShapeFunc,
1574  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
1575  const double Weight)
1576  {
1577  const unsigned int BlockSize = TDim + 1;
1578 
1579  // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
1580  array_1d<double, TNumNodes> AGradN, AGradNMod;
1581  this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1582 //G
1583  double AdvVelDiv = 0.0;
1584  this->GetAdvectiveVelDivergence(AdvVelDiv, rShapeDeriv);
1585  this->GetModifiedConvectionOperator(AGradNMod, rAdvVel, AdvVelDiv, rShapeFunc, rShapeDeriv); // Get a * grad(Ni) + div(a) * Ni
1586 //Z
1587  // Build the local matrix and RHS
1588  unsigned int FirstRow(0), FirstCol(0); // position of the first term of the local matrix that corresponds to each node combination
1589  double K, G, PDivV, L, qF; // Temporary results
1590 
1591  array_1d<double,3> BodyForce(3,0.0);
1592  this->EvaluateInPoint(BodyForce,BODY_FORCE,rShapeFunc);
1593  BodyForce *= Density;
1594 //G
1595  double PAlphaDivV, GAlpha, FluidFraction, FluidFractionRate;
1596  array_1d<double,3> FluidFractionGradient(3,0.0);
1597  this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
1598  this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
1599 
1600  for (unsigned int i = 0; i < TNumNodes; ++i) {
1601  this->GetGeometry()[i].FastGetSolutionStepValue(FLUID_FRACTION_GRADIENT) = FluidFractionGradient;
1602  }
1603 
1604  this->EvaluateInPoint(FluidFractionRate,FLUID_FRACTION_RATE,rShapeFunc);
1605 
1606  const double EpsilonInside = false;
1607 
1608  if (EpsilonInside){
1609  array_1d<double,3> rGradEpsOverEps;
1610  rGradEpsOverEps = 1.0 / FluidFraction * FluidFractionGradient ;
1611  }
1612 
1613 //Z
1614  for (unsigned int i = 0; i < TNumNodes; ++i) // iterate over rows
1615  {
1616  for (unsigned int j = 0; j < TNumNodes; ++j) // iterate over columns
1617  {
1618  // Calculate the part of the contributions that is constant for each node combination
1619 
1620  // Velocity block
1621  K = Density * rShapeFunc[i] * AGradN[j]; // Convective term: v * ( a * Grad(u) )
1622  //K = 0.5 * Density * (rShapeFunc[i] * AGradN[j] - AGradN[i] * rShapeFunc[j]); // Skew-symmetric convective term 1/2( v*grad(u)*u - grad(v) uu )
1623 //G
1624  K += TauOne * Density * AGradN[i] * Density * AGradN[j]; // Stabilization: (a * Grad(v)) * TauOne * (a * Grad(u))
1625  //K += TauOne * Density * AGradNMod[i] * Density * AGradN[j]; // Stabilization: (a * Grad(v) + Div(a) * v) * TauOne * (a * Grad(u))
1626 //Z
1627  K *= Weight;
1628 
1629  // q-p stabilization block (reset result)
1630  L = 0;
1631 
1632  for (unsigned int m = 0; m < TDim; ++m) // iterate over v components (vx,vy[,vz])
1633  {
1634  // Velocity block
1635  //K += Weight * Density * Viscosity * rShapeDeriv(i, m) * rShapeDeriv(j, m); // Diffusive term: Viscosity * Grad(v) * Grad(u)
1636 
1637  // v * Grad(p) block
1638 //G
1639  G = TauOne * Density * AGradN[i] * rShapeDeriv(j, m); // Stabilization: (a * Grad(v)) * TauOne * Grad(p)
1640  //G = TauOne * Density * AGradNMod[i] * rShapeDeriv(j, m); // Stabilization: (a * Grad(v) + Div(a) * v) * TauOne * Grad(p)
1641 
1642  GAlpha = TauOne * Density * AGradN[i] * (FluidFraction * rShapeDeriv(j, m)); // Stabilization: (a * Grad(u)) * TauOne * (alpha * Grad(q))
1643 
1644  PAlphaDivV = (FluidFraction * rShapeDeriv(i, m) + FluidFractionGradient[m] * rShapeFunc[i]) * rShapeFunc[j]; // alpha * q * Div(u) + q * Grad(alpha) * u
1645 //Z
1646  PDivV = rShapeDeriv(i, m) * rShapeFunc[j]; // Div(v) * p
1647  // Write v * Grad(p) component
1648  rDampingMatrix(FirstRow + m, FirstCol + TDim) += Weight * (G - PDivV);
1649  // Use symmetry to write the q * Div(u) component
1650 //G
1651  // rDampingMatrix(FirstCol + TDim, FirstRow + m) += Weight * (G + PDivV);
1652  rDampingMatrix(FirstCol + TDim, FirstRow + m) += Weight * (GAlpha + PAlphaDivV); // note that here PAlphaDivV includes G; do not look for it in GAlpha!
1653 
1654  // q-p stabilization block
1655  // L += rShapeDeriv(i, m) * rShapeDeriv(j, m); // Stabilization: Grad(q) * TauOne * Grad(p)
1656  L += FluidFraction * rShapeDeriv(i, m) * rShapeDeriv(j, m); // Stabilization: alpha * Grad(q) * TauOne * Grad(p)
1657 //Z
1658  for (unsigned int n = 0; n < TDim; ++n) // iterate over u components (ux,uy[,uz])
1659  {
1660  // Velocity block
1661 //G
1662  // rDampingMatrix(FirstRow + m, FirstCol + n) += Weight * TauTwo * rShapeDeriv(i, m) * rShapeDeriv(j, n); // Stabilization: Div(v) * TauTwo * Div(u)
1663  rDampingMatrix(FirstRow + m, FirstCol + n) += Weight * TauTwo * rShapeDeriv(i, m) * (FluidFraction * rShapeDeriv(j, n) + FluidFractionGradient[n] * rShapeFunc[j]); // Stabilization: Div(v) * TauTwo * (alpha * Div(u) + Grad(alpha) * u)
1664 //Z
1665  }
1666 
1667  }
1668 
1669  // Write remaining terms to velocity block
1670  for (unsigned int d = 0; d < TDim; ++d)
1671  rDampingMatrix(FirstRow + d, FirstCol + d) += K;
1672 
1673  // Write q-p stabilization block
1674  rDampingMatrix(FirstRow + TDim, FirstCol + TDim) += Weight * TauOne * L;
1675 
1676 
1677  // Update reference column index for next iteration
1678  FirstCol += BlockSize;
1679  }
1680 
1681  // Operate on RHS
1682  qF = 0.0;
1683 
1684  for (unsigned int d = 0; d < TDim; ++d)
1685  {
1686  //rDampRHS[FirstRow + d] += Weight * TauOne * Density * AGradN[i] * BodyForce[d]; // ( a * Grad(v) ) * TauOne * (Density * BodyForce)
1687 //G
1688 //GG rDampRHS[FirstRow + d] += Weight * (TauOne * Density * AGradNMod[i] * BodyForce[d] - TauTwo * rShapeDeriv(i, d) * FluidFractionRate); // ( a * Grad(v) ) * TauOne * (Density * BodyForce) + Div(v) * TauTwo * (- DAlphaDt)
1689  rDampRHS[FirstRow + d] += Weight * (TauOne * Density * AGradN[i] * BodyForce[d] - TauTwo * rShapeDeriv(i, d) * FluidFractionRate); // ( a * Grad(v) ) * TauOne * (Density * BodyForce) + Div(v) * TauTwo * (- DAlphaDt)
1690 
1691  // qF += rShapeDeriv(i, d) * BodyForce[d];
1692  qF += (FluidFraction * rShapeDeriv(i, d)) * BodyForce[d];
1693 //Z
1694  }
1695  rDampRHS[FirstRow + TDim] += Weight * TauOne * qF; // Grad(q) * TauOne * (Density * BodyForce)
1696 
1697  // Update reference indices
1698  FirstRow += BlockSize;
1699  FirstCol = 0;
1700  }
1701 
1702 // this->AddBTransCB(rDampingMatrix,rShapeDeriv,Viscosity*Coef);
1703  this->AddViscousTerm(rDampingMatrix,rShapeDeriv,Viscosity*Density*Weight);
1704  }
1705 
1706 /*
1708 // void AddIntegrationPointVelocityContribution(MatrixType& rDampingMatrix,
1709 // VectorType& rDampRHS,
1710 // const double Density,
1711 // const double Viscosity,
1712 // const array_1d< double, 3 > & rAdvVel,
1713 // const double TauOne,
1714 // const double TauTwo,
1715 // const array_1d< double, TNumNodes >& rShapeFunc,
1716 // const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
1717 // const double Weight,
1718 // const double DeltaTime)
1719 // {
1720 // const unsigned int BlockSize = TDim + 1;
1721 //
1722 // const double TauCoef = TauOne*TauTwo / DeltaTime;
1723 //
1724 // // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
1725 // array_1d<double, TNumNodes> AGradN;
1726 // this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1727 //
1728 // // Build the local matrix and RHS
1729 // unsigned int FirstRow(0), FirstCol(0); // position of the first term of the local matrix that corresponds to each node combination
1730 // double K, G, PDivV, L; // Temporary results
1731 // double Coef = Density * Weight;
1732 //
1733 // // Note that we iterate first over columns, then over rows to read the Body Force only once per node
1734 // for (unsigned int j = 0; j < TNumNodes; ++j) // iterate over colums
1735 // {
1736 // // Get Body Force
1737 // const array_1d<double, 3 > & rBodyForce = this->GetGeometry()[j].FastGetSolutionStepValue(BODY_FORCE);
1738 //
1739 // const array_1d<double,3>& OldVelocity = this->GetGeometry()[j].FastGetSolutionStepValue(VELOCITY,1);
1740 //
1741 // for (unsigned int i = 0; i < TNumNodes; ++i) // iterate over rows
1742 // {
1743 // // Calculate the part of the contributions that is constant for each node combination
1744 //
1745 // // Velocity block
1746 // K = rShapeFunc[i] * AGradN[j]; // Convective term: v * ( a * Grad(u) )
1747 // K += TauOne * AGradN[i] * AGradN[j]; // Stabilization: (a * Grad(v)) * TauOne * (a * Grad(u))
1748 //
1749 // // q-p stabilization block (reset result)
1750 // L = 0;
1751 //
1752 // for (unsigned int m = 0; m < TDim; ++m) // iterate over v components (vx,vy[,vz])
1753 // {
1754 // // Velocity block
1755 // // K += Viscosity * rShapeDeriv(i, m) * rShapeDeriv(j, m); // Diffusive term: Viscosity * Grad(v) * Grad(u)
1756 // // Note that we are usig kinematic viscosity, as we will multiply it by density later
1757 //
1758 // // v * Grad(p) block
1759 // G = TauOne * AGradN[i] * rShapeDeriv(j, m); // Stabilization: (a * Grad(v)) * TauOne * Grad(p)
1760 // PDivV = rShapeDeriv(i, m) * rShapeFunc[j]; // Div(v) * p
1761 //
1762 // // Write v * Grad(p) component
1763 // rDampingMatrix(FirstRow + m, FirstCol + TDim) += Weight * (G - PDivV);
1764 // // Use symmetry to write the q * rho * Div(u) component
1765 // rDampingMatrix(FirstCol + TDim, FirstRow + m) += Coef * (G + PDivV);
1766 //
1767 // // q-p stabilization block
1768 // L += rShapeDeriv(i, m) * rShapeDeriv(j, m); // Stabilization: Grad(q) * TauOne * Grad(p)
1769 //
1770 // for (unsigned int n = 0; n < TDim; ++n) // iterate over u components (ux,uy[,uz])
1771 // {
1772 // // Velocity block
1773 // rDampingMatrix(FirstRow + m, FirstCol + n) += Coef * (TauTwo + TauCoef) * rShapeDeriv(i, m) * rShapeDeriv(j, n); // Stabilization: Div(v) * TauTwo *( 1+TauOne/Dt) * Div(u)
1774 // rDampRHS[FirstRow + m] -= Coef * TauCoef * rShapeDeriv(i, m) * rShapeDeriv(j, n) * OldVelocity[n]; // Stabilization: Div(v) * TauTwo*TauOne/Dt * Div(u_old)
1775 // }
1776 //
1777 // }
1778 //
1779 // // Write remaining terms to velocity block
1780 // K *= Coef; // Weight by nodal area and density
1781 // for (unsigned int d = 0; d < TDim; ++d)
1782 // rDampingMatrix(FirstRow + d, FirstCol + d) += K;
1783 //
1784 // // Write q-p stabilization block
1785 // rDampingMatrix(FirstRow + TDim, FirstCol + TDim) += Weight * TauOne * L;
1786 //
1787 // // Operate on RHS
1788 // L = 0; // We reuse one of the temporary variables for the pressure RHS
1789 //
1790 // for (unsigned int d = 0; d < TDim; ++d)
1791 // {
1792 // rDampRHS[FirstRow + d] += Coef * TauOne * AGradN[i] * rShapeFunc[j] * rBodyForce[d]; // ( a * Grad(v) ) * TauOne * (Density * BodyForce)
1793 // L += rShapeDeriv(i, d) * rShapeFunc[j] * rBodyForce[d];
1794 // }
1795 // rDampRHS[FirstRow + TDim] += Coef * TauOne * L; // Grad(q) * TauOne * (Density * BodyForce)
1796 //
1797 // // Update reference row index for next iteration
1798 // FirstRow += BlockSize;
1799 // }
1800 //
1801 // // Update reference indices
1802 // FirstRow = 0;
1803 // FirstCol += BlockSize;
1804 // }
1805 //
1806 // // this->AddBTransCB(rDampingMatrix,rShapeDeriv,Viscosity*Coef);
1807 // this->AddViscousTerm(rDampingMatrix,rShapeDeriv,Viscosity*Coef);
1808 // }
1809 
1810  */
1812 
1820  const double Density,
1821  array_1d< double, 3 > & rElementalMomRes,
1822  double& rElementalMassRes,
1823  const array_1d< double, TNumNodes >& rShapeFunc,
1824  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
1825  const double Weight)
1826  {
1827  // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
1829  this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1830 //G
1831  double FluidFraction, FluidFractionRate;
1832  array_1d<double,3> FluidFractionGradient(3,0.0);
1833  this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
1834  this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
1835  this->EvaluateInPoint(FluidFractionRate,FLUID_FRACTION_RATE,rShapeFunc);
1836 //Z
1837  // Compute contribution to Kij * Uj, with Kij = Ni * Residual(Nj); Uj = (v,p)Node_j (column vector)
1838  for (unsigned int i = 0; i < TNumNodes; ++i) // Iterate over element nodes
1839  {
1840 
1841  // Variable references
1842  const array_1d< double, 3 > & rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
1843  const array_1d< double, 3 > & rBodyForce = this->GetGeometry()[i].FastGetSolutionStepValue(BODY_FORCE);
1844  const double& rPressure = this->GetGeometry()[i].FastGetSolutionStepValue(PRESSURE);
1845 
1846  // Compute this node's contribution to the residual (evaluated at inegration point)
1847  for (unsigned int d = 0; d < TDim; ++d)
1848  {
1849 
1850  rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * rBodyForce[d] - AGradN[i] * rVelocity[d]) - rShapeDeriv(i, d) * rPressure);
1851 //G
1852  // rElementalMassRes -= Weight * rShapeDeriv(i, d) * rVelocity[d];
1853  rElementalMassRes -= Weight * (FluidFraction * rShapeDeriv(i, d) * rVelocity[d] + FluidFractionGradient[d] * rShapeFunc[i] * rVelocity[d]);
1854  }
1855 
1856  }
1857 
1858  rElementalMassRes -= Weight * FluidFractionRate;
1859 //Z
1860  }
1861 
1863 
1873  const double Density,
1874  array_1d< double, 3 > & rElementalMomRes,
1875  const array_1d< double, TNumNodes >& rShapeFunc,
1876  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
1877  const double Weight)
1878  {
1879  // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
1881  this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1882  // Compute contribution to Kij * Uj, with Kij = Ni * Residual(Nj); Uj = (v,p)Node_j (column vector)
1883  for (unsigned int i = 0; i < TNumNodes; ++i) // Iterate over element nodes
1884  {
1885 
1886  // Variable references
1887  const array_1d< double, 3 > & rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
1888  const array_1d< double, 3 > & rAcceleration = this->GetGeometry()[i].FastGetSolutionStepValue(ACCELERATION);
1889  const array_1d< double, 3 > & rBodyForce = this->GetGeometry()[i].FastGetSolutionStepValue(BODY_FORCE);
1890  const double& rPressure = this->GetGeometry()[i].FastGetSolutionStepValue(PRESSURE);
1891 
1892  // Compute this node's contribution to the residual (evaluated at inegration point)
1893  for (unsigned int d = 0; d < TDim; ++d)
1894  {
1895  rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * (rBodyForce[d] - rAcceleration[d]) - AGradN[i] * rVelocity[d]) - rShapeDeriv(i, d) * rPressure);
1896  }
1897  }
1898  }
1899 
1901 
1911  const double Density,
1912  array_1d< double, 3 > & rElementalMomRes,
1913  const array_1d< double, TNumNodes >& rShapeFunc,
1914  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
1915  const double Weight)
1916  {
1917  // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
1919  this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
1920  // Compute contribution to Kij * Uj, with Kij = Ni * Residual(Nj); Uj = (v,p)Node_j (column vector)
1921  for (unsigned int i = 0; i < TNumNodes; ++i) // Iterate over element nodes
1922  {
1923 
1924  // Variable references
1925  const array_1d< double, 3 > & rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
1926  const array_1d< double, 3 > & rBodyForce = this->GetGeometry()[i].FastGetSolutionStepValue(BODY_FORCE);
1927  const array_1d< double, 3 > & rProjection = this->GetGeometry()[i].FastGetSolutionStepValue(ADVPROJ);
1928  const double& rPressure = this->GetGeometry()[i].FastGetSolutionStepValue(PRESSURE);
1929 
1930  // Compute this node's contribution to the residual (evaluated at inegration point)
1931  for (unsigned int d = 0; d < TDim; ++d)
1932  {
1933  rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * rBodyForce[d] - AGradN[i] * rVelocity[d]) - rShapeDeriv(i, d) * rPressure);
1934  rElementalMomRes[d] -= Weight * rShapeFunc[i] * rProjection[d];
1935 
1936  }
1937  }
1938  }
1939 
1940 
1942 
1950  virtual void GetEffectiveViscosity(const double Density,
1951  const double MolecularViscosity,
1952  const array_1d<double, TNumNodes>& rShapeFunc,
1953  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
1954  double& TotalViscosity,
1955  const ProcessInfo& rCurrentProcessInfo)
1956  {
1957  const double C = this->GetValue(C_SMAGORINSKY);
1958 
1959  TotalViscosity = MolecularViscosity;
1960  if (C != 0.0 )
1961  {
1962  // The filter width in Smagorinsky is typically the element size h. We will store the square of h, as the final formula involves the squared filter width
1963  const double FilterWidth = this->FilterWidth(rShapeDeriv);
1964 
1965  const double NormS = this->SymmetricGradientNorm(rShapeDeriv);
1966 
1967  // Total Viscosity
1968  TotalViscosity += 2.0 * C * C * FilterWidth * NormS;
1969  }
1970  }
1971 
1973 
1979  virtual void GetAdvectiveVel(array_1d< double, 3 > & rAdvVel,
1980  const array_1d< double, TNumNodes >& rShapeFunc)
1981  {
1982  // Compute the weighted value of the advective velocity in the (Gauss) Point
1983  rAdvVel = rShapeFunc[0] * (this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY));
1984  for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
1985  rAdvVel += rShapeFunc[iNode] * (this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY));
1986  }
1987 
1989 
1996  virtual void GetAdvectiveVel(array_1d< double, 3 > & rAdvVel,
1997  const array_1d< double, TNumNodes >& rShapeFunc,
1998  const std::size_t Step)
1999  {
2000  // Compute the weighted value of the advective velocity in the (Gauss) Point
2001  rAdvVel = rShapeFunc[0] * (this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY, Step));
2002  for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
2003  rAdvVel += rShapeFunc[iNode] * (this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY, Step));
2004  }
2005 
2006 //G
2008 
2015  virtual void GetAdvectiveVelDivergence(double & rAdvVelDiv,
2016  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
2017  {
2018  rAdvVelDiv = 0.0;
2019 
2020  for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode){ // loop over nodes
2021  const array_1d< double, 3 > vel_at_nodes = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY);
2022 
2023  for (unsigned int d = 1; d < TDim; ++d){
2024  // loop over components
2025  rAdvVelDiv += vel_at_nodes[d] * rShapeDeriv(iNode, d);
2026  }
2027 
2028  }
2029 
2030  }
2031 
2032  virtual void GetAdvectiveVelDivergence(double & rAdvVelDiv,
2033  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
2034  const std::size_t Step)
2035  {
2036  rAdvVelDiv = 0.0;
2037 
2038  for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode){ // loop over nodes
2039  array_1d< double, 3 > vel_at_nodes = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY, Step);
2040 
2041  for (unsigned int d = 1; d < TDim; ++d){
2042  // loop over components
2043  rAdvVelDiv += vel_at_nodes[d] * rShapeDeriv(iNode, d);
2044  }
2045 
2046  }
2047 
2048  }
2049 //Z
2050 
2052 
2060  const array_1d< double, 3 > & rVelocity,
2061  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
2062  {
2063  // Evaluate (and weight) the a * Grad(Ni) operator in the integration point, for each node i
2064  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode) // Loop over nodes
2065  {
2066  // Initialize result
2067  rResult[iNode] = rVelocity[0] * rShapeDeriv(iNode, 0);
2068  for (unsigned int d = 1; d < TDim; ++d) // loop over components
2069  rResult[iNode] += rVelocity[d] * rShapeDeriv(iNode, d);
2070  }
2071  }
2072 
2074 
2087 
2095  const array_1d< double, 3 > & rVelocity,
2096  const double & rVelocityDiv,
2097  const array_1d< double, TNumNodes >& rShapeFunc,
2098  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
2099  {
2100  // Evaluate (and weight) the a * Grad(Ni) + div(a) * Ni operator in the integration point, for each node i
2101  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode){ // Loop over nodes{
2102  // Initialize result
2103  rResult[iNode] = rVelocityDiv * rShapeFunc[iNode];
2104 
2105  for (unsigned int d = 0; d < TDim; ++d){ // loop over components
2106  rResult[iNode] += rVelocity[d] * rShapeDeriv(iNode, d);
2107  }
2108 
2109  }
2110  }
2111 
2112  virtual void AddPointContribution(double& rResult,
2113  const Variable< double >& rVariable,
2114  const array_1d< double, TNumNodes >& rShapeFunc,
2115  const double Weight = 1.0)
2116  {
2117  // Compute the weighted value of the nodal variable in the (Gauss) Point
2118  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
2119  rResult += Weight * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2120  }
2121 
2123 
2132  virtual void EvaluateInPoint(double& rResult,
2133  const Variable< double >& rVariable,
2134  const array_1d< double, TNumNodes >& rShapeFunc)
2135  {
2136  // Compute the weighted value of the nodal variable in the (Gauss) Point
2137  rResult = rShapeFunc[0] * this->GetGeometry()[0].FastGetSolutionStepValue(rVariable);
2138 
2139  for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
2140  rResult += rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2141  }
2142 
2143 //G
2144 
2146  const Variable< double >& rVariable,
2147  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
2148  {
2149 
2150  for (unsigned int i = 0; i < TNumNodes; ++i) {
2151  double& scalar = this->GetGeometry()[i].FastGetSolutionStepValue(rVariable);
2152 
2153  for (unsigned int d = 0; d < TDim; ++d){
2154  rResult[d] += rShapeDeriv(i, d) * scalar;
2155  }
2156  }
2157  }
2158 
2160  const Variable< array_1d<double, 3 > >& rVariable,
2161  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
2162  {
2163  for (unsigned int d1 = 0; d1 < TDim; ++d1) {
2164 
2165  for (unsigned int i = 0; i < TNumNodes; ++i) {
2166  array_1d<double, 3 >& vector = this->GetGeometry()[i].FastGetSolutionStepValue(rVariable);
2167 
2168  for (unsigned int d2 = 0; d2 < TDim; ++d2) {
2169  rResult(d1, d2) += rShapeDeriv(i, d2) * vector[d1];
2170  }
2171  }
2172  }
2173  }
2174 
2175  virtual void EvaluateTimeDerivativeInPoint(double& rResult,
2176  const Variable< double >& rVariable,
2177  const array_1d< double, TNumNodes >& rShapeFunc,
2178  const double& DeltaTime,
2179  const std::vector<double>& rSchemeWeigths)
2180  {
2181  // Compute the time derivative of a nodal variable as a liner contribution of weighted value of the nodal variable in the (Gauss) Point
2182 
2183  if (rVariable == FLUID_FRACTION_RATE){
2184 // int index = 0;
2185  double delta_time_inv = 1.0 / DeltaTime;
2186 
2187 // for (unsigned int iWeight = 0; iWeight < rSchemeWeigths.size(); ++iWeight){
2188 
2189 // for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
2190 
2191 // rResult += rSchemeWeigths[iWeight] * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION, index);
2192 // }
2193 
2194 // ++index;
2195 
2196 // }
2197 
2198 
2199  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
2200  double rate = delta_time_inv * (this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION) - this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION_OLD));
2201  this->GetGeometry()[iNode].SetLock();
2202  this->GetGeometry()[iNode].FastGetSolutionStepValue(FLUID_FRACTION_RATE) = rate;
2203  this->GetGeometry()[iNode].UnSetLock();
2204  rResult += rShapeFunc[iNode] * rate;
2205  }
2206 
2207  }
2208 
2209  }
2210 //Z
2211 
2213 
2224  const Variable< array_1d< double, 3 > >& rVariable,
2225  const array_1d< double, TNumNodes>& rShapeFunc,
2226  const double Weight = 1.0)
2227  {
2228  // Compute the weighted value of the nodal variable in the (Gauss) Point
2229  for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
2230  rResult += Weight * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2231  }
2232 
2234 
2242  virtual void EvaluateInPoint(array_1d< double, 3 > & rResult,
2243  const Variable< array_1d< double, 3 > >& rVariable,
2244  const array_1d< double, TNumNodes >& rShapeFunc)
2245  {
2246  // Compute the weighted value of the nodal variable in the (Gauss) Point
2247  rResult = rShapeFunc[0] * this->GetGeometry()[0].FastGetSolutionStepValue(rVariable);
2248  for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
2249  rResult += rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
2250  }
2251 
2253 
2260  double ElementSize(const double);
2261 
2263  double FilterWidth();
2264 
2267 
2269 
2274  {
2275  const unsigned int GradientSize = (TDim*(TDim+1))/2; // Number of different terms in the symmetric gradient matrix
2276  array_1d<double,GradientSize> GradientVector( GradientSize, 0.0 );
2277  unsigned int Index;
2278 
2279  // Compute Symmetric Grad(u). Note that only the lower half of the matrix is calculated
2280  for (unsigned int k = 0; k < TNumNodes; ++k)
2281  {
2282  const array_1d< double, 3 > & rNodeVel = this->GetGeometry()[k].FastGetSolutionStepValue(VELOCITY);
2283  Index = 0;
2284  for (unsigned int i = 0; i < TDim; ++i)
2285  {
2286  for (unsigned int j = 0; j < i; ++j) // Off-diagonal
2287  GradientVector[Index++] += 0.5 * (rShapeDeriv(k, j) * rNodeVel[i] + rShapeDeriv(k, i) * rNodeVel[j]);
2288  GradientVector[Index++] += rShapeDeriv(k, i) * rNodeVel[i]; // Diagonal
2289  }
2290  }
2291 
2292  // Norm[ Symmetric Grad(u) ] = ( 2 * Sij * Sij )^(1/2)
2293  Index = 0;
2294  double NormS(0.0);
2295  for (unsigned int i = 0; i < TDim; ++i)
2296  {
2297  for (unsigned int j = 0; j < i; ++j)
2298  {
2299  NormS += 2.0 * GradientVector[Index] * GradientVector[Index]; // Using symmetry, lower half terms of the matrix are added twice
2300  ++Index;
2301  }
2302  NormS += GradientVector[Index] * GradientVector[Index]; // Diagonal terms
2303  ++Index; // Diagonal terms
2304  }
2305 
2306  NormS = sqrt( 2.0 * NormS );
2307  return NormS;
2308  }
2309 
2311 
2317  virtual void AddViscousTerm(MatrixType& rDampingMatrix,
2318  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
2319  const double Weight);
2320 
2322 
2332  void AddBTransCB(MatrixType& rDampingMatrix,
2333  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
2334  const double Weight)
2335  {
2336  BoundedMatrix<double, (TDim * TNumNodes)/2, TDim*TNumNodes > B;
2337  BoundedMatrix<double, (TDim * TNumNodes)/2, (TDim*TNumNodes)/2 > C;
2338  this->CalculateB(B,rShapeDeriv);
2339  this->CalculateC(C,Weight);
2340 
2341  const unsigned int BlockSize = TDim + 1;
2342  const unsigned int StrainSize = (TDim*TNumNodes)/2;
2343 
2344  DenseVector<unsigned int> aux(TDim*TNumNodes);
2345  for(unsigned int i=0; i<TNumNodes; i++)
2346  {
2347  int base_index = TDim*i;
2348  int aux_index = BlockSize*i;
2349  for(unsigned int j=0; j<TDim; j++)
2350  {
2351  aux[base_index+j] = aux_index+j;
2352  }
2353  }
2354 
2355  for(unsigned int k=0; k< StrainSize; k++)
2356  {
2357  for(unsigned int l=0; l< StrainSize; l++)
2358  {
2359  const double Ckl = C(k,l);
2360  for (unsigned int i = 0; i < TDim*TNumNodes; ++i) // iterate over v components (vx,vy[,vz])
2361  {
2362  const double Bki=B(k,i);
2363  for (unsigned int j = 0; j < TDim*TNumNodes; ++j) // iterate over u components (ux,uy[,uz])
2364  {
2365  rDampingMatrix(aux[i],aux[j]) += Bki*Ckl*B(l,j);
2366  }
2367 
2368  }
2369 
2370  }
2371  }
2372  }
2373 
2376  const double Weight)
2377  {
2378  const GeometryType& rGeom = this->GetGeometry();
2379 
2380  // Velocity gradient
2381  MatrixType GradU = ZeroMatrix(TDim,TDim);
2382  for (unsigned int n = 0; n < TNumNodes; n++)
2383  {
2384  const array_1d<double,3>& rVel = this->GetGeometry()[n].FastGetSolutionStepValue(VELOCITY);
2385  for (unsigned int i = 0; i < TDim; i++)
2386  for (unsigned int j = 0; j < TDim; j++)
2387  GradU(i,j) += rDN_DX(n,j)*rVel[i];
2388  }
2389 
2390  // Element lengths
2391  array_1d<double,3> Delta(3,0.0);
2392  Delta[0] = fabs(rGeom[TNumNodes-1].X()-rGeom[0].X());
2393  Delta[1] = fabs(rGeom[TNumNodes-1].Y()-rGeom[0].Y());
2394  Delta[2] = fabs(rGeom[TNumNodes-1].Z()-rGeom[0].Z());
2395 
2396  for (unsigned int n = 1; n < TNumNodes; n++)
2397  {
2398  double hx = fabs(rGeom[n].X()-rGeom[n-1].X());
2399  if (hx > Delta[0]) Delta[0] = hx;
2400  double hy = fabs(rGeom[n].Y()-rGeom[n-1].Y());
2401  if (hy > Delta[1]) Delta[1] = hy;
2402  double hz = fabs(rGeom[n].Z()-rGeom[n-1].Z());
2403  if (hz > Delta[2]) Delta[2] = hz;
2404  }
2405 
2406  double AvgDeltaSq = Delta[0];
2407  for (unsigned int d = 1; d < TDim; d++)
2408  AvgDeltaSq *= Delta[d];
2409  AvgDeltaSq = std::pow(AvgDeltaSq,2./TDim);
2410 
2411  Delta[0] = Delta[0]*Delta[0]/12.0;
2412  Delta[1] = Delta[1]*Delta[1]/12.0;
2413  Delta[2] = Delta[2]*Delta[2]/12.0;
2414 
2415  // Gij
2416  MatrixType G = ZeroMatrix(TDim,TDim);
2417  for (unsigned int i = 0; i < TDim; i++)
2418  for (unsigned int j = 0; j < TDim; j++)
2419  for (unsigned int d = 0; d < TDim; d++)
2420  G(i,j) += Delta[d]*GradU(i,d)*GradU(j,d);
2421 
2422  // Gij:Sij
2423  double GijSij = 0.0;
2424  for (unsigned int i = 0; i < TDim; i++)
2425  for (unsigned int j = 0; j < TDim; j++)
2426  GijSij += 0.5*G(i,j)*( GradU(i,j) + GradU(j,i) );
2427 
2428  if (GijSij < 0.0) // Otherwise model term is clipped
2429  {
2430  // Gkk
2431  double Gkk = G(0,0);
2432  for (unsigned int d = 1; d < TDim; d++)
2433  Gkk += G(d,d);
2434 
2435  // C_epsilon
2436  const double Ce = 1.0;
2437 
2438  // ksgs
2439  double ksgs = -4*AvgDeltaSq*GijSij/(Ce*Ce*Gkk);
2440 
2441  // Assembly of model term
2442  unsigned int RowIndex = 0;
2443  unsigned int ColIndex = 0;
2444 
2445  for (unsigned int i = 0; i < TNumNodes; i++)
2446  {
2447  for (unsigned int j = 0; j < TNumNodes; j++)
2448  {
2449  for (unsigned int d = 0; d < TDim; d++)
2450  {
2451  double Aux = rDN_DX(i,d) * Delta[0] * G(d,0)*rDN_DX(j,0);
2452  for (unsigned int k = 1; k < TDim; k++)
2453  Aux += rDN_DX(i,d) *Delta[k] * G(d,k)*rDN_DX(j,k);
2454  rDampingMatrix(RowIndex+d,ColIndex+d) += Weight * 2.0*ksgs * Aux;
2455  }
2456 
2457  ColIndex += TDim;
2458  }
2459  RowIndex += TDim;
2460  ColIndex = 0;
2461  }
2462  }
2463 
2464  }
2465 
2467 
2472  void CalculateB( BoundedMatrix<double, (TDim * TNumNodes) / 2, TDim * TNumNodes >& rB,
2473  const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv);
2474 
2476 
2483  virtual void CalculateC( BoundedMatrix<double, (TDim * TNumNodes) / 2, (TDim * TNumNodes) / 2 >& rC,
2484  const double Viscosity);
2485 
2486  double ConsistentMassCoef(const double Area);
2487 
2491 
2492 
2496 
2497 
2501 
2502 
2504 
2505 private:
2508 
2509 
2513 
2517 
2518  friend class Serializer;
2519 
2520  virtual void save(Serializer& rSerializer) const override
2521  {
2523  }
2524 
2525  virtual void load(Serializer& rSerializer) override
2526  {
2528  }
2529 
2533 
2534 
2538 
2539 
2543 
2544 
2548 
2549 
2553 
2555  MonolithicDEMCoupled & operator=(MonolithicDEMCoupled const& rOther);
2556 
2558  MonolithicDEMCoupled(MonolithicDEMCoupled const& rOther);
2559 
2561 
2562 }; // Class MonolithicDEMCoupled
2563 
2565 
2568 
2569 
2573 
2574 
2576 template< unsigned int TDim,
2577  unsigned int TNumNodes >
2578 inline std::istream& operator >>(std::istream& rIStream,
2580 {
2581  return rIStream;
2582 }
2583 
2585 template< unsigned int TDim,
2586  unsigned int TNumNodes >
2587 inline std::ostream& operator <<(std::ostream& rOStream,
2589 {
2590  rThis.PrintInfo(rOStream);
2591  rOStream << std::endl;
2592  rThis.PrintData(rOStream);
2593 
2594  return rOStream;
2595 }
2597 
2599 
2600 } // namespace Kratos.
2601 
2602 #endif // MONOLITHIC_DEM_COUPLED_H
PeriodicInterfaceProcess & operator=(const PeriodicInterfaceProcess &)=delete
Base class for all Elements.
Definition: element.h:60
std::vector< DofType::Pointer > DofsVectorType
Definition: element.h:100
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: element.h:1135
virtual int Check(const ProcessInfo &rCurrentProcessInfo) const
Definition: element.h:904
std::vector< std::size_t > EquationIdVectorType
Definition: element.h:98
std::size_t IndexType
Definition: flags.h:74
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometrical_object.h:248
Geometry base class.
Definition: geometry.h:71
static void CalculateGeometryData(const GeometryType &rGeometry, BoundedMatrix< double, 4, 3 > &rDN_DX, array_1d< double, 4 > &rN, double &rVolume)
This function is designed to compute the shape function derivatives, shape functions and volume in 3D...
Definition: geometry_utilities.h:176
This object defines an indexed object.
Definition: indexed_object.h:54
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
A stabilized element for the incompressible Navier-Stokes equations.
Definition: monolithic_dem_coupled.h:137
void AddBTransCB(MatrixType &rDampingMatrix, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Adds the contribution of the viscous term to the momentum equation (alternate).
Definition: monolithic_dem_coupled.h:2332
virtual void AddPointContribution(double &rResult, const Variable< double > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight=1.0)
Definition: monolithic_dem_coupled.h:2112
virtual void EvaluateGradientOfVectorInPoint(BoundedMatrix< double, TDim, TDim > &rResult, const Variable< array_1d< double, 3 > > &rVariable, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Definition: monolithic_dem_coupled.h:2159
void GetConvectionOperator(array_1d< double, TNumNodes > &rResult, const array_1d< double, 3 > &rVelocity, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Write the convective operator evaluated at this point (for each nodal funciton) to an array.
Definition: monolithic_dem_coupled.h:2059
IndexedObject BaseType
base type: an IndexedObject that automatically has a unique number
Definition: monolithic_dem_coupled.h:146
virtual void CalculateLeftHandSide(MatrixType &rLeftHandSideMatrix, const ProcessInfo &rCurrentProcessInfo) override
Returns a zero matrix of appropiate size (provided for compatibility with scheme)
Definition: monolithic_dem_coupled.h:304
virtual void GetEffectiveViscosity(const double Density, const double MolecularViscosity, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, double &TotalViscosity, const ProcessInfo &rCurrentProcessInfo)
Add the contribution from Smagorinsky model to the (kinematic) viscosity.
Definition: monolithic_dem_coupled.h:1950
void CalculateB(BoundedMatrix< double,(TDim *TNumNodes)/2, TDim *TNumNodes > &rB, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Calculate the strain rate matrix.
virtual void CalculateMassMatrix(MatrixType &rMassMatrix, const ProcessInfo &rCurrentProcessInfo) override
Computes local contributions to the mass matrix.
Definition: monolithic_dem_coupled.h:467
virtual void FinalizeNonLinearIteration(const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:662
virtual void CalculateOnIntegrationPoints(const Variable< array_1d< double, 6 > > &rVariable, std::vector< array_1d< double, 6 > > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:1116
virtual void CalculateLocalVelocityContribution(MatrixType &rDampingMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Computes the local contribution associated to 'new' velocity and pressure values.
Definition: monolithic_dem_coupled.h:576
virtual ~MonolithicDEMCoupled()
Destructor.
Definition: monolithic_dem_coupled.h:230
void FinalizeSolutionStep(const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:707
Geometry< NodeType >::PointsArrayType NodesArrayType
definition of nodes container type, redefined from GeometryType
Definition: monolithic_dem_coupled.h:161
virtual void AddRHSLaplacian(VectorType &F, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Definition: monolithic_dem_coupled.h:1336
MonolithicDEMCoupled(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
Constuctor using geometry and properties.
Definition: monolithic_dem_coupled.h:225
virtual void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const override
Provides the global indices for each one of this element's local rows.
virtual void AddViscousTerm(MatrixType &rDampingMatrix, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Adds the contribution of the viscous term to the momentum equation.
virtual void AddMomentumRHS(VectorType &F, const double Density, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight)
Add the momentum equation contribution to the RHS (body forces)
Definition: monolithic_dem_coupled.h:1313
virtual void EvaluateInPoint(double &rResult, const Variable< double > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc)
Write the value of a variable at a point inside the element to a double.
Definition: monolithic_dem_coupled.h:2132
virtual void AddMassRHS(VectorType &F, const double Density, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight, const std::vector< double > &TimeSchemeWeights, const double &DeltaTime)
Definition: monolithic_dem_coupled.h:1357
MonolithicDEMCoupled(IndexType NewId, GeometryType::Pointer pGeometry)
Constructor using a geometry object.
Definition: monolithic_dem_coupled.h:215
virtual void EvaluateGradientOfScalarInPoint(array_1d< double, 3 > &rResult, const Variable< double > &rVariable, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Definition: monolithic_dem_coupled.h:2145
std::vector< Dof< double >::Pointer > DofsVectorType
Definition: monolithic_dem_coupled.h:173
void CalculateWeights(ShapeFunctionDerivativesArrayType &rDN_DX, Matrix &rNContainer, Vector &rGaussWeights)
Determine integration point weights and shape funcition derivatives from the element's geometry.
virtual void AddProjectionToRHS(VectorType &RHS, const array_1d< double, 3 > &rAdvVel, const double Density, const double TauOne, const double TauTwo, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight, const double DeltaTime=1.0)
Add OSS projection terms to the RHS.
Definition: monolithic_dem_coupled.h:1379
virtual void Calculate(const Variable< double > &rVariable, double &rOutput, const ProcessInfo &rCurrentProcessInfo) override
Implementation of Calculate to compute an error estimate.
Definition: monolithic_dem_coupled.h:729
void AddMassStabTerms(MatrixType &rLHSMatrix, const double Density, const array_1d< double, 3 > &rAdvVel, const double TauOne, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Add mass-like stabilization terms to LHS.
Definition: monolithic_dem_coupled.h:1510
void AddConsistentMassMatrixContribution(MatrixType &rLHSMatrix, const array_1d< double, TNumNodes > &rShapeFunc, const double Density, const double Weight)
Definition: monolithic_dem_coupled.h:1463
double SymmetricGradientNorm(const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Compute the norm of the symmetric gradient of velocity.
Definition: monolithic_dem_coupled.h:2273
virtual void CalculateOnIntegrationPoints(const Variable< array_1d< double, 3 > > &rVariable, std::vector< array_1d< double, 3 > > &rOutput, const ProcessInfo &rCurrentProcessInfo) override
Obtain an array_1d<double,3> elemental variable, evaluated on gauss points.
void GetModifiedConvectionOperator(array_1d< double, TNumNodes > &rResult, const array_1d< double, 3 > &rVelocity, const double &rVelocityDiv, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Add the weighted value of a variable at a point inside the element to a double.
Definition: monolithic_dem_coupled.h:2094
virtual void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: monolithic_dem_coupled.h:1205
virtual void GetSecondDerivativesVector(Vector &Values, int Step=0) const override
Returns ACCELERATION_X, ACCELERATION_Y, (ACCELERATION_Z,) 0 for each node.
virtual void GetAdvectiveVel(array_1d< double, 3 > &rAdvVel, const array_1d< double, TNumNodes > &rShapeFunc)
Write the advective velocity evaluated at this point to an array.
Definition: monolithic_dem_coupled.h:1979
Geometry< NodeType > GeometryType
definition of the geometry type with given NodeType
Definition: monolithic_dem_coupled.h:158
virtual void AddPointContribution(array_1d< double, 3 > &rResult, const Variable< array_1d< double, 3 > > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc, const double Weight=1.0)
Add the weighted value of a variable at a point inside the element to a vector.
Definition: monolithic_dem_coupled.h:2223
void ModulatedGradientDiffusion(MatrixType &rDampingMatrix, const BoundedMatrix< double, TNumNodes, TDim > &rDN_DX, const double Weight)
Definition: monolithic_dem_coupled.h:2374
Vector VectorType
Definition: monolithic_dem_coupled.h:163
virtual std::string Info() const override
Turn back information as a string.
Definition: monolithic_dem_coupled.h:1197
std::vector< std::size_t > EquationIdVectorType
Definition: monolithic_dem_coupled.h:171
GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType
Type for an array of shape function gradient matrices.
Definition: monolithic_dem_coupled.h:185
void ASGSMomResidual(const array_1d< double, 3 > &rAdvVel, const double Density, array_1d< double, 3 > &rElementalMomRes, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Assemble the contribution from an integration point to the element's residual.
Definition: monolithic_dem_coupled.h:1872
void OSSMomResidual(const array_1d< double, 3 > &rAdvVel, const double Density, array_1d< double, 3 > &rElementalMomRes, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Assemble the contribution from an integration point to the element's residual.
Definition: monolithic_dem_coupled.h:1910
virtual void GetFirstDerivativesVector(Vector &Values, int Step=0) const override
Returns VELOCITY_X, VELOCITY_Y, (VELOCITY_Z,) PRESSURE for each node.
Element::Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const override
Create a new element of this type.
Definition: monolithic_dem_coupled.h:251
virtual void GetAdvectiveVelDivergence(double &rAdvVelDiv, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const std::size_t Step)
Definition: monolithic_dem_coupled.h:2032
double FilterWidth()
Return the filter width for the smagorinsky model (Delta squared)
virtual void GetDofList(DofsVectorType &rElementalDofList, const ProcessInfo &rCurrentProcessInfo) const override
Returns a list of the element's Dofs.
virtual void CalculateTau(double &TauOne, double &TauTwo, const array_1d< double, 3 > &rAdvVel, const double Area, const double Density, const double KinViscosity, const ProcessInfo &rCurrentProcessInfo)
Calculate Stabilization parameters.
Definition: monolithic_dem_coupled.h:1256
Matrix MatrixType
Definition: monolithic_dem_coupled.h:165
virtual void CalculateRightHandSide(VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Provides local contributions from body forces and projections to the RHS.
Definition: monolithic_dem_coupled.h:335
MonolithicDEMCoupled(IndexType NewId, const NodesArrayType &ThisNodes)
Constructor using an array of nodes.
Definition: monolithic_dem_coupled.h:206
virtual void Calculate(const Variable< array_1d< double, 3 > > &rVariable, array_1d< double, 3 > &rOutput, const ProcessInfo &rCurrentProcessInfo) override
Implementation of Calculate to compute the local OSS projections.
Definition: monolithic_dem_coupled.h:835
KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(MonolithicDEMCoupled)
Pointer definition of MonolithicDEMCoupled.
double ConsistentMassCoef(const double Area)
std::size_t SizeType
Definition: monolithic_dem_coupled.h:169
virtual int Check(const ProcessInfo &rCurrentProcessInfo) const override
Checks the input and that all required Kratos variables have been registered.
Definition: monolithic_dem_coupled.h:1148
PointerVectorSet< Dof< double >, IndexedObject > DofsArrayType
Definition: monolithic_dem_coupled.h:175
void AddProjectionResidualContribution(const array_1d< double, 3 > &rAdvVel, const double Density, array_1d< double, 3 > &rElementalMomRes, double &rElementalMassRes, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Assemble the contribution from an integration point to the element's residual.
Definition: monolithic_dem_coupled.h:1819
MonolithicDEMCoupled(IndexType NewId=0)
Default constuctor.
Definition: monolithic_dem_coupled.h:197
void CalculateLumpedMassMatrix(MatrixType &rLHSMatrix, const double Mass)
Add lumped mass matrix.
Definition: monolithic_dem_coupled.h:1432
void CalculateLaplacianLumpedMassMatrix(MatrixType &rLHSMatrix, const double Mass)
Definition: monolithic_dem_coupled.h:1448
virtual void EvaluateInPoint(array_1d< double, 3 > &rResult, const Variable< array_1d< double, 3 > > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc)
Write the value of a variable at a point inside the element to a double.
Definition: monolithic_dem_coupled.h:2242
Properties PropertiesType
Definition: monolithic_dem_coupled.h:155
Kratos::Matrix ShapeFunctionDerivativesType
Type for a matrix containing the shape function gradients.
Definition: monolithic_dem_coupled.h:182
virtual void CalculateOnIntegrationPoints(const Variable< Matrix > &rVariable, std::vector< Matrix > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:1126
virtual void CalculateOnIntegrationPoints(const Variable< double > &rVariable, std::vector< double > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Obtain a double elemental variable, evaluated on gauss points.
Definition: monolithic_dem_coupled.h:1001
Node NodeType
definition of node type (default is: Node)
Definition: monolithic_dem_coupled.h:149
virtual void CalculateOnIntegrationPoints(const Variable< Vector > &rVariable, std::vector< Vector > &rValues, const ProcessInfo &rCurrentProcessInfo) override
Definition: monolithic_dem_coupled.h:1121
virtual void GetAdvectiveVel(array_1d< double, 3 > &rAdvVel, const array_1d< double, TNumNodes > &rShapeFunc, const std::size_t Step)
Write the advective velocity evaluated at this point to an array.
Definition: monolithic_dem_coupled.h:1996
std::size_t IndexType
Definition: monolithic_dem_coupled.h:167
virtual void EvaluateTimeDerivativeInPoint(double &rResult, const Variable< double > &rVariable, const array_1d< double, TNumNodes > &rShapeFunc, const double &DeltaTime, const std::vector< double > &rSchemeWeigths)
Definition: monolithic_dem_coupled.h:2175
double FilterWidth(const BoundedMatrix< double, TNumNodes, TDim > &DN_DX)
Return the filter width for the smagorinsky model (Delta squared)
Kratos::Vector ShapeFunctionsType
Type for shape function values container.
Definition: monolithic_dem_coupled.h:179
virtual void CalculateC(BoundedMatrix< double,(TDim *TNumNodes)/2,(TDim *TNumNodes)/2 > &rC, const double Viscosity)
Calculate a matrix that provides the stress given the strain rate.
virtual void CalculateStaticTau(double &TauOne, const array_1d< double, 3 > &rAdvVel, const double Area, const double Density, const double KinViscosity)
Calculate momentum stabilization parameter (without time term).
Definition: monolithic_dem_coupled.h:1294
virtual void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Provides local contributions from body forces and OSS projection terms.
Definition: monolithic_dem_coupled.h:269
double ElementSize(const double)
Return an estimate for the element size h, used to calculate the stabilization parameters.
virtual void CalculateLaplacianMassMatrix(MatrixType &rMassMatrix, const ProcessInfo &rCurrentProcessInfo)
Definition: monolithic_dem_coupled.h:542
void AddIntegrationPointVelocityContribution(MatrixType &rDampingMatrix, VectorType &rDampRHS, const double Density, const double Viscosity, const array_1d< double, 3 > &rAdvVel, const double TauOne, const double TauTwo, const array_1d< double, TNumNodes > &rShapeFunc, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv, const double Weight)
Add a the contribution from a single integration point to the velocity contribution.
Definition: monolithic_dem_coupled.h:1566
virtual void GetAdvectiveVelDivergence(double &rAdvVelDiv, const BoundedMatrix< double, TNumNodes, TDim > &rShapeDeriv)
Write the divergence of the advective velocity evaluated at this point to an array.
Definition: monolithic_dem_coupled.h:2015
This class defines the node.
Definition: node.h:65
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
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
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Properties encapsulates data shared by different Elements or Conditions. It can store any type of dat...
Definition: properties.h:69
The serialization consists in storing the state of an object into a storage format like data file or ...
Definition: serializer.h:123
#define KRATOS_SERIALIZE_SAVE_BASE_CLASS(Serializer, BaseType)
Definition: define.h:812
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_TRY
Definition: define.h:109
#define KRATOS_SERIALIZE_LOAD_BASE_CLASS(Serializer, BaseType)
Definition: define.h:815
#define KRATOS_CHECK_DOF_IN_NODE(TheVariable, TheNode)
Definition: checks.h:176
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
#define KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(TheVariable, TheNode)
Definition: checks.h:171
Modeler::Pointer Create(const std::string &ModelerName, Model &rModel, const Parameters ModelParameters)
Checks if the modeler is registered.
Definition: modeler_factory.cpp:30
Parameters GetValue(Parameters &rParameters, const std::string &rEntry)
Definition: add_kratos_parameters_to_python.cpp:53
void CalculateB(const GeometricalObject &rElement, const TMatrixType1 &rDN_DX, TMatrixType2 &rB)
This method computes the deformation tensor B (for small deformation solid elements)
Definition: structural_mechanics_element_utilities.h:105
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
TABLE_NUMBER_ANGULAR_VELOCITY TABLE_NUMBER_MOMENT I33 BEAM_INERTIA_ROT_UNIT_LENGHT_Y KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BEAM_INERTIA_ROT_UNIT_LENGHT_Z) typedef std double
Definition: DEM_application_variables.h:182
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
def SetValue(entity, variable, value)
Definition: coupling_interface_data.py:256
X2
Definition: generate_frictional_mortar_condition.py:120
X1
Definition: generate_frictional_mortar_condition.py:119
int C
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:27
def Index()
Definition: hdf5_io_tools.py:38
F
Definition: hinsberg_optimization.py:144
def load(f)
Definition: ode_solve.py:307
int L
Definition: ode_solve.py:390
int d
Definition: ode_solve.py:397
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
int m
Definition: run_marine_rain_substepping.py:8
J
Definition: sensitivityMatrix.py:58
N
Definition: sensitivityMatrix.py:29
K
Definition: sensitivityMatrix.py:73
B
Definition: sensitivityMatrix.py:76
integer i
Definition: TensorModule.f:17
subroutine d1(DSTRAN, D, dtime, NDI, NSHR, NTENS)
Definition: TensorModule.f:594
integer l
Definition: TensorModule.f:17