1 #if !defined(KRATOS_CALCULATE_WATER_FRACTION_UTILITY_INCLUDED )
2 #define KRATOS_CALCULATE_WATER_FRACTION_UTILITY_INCLUDED
13 #include "utilities/geometry_utilities.h"
24 template<
unsigned int TDim>
73 double sum_areas=1.0e-100;
77 ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.
ElementsBegin();
80 vector<unsigned int> element_partition;
82 int number_of_threads = omp_get_max_threads();
84 int number_of_threads = 1;
86 OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.
Elements().size(), element_partition);
88 #pragma omp parallel for reduction(+:sum_areas)
89 for(
int kkk=0; kkk<number_of_threads; kkk++)
91 double thread_sum_areas=0.0;
92 for(
unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
94 ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
101 int negative_nodes=0;
102 int positive_nodes=0;
103 for (
unsigned int k = 0;
k < (TDim+1);
k++)
105 if(geom[
k].FastGetSolutionStepValue(DISTANCE)<0.0)
110 if (negative_nodes==(TDim+1))
111 thread_sum_areas+=Area;
113 else if (negative_nodes>0)
116 for (
unsigned int i = 0;
i < (TDim+1);
i++)
118 distances[
i] = geom[
i].FastGetSolutionStepValue(DISTANCE);
127 std::vector< Matrix > gauss_gradients(3*(TDim-1));
131 for (
unsigned int i = 0;
i < (TDim+1);
i++)
134 for (
unsigned int j = 0;
j < TDim;
j++)
137 for (
unsigned int i = 0;
i < 3*(TDim-1);
i++)
138 gauss_gradients[
i].resize(2, TDim,
false);
140 for (
unsigned int i=0;
i!=ndivisions;
i++)
141 if (signs(
i)<0.0) thread_sum_areas+=volumes(
i);
146 sum_areas = thread_sum_areas;
158 double all_threads_water_height=-100000000.0;
159 const double tolerance=0.001;
160 const double upper_limit=x_position+tolerance;
161 const double lower_limit=x_position-tolerance;
162 ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.
NodesBegin();
164 vector<unsigned int> node_partition;
166 int number_of_threads = omp_get_max_threads();
168 int number_of_threads = 1;
170 OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.
Nodes().size(), node_partition);
172 #pragma omp parallel for
173 for(
int kkk=0; kkk<number_of_threads; kkk++)
175 double local_thread_water_height=-100000000.0;
176 for(
unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
178 ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
179 double & distance = (inode->FastGetSolutionStepValue(DISTANCE));
182 if ((inode->X())<upper_limit && (inode->X())>lower_limit && (inode->Y())>local_thread_water_height)
183 local_thread_water_height = (inode->Y());
190 if ( local_thread_water_height > all_threads_water_height )
194 if ( local_thread_water_height > all_threads_water_height ) all_threads_water_height = local_thread_water_height;
199 return all_threads_water_height;
210 const double delta_t = CurrentProcessInfo[DELTA_TIME];
212 double all_threads_max_courant = 0.0;
213 ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.
ElementsBegin();
215 vector<unsigned int> element_partition;
217 int number_of_threads = omp_get_max_threads();
219 int number_of_threads = 1;
221 OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.
Elements().size(), element_partition);
223 #pragma omp parallel for
224 for(
int kkk=0; kkk<number_of_threads; kkk++)
226 double local_thread_max_courant = 0.0;
227 for(
unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
229 ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
231 if ((ielem->GetValue(VELOCITY_OVER_ELEM_SIZE))>local_thread_max_courant)
232 local_thread_max_courant = (ielem->GetValue(VELOCITY_OVER_ELEM_SIZE));
236 if ( local_thread_max_courant > all_threads_max_courant )
240 if ( local_thread_max_courant > all_threads_max_courant ) all_threads_max_courant = local_thread_max_courant;
245 all_threads_max_courant *=
delta_t * 1.414;
247 return all_threads_max_courant;
258 const double delta_t = CurrentProcessInfo[DELTA_TIME];
260 double all_threads_max_courant = 0.0;
261 ModelPart::NodesContainerType::iterator inodebegin = mr_model_part.
NodesBegin();
263 vector<unsigned int> node_partition;
265 int number_of_threads = omp_get_max_threads();
267 int number_of_threads = 1;
269 OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.
Nodes().size(), node_partition);
271 #pragma omp parallel for
272 for(
int kkk=0; kkk<number_of_threads; kkk++)
274 double local_thread_max_courant = 0.0;
275 for(
unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
277 ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
278 const double & distance = (inode->FastGetSolutionStepValue(DISTANCE));
279 const double velocity = sqrt(pow(inode->FastGetSolutionStepValue(VELOCITY_X),2)+pow(inode->FastGetSolutionStepValue(VELOCITY_Y),2)+pow(inode->FastGetSolutionStepValue(VELOCITY_Z),2));
280 const double nodal_courant = (
velocity*
delta_t/inode->FastGetSolutionStepValue(MEAN_SIZE));
281 if(nodal_courant>local_thread_max_courant && distance < 0.0)
282 local_thread_max_courant = nodal_courant;
286 if ( local_thread_max_courant > all_threads_max_courant )
290 if ( local_thread_max_courant > all_threads_max_courant ) all_threads_max_courant = local_thread_max_courant;
297 return all_threads_max_courant;
307 const double delta_t = CurrentProcessInfo[DELTA_TIME];
311 double sum_courant=0.0;
313 ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.
ElementsBegin();
315 vector<unsigned int> element_partition;
317 int number_of_threads = omp_get_max_threads();
319 int number_of_threads = 1;
321 OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.
Elements().size(), element_partition);
323 #pragma omp parallel for reduction(+:sum_courant)
324 for(
int kkk=0; kkk<number_of_threads; kkk++)
326 double thread_sum_courant=0.0;
327 for(
unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
329 ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
331 if ((ielem->GetValue(VELOCITY_OVER_ELEM_SIZE))>0.0)
332 thread_sum_courant += ielem->GetValue(VELOCITY_OVER_ELEM_SIZE);
334 sum_courant += thread_sum_courant;
350 double viscosity = CurrentProcessInfo[VISCOSITY];
356 double nodal_weight = 1.0/ (
double (TDim) );
360 ModelPart::ElementsContainerType::iterator ielembegin = mr_model_part.
ElementsBegin();
363 vector<unsigned int> element_partition;
365 int number_of_threads = omp_get_max_threads();
367 int number_of_threads = 1;
369 OpenMPUtils::CreatePartition(number_of_threads, mr_model_part.
Elements().size(), element_partition);
371 #pragma omp parallel for reduction(+:force)
372 for(
int kkk=0; kkk<number_of_threads; kkk++)
374 double thread_force=0.0;
375 for(
unsigned int ii=element_partition[kkk]; ii<element_partition[kkk+1]; ii++)
377 ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
378 if (ielem->Is(ACTIVE))
385 unsigned int number_of_fixed_nodes=0;
389 for (
unsigned int i = 0;
i < (TDim+1);
i++)
392 for (
unsigned int j = 0;
j < (TDim);
j++)
395 if constexpr (TDim==2)
397 if (geom[
i].IsFixed(FRACT_VEL_X) && geom[
i].IsFixed(FRACT_VEL_Y))
399 fixed_nodes[number_of_fixed_nodes]=
i;
400 is_node_fixed[
i]=
true;
401 number_of_fixed_nodes++;
404 is_node_fixed[
i]=
false;
408 if (geom[
i].IsFixed(FRACT_VEL_X) && geom[
i].IsFixed(FRACT_VEL_Y) && geom[
i].IsFixed(FRACT_VEL_Z) )
410 fixed_nodes[number_of_fixed_nodes]=
i;
411 number_of_fixed_nodes++;
412 is_node_fixed[
i]=
true;
415 is_node_fixed[
i]=
false;
419 double plane_point_distance=1.0;
420 double fixed_face_area_or_lenght=0.0;
422 if (number_of_fixed_nodes==TDim)
426 unsigned int free_node=0;
427 if constexpr (TDim==2)
429 fixed_face_area_or_lenght = fabs(sqrt(pow((geom[fixed_nodes[1]].
Y()-geom[fixed_nodes[0]].
Y()),2 ) + pow( (geom[fixed_nodes[1]].
X()-geom[fixed_nodes[0]].
X() ),2 ) ) );
430 normal[0] = geom[fixed_nodes[1]].Y()-geom[fixed_nodes[0]].Y();
431 normal[1] = - ( geom[fixed_nodes[1]].X()-geom[fixed_nodes[0]].X() );
433 normal /= sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
435 if (fixed_nodes[0]==0)
437 if (fixed_nodes[1]==1)
446 plane_point_distance =
inner_prod( (geom[free_node].Coordinates()-geom[fixed_nodes[0]].Coordinates()) , normal);
448 if (plane_point_distance<0.0)
450 plane_point_distance*=-1.0;
457 MathUtils<double>::CrossProduct(normal, geom[fixed_nodes[1]].Coordinates() - geom[fixed_nodes[0]].Coordinates(), geom[fixed_nodes[2]].Coordinates() - geom[fixed_nodes[0]].Coordinates() );
458 fixed_face_area_or_lenght = 0.5 * sqrt( pow(normal[0],2) + pow(normal[1],2) + pow(normal[2],2) );
459 normal /= 2.0 * fixed_face_area_or_lenght;
463 for (
unsigned int j=0;
j!=(TDim+1);
j++)
465 if (is_node_fixed[
j]==
false)
473 plane_point_distance =
inner_prod( (geom[free_node].Coordinates()-geom[fixed_nodes[0]].Coordinates()) , normal);
474 if (plane_point_distance<0.0)
477 plane_point_distance*=-1.0;
483 boundary_stress = - geom[free_node].FastGetSolutionStepValue(VELOCITY)*
viscosity/(fabs(plane_point_distance));
489 thread_force += boundary_stress[direction]*fixed_face_area_or_lenght;
493 for (
unsigned int j=0;
j!=(TDim);
j++)
501 thread_force +=nodal_weight*(geom[fixed_nodes[
j]].FastGetSolutionStepValue(PRESSURE))*fixed_face_area_or_lenght*normal[direction];
Definition: calculate_water_fraction.h:26
~CalculateWaterFraction()
Definition: calculate_water_fraction.h:40
double CalculateForce(int direction)
Definition: calculate_water_fraction.h:345
KRATOS_CLASS_POINTER_DEFINITION(CalculateWaterFraction)
double CalculateMaxCourant()
Definition: calculate_water_fraction.h:205
double Calculate()
Definition: calculate_water_fraction.h:68
double CalculateMeanCourant()
Definition: calculate_water_fraction.h:302
double CalculateWaterHeight(double x_position)
Definition: calculate_water_fraction.h:154
CalculateWaterFraction(ModelPart &model_part)
Definition: calculate_water_fraction.h:31
double CalculateMaxCourantInNegativeElements()
Definition: calculate_water_fraction.h:252
static int CalculateEnrichedShapeFuncions(BoundedMatrix< double, 4, 3 > &rPoints, BoundedMatrix< double, 4, 3 > &DN_DX, array_1d< double, 4 > &rDistances, array_1d< double, 6 > &rVolumes, BoundedMatrix< double, 6, 4 > &rShapeFunctionValues, array_1d< double, 6 > &rPartitionsSign, std::vector< Matrix > &rGradientsValue, BoundedMatrix< double, 6, 2 > &NEnriched)
Definition: enrichment_utilities.h:500
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
Definition: amatrix_interface.h:41
static T CrossProduct(const T &a, const T &b)
Performs the vector product of the two input vectors a,b.
Definition: math_utils.h:762
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Short class definition.
Definition: array_1d.h:61
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_TRY
Definition: define.h:109
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
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
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
float velocity
Definition: PecletTest.py:54
float viscosity
Definition: edgebased_var.py:8
model_part
Definition: face_heat.py:14
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
float delta_t
Definition: rotatingcone_PureConvectionBenchmarking.py:129
N
Definition: sensitivityMatrix.py:29
integer i
Definition: TensorModule.f:17