12 #if !defined(KRATOS_PURE_CONVECTION_UTILITIES_INCLUDED )
13 #define KRATOS_PURE_CONVECTION_UTILITIES_INCLUDED
29 #include "utilities/geometry_utilities.h"
35 template<
int TDim,
class TSparseSpace,
class TLinearSolver>
65 if( (it->GetValue(NEIGHBOUR_NODES)).size() != 0 )
67 mDofSet.
push_back( it->pGetDof(rScalarVar) );
69 tot_nnz += (it->GetValue(NEIGHBOUR_NODES)).size()+1;
76 int fix_id = mDofSet.
size();
79 if (dof_iterator->IsFixed())
80 dof_iterator->SetEquationId(--fix_id);
82 dof_iterator->SetEquationId(free_id++);
84 mEquationSystemSize = fix_id;
86 std::vector< int > work_array;
87 work_array.reserve(1000);
90 mA.resize(mEquationSystemSize, mEquationSystemSize, tot_nnz);
99 unsigned int index_i = it->GetDof(rScalarVar).EquationId();
101 if(index_i < mEquationSystemSize && (it->GetValue(NEIGHBOUR_NODES)).size() != 0)
109 unsigned int index_j =
i->GetDof(rScalarVar).EquationId();
110 if(index_j < mEquationSystemSize)
111 work_array.push_back(index_j);
115 std::sort(work_array.begin(),work_array.end());
116 typename std::vector<int>::iterator new_end = std::unique(work_array.begin(),work_array.end());
117 unsigned int number_of_entries = new_end - work_array.begin();
120 for(
unsigned int j=0;
j<number_of_entries;
j++)
122 mA.push_back(index_i,work_array[
j] , 0.00);
126 work_array.erase(work_array.begin(),work_array.end());
135 mDx.resize(mA.size1(),
false);
136 mb.resize(mA.size1(),
false);
143 typename TLinearSolver::Pointer linear_solver,
152 TSparseSpace::SetToZero(mA);
153 TSparseSpace::SetToZero(mb);
154 TSparseSpace::SetToZero(mDx);
157 double dt = CurrentProcessInfo[DELTA_TIME];
166 BDFcoeffs[0] = 1.5 /
dt;
167 BDFcoeffs[1] = -2.0 /
dt;
168 BDFcoeffs[2] = 0.5 /
dt;
169 KRATOS_WATCH(
"CONVECTION: BDF_2.......................................................................................");
173 BDFcoeffs[0] = 1.0 /
dt;
174 BDFcoeffs[1] = -1.0 /
dt;
175 KRATOS_WATCH(
"CONVECTION: BDF_1.......................................................................................");
195 for(
unsigned int i=0;
i<TDim+1;
i++)
196 MassFactors[
i] = 1.0/
double(TDim+1);
201 double lumping_factor = 1.0/
double(TDim+1.0);
202 const unsigned int number_of_points = TDim+1;
205 for(ModelPart::ElementsContainerType::iterator
i =
model_part.ElementsBegin();
216 for(
int ii = 0; ii<TDim+1; ii++)
218 local_indices[ii] = geom[ii].GetDof(rScalarVar).EquationId();
227 for(
unsigned int i = 0;
i<number_of_points;
i++)
231 proj += geom[
i].FastGetSolutionStepValue(rProjVar);
232 for(
unsigned int j = 0;
j<TDim;
j++)
234 vel_gauss[
j] +=
v[
j] -
w[
j];
237 vel_gauss *= lumping_factor;
238 proj *= lumping_factor;
246 double tau = CalculateTAU(vel_gauss,Volume);
250 for(
unsigned int iii = 0; iii<number_of_points; iii++)
251 lhs_contribution(iii,iii) += BDFcoeffs[0] * MassFactors[iii];
261 for(
unsigned int iii = 0; iii<number_of_points; iii++)
263 temp_vec_np[iii] = BDFcoeffs[1]*geom[iii].FastGetSolutionStepValue(rScalarVar,1);
265 for(
unsigned int step = 2;
step<BDFcoeffs.size();
step++)
267 for(
unsigned int iii = 0; iii<number_of_points; iii++)
268 temp_vec_np[iii] += BDFcoeffs[
step]*geom[iii].FastGetSolutionStepValue(rScalarVar,
step);
270 for(
unsigned int iii = 0; iii<number_of_points; iii++)
271 rhs_contribution[iii] -= MassFactors[iii]*temp_vec_np[iii];
277 for(
unsigned int iii = 0; iii<number_of_points; iii++)
278 temp_vec_np[iii] = geom[iii].FastGetSolutionStepValue(rScalarVar);
279 noalias(rhs_contribution) -=
prod(lhs_contribution,temp_vec_np);
282 rhs_contribution *= Volume;
283 lhs_contribution *= Volume;
285 AssembleLHS(mA,lhs_contribution,local_indices);
286 AssembleRHS(mb,rhs_contribution,local_indices);
300 linear_solver->Solve(mA,mDx,mb);
301 std::cout << *(linear_solver) << std::endl;
311 i_dof->GetSolutionStepValue() += mDx[i_dof->EquationId()];
343 double lumping_factor = 1.0/
double(TDim+1.0);
344 const unsigned int number_of_points = TDim+1;
347 for(ModelPart::NodesContainerType::iterator
i =
model_part.NodesBegin();
350 i->FastGetSolutionStepValue(rNodalArea) = 0.0;
351 i->FastGetSolutionStepValue(rProjVar) = 0.0;
356 for(ModelPart::ElementsContainerType::iterator
i =
model_part.ElementsBegin();
367 double nodal_area = Volume*lumping_factor;
369 for(
unsigned int i = 0;
i<number_of_points;
i++)
375 geom[
i].FastGetSolutionStepValue(rNodalArea) += nodal_area;
377 for(
unsigned int j = 0;
j<TDim;
j++)
379 vel_gauss[
j] +=
v[
j] -
w[
j];
382 vel_gauss *= lumping_factor;
385 for(
unsigned int iii = 0; iii<number_of_points; iii++)
386 temp_vec_np[iii] = geom[iii].FastGetSolutionStepValue(rScalarVar);
390 double conv_proj =
inner_prod( u_DN , temp_vec_np);
391 conv_proj *= Volume * lumping_factor;
393 for(
unsigned int iii = 0; iii<number_of_points; iii++)
394 geom[iii].FastGetSolutionStepValue(rProjVar) += conv_proj;
398 for(ModelPart::NodesContainerType::iterator
i =
model_part.NodesBegin();
401 i->FastGetSolutionStepValue(rProjVar) /=
i->FastGetSolutionStepValue(rNodalArea);
414 mA.resize(0,0,
false);
424 unsigned int mEquationSystemSize;
441 unsigned int local_size = LHS_Contribution.size1();
443 for (
unsigned int i_local=0; i_local<
local_size; i_local++)
445 unsigned int i_global=EquationId[i_local];
446 if ( i_global < mEquationSystemSize )
448 for (
unsigned int j_local=0; j_local<
local_size; j_local++)
450 unsigned int j_global=EquationId[j_local];
451 if ( j_global < mEquationSystemSize )
452 A(i_global,j_global) += LHS_Contribution(i_local,j_local);
463 const array_1d<double,TDim+1>& RHS_Contribution,
464 const array_1d<unsigned int,TDim+1>& EquationId
467 unsigned int local_size = RHS_Contribution.size();
469 for (
unsigned int i_local=0; i_local<
local_size; i_local++)
471 unsigned int i_global=EquationId[i_local];
472 if ( i_global < mEquationSystemSize )
475 b[i_global] += RHS_Contribution[i_local];
481 double CalculateTAU(
const array_1d<double,2>&
vel,
const double& Volume)
486 double h = sqrt(2.00*Volume);
489 return h / ( c2*norm_u );
493 double CalculateTAU(
const array_1d<double,3>&
vel,
const double& Volume)
498 double h = pow(6.00*Volume,0.3333333);
501 return h / ( c2*norm_u );
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 class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: global_pointers_vector.h:79
void push_back(TPointerType x)
Definition: global_pointers_vector.h:322
iterator begin()
Definition: global_pointers_vector.h:221
iterator end()
Definition: global_pointers_vector.h:229
This object defines an indexed object.
Definition: indexed_object.h:54
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
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
void clear()
Clear the set, removing all elements.
Definition: pointer_vector_set.h:663
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: pointer_vector_set.h:95
size_type size() const
Returns the number of elements in the container.
Definition: pointer_vector_set.h:502
void push_back(TPointerType x)
Adds a pointer to the end of the set.
Definition: pointer_vector_set.h:544
iterator begin()
Returns an iterator pointing to the beginning of the container.
Definition: pointer_vector_set.h:278
void reserve(int reservedsize)
Reserves memory for a specified number of elements.
Definition: pointer_vector_set.h:733
iterator end()
Returns an iterator pointing to the end of the container.
Definition: pointer_vector_set.h:314
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Definition: pure_convection_tools.h:37
void ClearSystem()
Definition: pure_convection_tools.h:407
PureConvectionUtilities()
Definition: pure_convection_tools.h:45
void ConvectScalarVar(ModelPart &model_part, typename TLinearSolver::Pointer linear_solver, Variable< double > &rScalarVar, const Variable< array_1d< double, 3 > > &rTransportVel, const Variable< array_1d< double, 3 > > &rMeshVel, const Variable< double > &rProjVar, unsigned int time_order)
Definition: pure_convection_tools.h:142
void ConstructSystem(ModelPart &model_part, Variable< double > &rScalarVar, Variable< array_1d< double, 3 > > &rTransportVel, Variable< array_1d< double, 3 > > &rMeshVel)
Definition: pure_convection_tools.h:52
TSparseSpace::VectorType TSystemVectorType
Definition: pure_convection_tools.h:41
PointerVectorSet< Dof< double >, IndexedObject > DofsArrayType
Definition: pure_convection_tools.h:42
TSparseSpace::MatrixType TSystemMatrixType
Definition: pure_convection_tools.h:40
~PureConvectionUtilities()
Definition: pure_convection_tools.h:46
void CalculateProjection(ModelPart &model_part, const Variable< double > &rScalarVar, const Variable< double > &rNodalArea, const Variable< array_1d< double, 3 > > &rTransportVel, const Variable< array_1d< double, 3 > > &rMeshVel, Variable< double > &rProjVar)
Definition: pure_convection_tools.h:324
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_WATCH(variable)
Definition: define.h:806
#define KRATOS_TRY
Definition: define.h:109
dt
Definition: DEM_benchmarks.py:173
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
AMatrix::VectorOuterProductExpression< TExpression1Type, TExpression2Type > outer_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:582
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
int time_order
Definition: ProjectParameters.py:51
int step
Definition: face_heat.py:88
model_part
Definition: face_heat.py:14
v
Definition: generate_convection_diffusion_explicit_element.py:114
w
Definition: generate_convection_diffusion_explicit_element.py:108
tau
Definition: generate_convection_diffusion_explicit_element.py:115
h
Definition: generate_droplet_dynamics.py:91
int local_size
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:17
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
vel
Definition: pure_conduction.py:76
int j
Definition: quadrature.py:648
A
Definition: sensitivityMatrix.py:70
N
Definition: sensitivityMatrix.py:29
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31