/*============================================================================ * General-purpose user-defined functions called before time stepping, at * the end of each time step, and after time-stepping. * * These can be used for operations which do not fit naturally in any other * dedicated user function. *============================================================================*/ /* code_saturne version 8.0-beta */ /* This file is part of code_saturne, a general-purpose CFD tool. Copyright (C) 1998-2023 EDF S.A. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ /*----------------------------------------------------------------------------*/ #include "cs_defs.h" /*---------------------------------------------------------------------------- * Standard C library headers *----------------------------------------------------------------------------*/ #include #include #include #if defined(HAVE_MPI) #include #endif /*---------------------------------------------------------------------------- * PLE library headers *----------------------------------------------------------------------------*/ #include /*---------------------------------------------------------------------------- * Local headers *----------------------------------------------------------------------------*/ #include "cs_headers.h" /*----------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------- * Header for the current file *----------------------------------------------------------------------------*/ #include "cs_prototypes.h" #include #include // For fluid properties #include // For turbulence properties BEGIN_C_DECLS /*----------------------------------------------------------------------------*/ /*! * \file cs_user_extra_operations-boundary_forces.c * * \brief This function is called at the end of each time step, and has a very * general purpose (i.e. anything that does not have another dedicated * user function). * * This is an example of cs_user_extra_operations.c which computes the total * force on a boundary zone. * */ /*----------------------------------------------------------------------------*/ /*============================================================================ * User function definitions *============================================================================*/ /*----------------------------------------------------------------------------*/ /*! * \brief This function is called at the end of each time step. * * It has a very general purpose, although it is recommended to handle * mainly postprocessing or data-extraction type operations. * * \param[in, out] domain pointer to a cs_domain_t structure */ /*----------------------------------------------------------------------------*/ // Skin friction coefficient calculation (Fig. 29 Storms) cs_lnum_t lst_faces_skfr[26]={0}; cs_lnum_t rank_faces_skfr[26]={0}; cs_real_t pos_faces_skfr[78]={0}; cs_real_t area_faces_skfr[78]={0}; int flag_skfr=0; void cs_user_extra_operations(cs_domain_t *domain) { // Net force calculation // Number of boundary faces const cs_lnum_t n_b_faces = domain->mesh->n_b_faces; // Return a pointer to a field based on its name if present. // If no field of the given name is defined, NULL is returned. // All forces of all boundary faces, in an array cs_field_t *b_forces = cs_field_by_name_try("boundary_forces"); //if (b_forces != NULL) { // Initilization of total force vector (total_b_forces) cs_real_3_t total_b_forces = {0., 0., 0.}; // Number and array of boundary faces (to use only certain parts, // containing boundary faces that meet a certain criteria) cs_lnum_t n_elts, *lst_elts; // Allocate memory for n_b_faces elements of type cs_lnum_t, // in lst_elts pointer (to fill a list of boundary faces) BFT_MALLOC(lst_elts, n_b_faces, cs_lnum_t); // Fill a list of boundary faces verifying a given selection criteria // For this case, using a group name of mesh entities, only faces on GTS // (name set at snappyHexMeshDict) // Creo que se pueden usar operadores lógicos cs_selector_get_b_face_list("gtsModel or gtsStrutsPads", &n_elts, lst_elts); //cs_selector_get_b_face_list("gtsModel or gtsStruts or gtsStrutsPads", &n_elts, lst_elts); // Another example: // cs_selector_get_b_face_list // ("normal[0, -1, 0, 0.1] and box[-1000, -1000, -1000, 1000, 0.01, 1000]", // &nlelt, lstelt); // Loop over GTS boundary faces for (cs_lnum_t i_elt = 0; i_elt < n_elts; i_elt++) { cs_lnum_t face_id = lst_elts[i_elt]; // Loop over dimensions for (int ii = 0; ii < 3; ii++) // Sum of forces total_b_forces[ii] += b_forces->val[face_id*3 + ii]; } // Free resources BFT_FREE(lst_elts); /* parallel sum */ cs_parall_sum(3, CS_DOUBLE, total_b_forces); //} // Force coefficients (drag, lift, side force) // Get user parameters defined in the GUI notebook (dynamic pressure and yaw angle) cs_real_t q_inf = cs_notebook_parameter_value_by_name("pressure_dynamic"); cs_real_t yaw_angle = cs_notebook_parameter_value_by_name("yaw_angle"); // Projected area of the GTS in the movement direction // (to be modified later according to the yaw angle) cs_real_t proj_area = -1.0; // Whether struts are considered (or not) in the calculation of the projected area bool strutsInForce = false; // Tolerance used to match floating-point values written in notebook environment double tol_read=1e-6; // WARNING: in case of changing this array, modify accordingly the // structure that takes the value of the projected area (proj_area) cs_real_t yaw_angle_values[]={0.0,2.5,5.0,7.5,10.0,12.5,14.0,-10.0}; int n_yaw_angle_values = sizeof(yaw_angle_values)/sizeof(yaw_angle_values[0]); for (int i=0;i(-tol_read) ){ if(strutsInForce){ switch(i){ case 0: proj_area = cs_notebook_parameter_value_by_name("proj_area_0_0");break; case 1: proj_area = cs_notebook_parameter_value_by_name("proj_area_2_5");break; case 2: proj_area = cs_notebook_parameter_value_by_name("proj_area_5_0");break; case 3: proj_area = cs_notebook_parameter_value_by_name("proj_area_7_5");break; case 4: proj_area = cs_notebook_parameter_value_by_name("proj_area_10_0");break; case 5: proj_area = cs_notebook_parameter_value_by_name("proj_area_12_5");break; case 6: proj_area = cs_notebook_parameter_value_by_name("proj_area_14_0");break; case 7: proj_area = cs_notebook_parameter_value_by_name("proj_area_10_0");break; default: {fprintf(stderr, "Projected area not found. Exiting...\n"); exit(-1);} } } else{ switch(i){ case 0: proj_area = cs_notebook_parameter_value_by_name("proj_area_0_0_wostruts");break; case 1: proj_area = cs_notebook_parameter_value_by_name("proj_area_2_5_wostruts");break; case 2: proj_area = cs_notebook_parameter_value_by_name("proj_area_5_0_wostruts");break; case 3: proj_area = cs_notebook_parameter_value_by_name("proj_area_7_5_wostruts");break; case 4: proj_area = cs_notebook_parameter_value_by_name("proj_area_10_0_wostruts");break; case 5: proj_area = cs_notebook_parameter_value_by_name("proj_area_12_5_wostruts");break; case 6: proj_area = cs_notebook_parameter_value_by_name("proj_area_14_0_wostruts");break; case 7: proj_area = cs_notebook_parameter_value_by_name("proj_area_10_0_wostruts");break; default: {fprintf(stderr, "Projected area not found. Exiting...\n"); exit(-1);} } } break; } } if( (proj_area-(-1.0))(-tol_read) ) {fprintf(stderr, "Projected area not found. Exiting...\n"); exit(-1);} // Force coeficients (WARNING: referred to wind tunnel axis, NOT body axis) cs_real_t c_d = total_b_forces[0]/(q_inf*proj_area); cs_real_t c_s = total_b_forces[1]/(q_inf*proj_area); cs_real_t c_l = total_b_forces[2]/(q_inf*proj_area); // Net static pressure force calculation // Number of boundary faces //const cs_lnum_t n_b_faces = domain->mesh->n_b_faces; // Already defined // Pointer to array containing surface normals of all boundary faces // (L2-norm equals area of the face) const cs_real_t *b_f_face_normal = domain->mesh_quantities->b_f_face_normal; // Initilization of net pressure force vector cs_real_3_t total_b_p_forces = {0., 0., 0.}; // Number and array of boundary faces (to use only certain parts, // containing boundary faces that meet a certain criteria) cs_lnum_t n_elts_p, *lst_elts_p; // Allocate memory for n_b_faces elements of type cs_lnum_t, // in lst_elts pointer (to fill a list of boundary faces) BFT_MALLOC(lst_elts_p, n_b_faces, cs_lnum_t); // Fill a list of boundary faces verifying a given selection criteria // For this case, using a group name of mesh entities, only faces on GTS // (name set at snappyHexMeshDict) // Creo que se pueden usar operadores logicos cs_selector_get_b_face_list("gtsModel or gtsStrutsPads", &n_elts_p, lst_elts_p); //cs_selector_get_b_face_list("gtsModel or gtsStruts or gtsStrutsPads", &n_elts_p, lst_elts_p); /* compute static pressure on selected boundary faces */ // Array to store static pressure values on selected boundary faces cs_real_t *p_b_val; BFT_MALLOC(p_b_val, n_elts_p, cs_real_t); // Compute pressure on a specific boundary region // (A value for each face) cs_post_b_pressure(n_elts_p, lst_elts_p, p_b_val); // Loop over GTS boundary faces for (cs_lnum_t i_elt = 0; i_elt < n_elts_p; i_elt++) { cs_lnum_t face_id = lst_elts_p[i_elt]; for (int ii = 0; ii < 3; ii++) // Sum of the values of forces due to static pressure on each face of // the selected boundary total_b_p_forces[ii] += p_b_val[i_elt]*b_f_face_normal[face_id*3+ii]; } // Free resources BFT_FREE(lst_elts_p); BFT_FREE(p_b_val); /* parallel sum */ cs_parall_sum(3, CS_DOUBLE, total_b_p_forces); // Get current physical time cs_real_t time = cs_glob_time_step->t_cur; // Get current iteration number int nt = cs_glob_time_step->nt_cur; // Manage file to write in FILE *file = NULL; if(cs_glob_rank_id < 1) { file = fopen("user_data_forces.csv", "a"); if(nt==1){ fprintf(file,"iteration, time, drag coeff, lift coeff, side force coeff, net force [x], net force [y], net force [z], net static pressure force [x], net static pressure force [y], net static pressure force [z] \n"); } fprintf(file,"%d, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g \n",nt,time,c_d,c_l,c_s, total_b_forces[0],total_b_forces[1],total_b_forces[2], total_b_p_forces[0],total_b_p_forces[1],total_b_p_forces[2]); if (file != NULL) fclose(file); } // // yplus control over GTS surface // // yplus is a face scalar field cs_field_t *yplus_b = cs_field_by_name_try("yplus"); // Number of boundary faces //const cs_lnum_t n_b_faces = domain->mesh->n_b_faces; // Already defined // Field of boundary faces center of gravity (cog) const cs_real_3_t *cdgfbo = (const cs_real_3_t *)domain->mesh_quantities->b_face_cog; // Field of boundary faces surface (surf) const cs_real_t *surfbn = cs_glob_mesh_quantities->b_face_surf; // Number and array of boundary faces (to use only certain parts, // containing boundary faces that meet a certain criteria) cs_lnum_t n_elts_y, *lst_elts_y; // Allocate memory for n_b_faces elements of type cs_lnum_t, // in lst_elts pointer (to fill a list of boundary faces) BFT_MALLOC(lst_elts_y, n_b_faces, cs_lnum_t); // Fill a list of boundary faces verifying a given selection criteria // For this case, using a group name of mesh entities, only faces on GTS // (name set at snappyHexMeshDict) // Creo que se pueden usar operadores logicos cs_selector_get_b_face_list("gtsModel or gtsStruts or gtsStrutsPads", &n_elts_y, lst_elts_y); // yplus variables definition cs_real_t yplus_min=INFINITY, yplus_max=-INFINITY, yplus_prom=0, yplus_max_cog[3]={-INFINITY,-INFINITY,-INFINITY}, gts_area=0; // Loop over GTS boundary faces for (cs_lnum_t i_elt = 0; i_elt < n_elts_y; i_elt++) { cs_lnum_t face_id = lst_elts_y[i_elt]; // Compute yplus_min yplus_min = cs_math_fmin(yplus_min,yplus_b->val[face_id]); // Compute yplus_max and save position on yplus_max_cog if(yplus_b->val[face_id]>yplus_max) { yplus_max = yplus_b->val[face_id]; for(int j=0;j<3;j++) yplus_max_cog[j]=cdgfbo[face_id][j]; } // Compute yplus_prom yplus_prom += yplus_b->val[face_id]*surfbn[face_id]; gts_area += surfbn[face_id]; } // Reduce variables over procceses cs_parall_min(1,CS_DOUBLE,&yplus_min); cs_parall_max_loc_vals(3,&yplus_max,yplus_max_cog); cs_parall_sum(1,CS_DOUBLE,&yplus_prom); cs_parall_sum(1,CS_DOUBLE,>s_area); yplus_prom = yplus_prom/gts_area; // Manage file to write in FILE *file_yplus = NULL; if(cs_glob_rank_id < 1) { file_yplus = fopen("user_data_yplus.csv", "a"); if(nt==1){ fprintf(file_yplus,"iteration, time, yplus_min, yplus_max, yplus_max_pos_x, yplus_max_pos_y, yplus_max_pos_z, yplus_prom \n"); } fprintf(file_yplus,"%d, %g, %g, %g, %g, %g, %g, %g \n",nt,time,yplus_min,yplus_max, yplus_max_cog[0],yplus_max_cog[1],yplus_max_cog[2],yplus_prom); if (file_yplus != NULL) fclose(file_yplus); } // Free resources BFT_FREE(lst_elts_y); /* // Detecting p0 prescribed position // Number of boundary faces const cs_lnum_t n_b_faces_3 = domain->mesh->n_b_faces; // pressure is a face scalar field cs_field_t *pressure_b = cs_field_by_name_try("pressure"); // Number and array of boundary faces (to use only certain parts, // containing boundary faces that meet a certain criteria) cs_lnum_t n_elts_3, *lst_elts_3; // Allocate memory for n_b_faces elements of type cs_lnum_t, // in lst_elts pointer (to fill a list of boundary faces) BFT_MALLOC(lst_elts_3, n_b_faces_3, cs_lnum_t); // Fill a list of boundary faces verifying a given selection criteria // For this case, using a group name of mesh entities, only faces on GTS // (name set at snappyHexMeshDict) // Creo que se pueden usar operadores lógicos cs_selector_get_b_face_list("outlet", &n_elts_3, lst_elts_3); //cs_selector_get_b_face_list("gtsModel or gtsStruts or gtsStrutsPads", &n_elts, lst_elts); // Another example: // cs_selector_get_b_face_list // ("normal[0, -1, 0, 0.1] and box[-1000, -1000, -1000, 1000, 0.01, 1000]", // &nlelt, lstelt); // Loop over GTS boundary faces for (cs_lnum_t i_elt_3 = 0; i_elt_3 < n_elts_3; i_elt_3++) { cs_lnum_t face_id = lst_elts_3[i_elt_3]; cs_lnum_t tol = 1e-4; if (pressure_b->val[face_id] > -tol && pressure_b->val[face_id] < tol) printf("face id 0 Pa: %d \n",face_id); if (pressure_b->val[face_id] < -5000) printf("face id -5000 Pa: %d \n",face_id); } // Free resources BFT_FREE(lst_elts_3); */ // Skin friction coefficient calculation (Fig. 29 Storms) if (!flag_skfr) { flag_skfr=1; // Llenas faces... cs_real_t input_pos_skfr[]={-2.249733,-2.167852,-2.101718,-2.013538, -1.837179,-1.714357,-1.600984,-1.500207, -1.421475,-1.292355,-1.106548,-1.046712, -0.939637,-0.864054,-0.738084,-0.719188, -0.646755,-0.539679,-0.429455,-0.388514, -0.294036,-0.271991,-0.224752,-0.155468, -0.057841,-0.023199}; cs_lnum_t proc_face_skfr[26]={0}; int size_input_pos_skfr = sizeof(input_pos_skfr)/sizeof(cs_real_t); cs_lnum_t n_elts_skfr, *lst_elts_skfr; BFT_MALLOC(lst_elts_skfr, n_b_faces, cs_lnum_t); cs_selector_get_b_face_list("gtsModel or gtsStrutsPads", &n_elts_skfr, lst_elts_skfr); //printf("point 0"); // Loop over input positions for (cs_lnum_t pos_ind=0; pos_ind < size_input_pos_skfr; pos_ind++){ cs_real_t min_dist=INFINITY; // min_dist_info: face_id, proc, x, y, z cs_real_t min_dist_info[5]={-1,-1,0,0,0}; cs_real_t pos[]={input_pos_skfr[pos_ind],0.0,0.4508}; //if(pos_ind==0) printf("point 1"); // Loop over GTS faces for (cs_lnum_t i_elt_skfr = 0; i_elt_skfr < n_elts_skfr; i_elt_skfr++) { cs_lnum_t face_id = lst_elts_skfr[i_elt_skfr]; //if (i_elt_skfr==5) printf("point 2"); cs_real_t act_dist = (pos[0]-cdgfbo[face_id][0])*(pos[0]-cdgfbo[face_id][0]) + (pos[1]-cdgfbo[face_id][1])*(pos[1]-cdgfbo[face_id][1]) + (pos[2]-cdgfbo[face_id][2])*(pos[2]-cdgfbo[face_id][2]); //if (i_elt_skfr==5) printf("point 3"); if(act_dist < min_dist) { min_dist = act_dist; min_dist_info[0] = face_id; min_dist_info[1] = cs_glob_rank_id; min_dist_info[2] = cdgfbo[face_id][0]; min_dist_info[3] = cdgfbo[face_id][1]; min_dist_info[4] = cdgfbo[face_id][2]; } //if (i_elt_skfr==5) printf("point 4"); } //printf("point 5"); cs_parall_min_loc_vals(5,&min_dist,min_dist_info); //cs_parall_max(26,CS_INT32,proc_face_skfr); //printf("point 6"); if(cs_glob_rank_id < 1) { lst_faces_skfr[pos_ind] = (cs_lnum_t)min_dist_info[0]; rank_faces_skfr[pos_ind] = (cs_lnum_t)min_dist_info[1]; pos_faces_skfr[pos_ind*3+0] = min_dist_info[2]; pos_faces_skfr[pos_ind*3+1] = min_dist_info[3]; pos_faces_skfr[pos_ind*3+2] = min_dist_info[4]; } cs_parall_bcast(0,26,CS_INT32,lst_faces_skfr); cs_parall_bcast(0,26,CS_INT32,rank_faces_skfr); //printf("point 7"); } //printf("point 8"); //printf("point 9"); //fprintf("pos0 %d, pos5 %d, pos20 %d \n",proc_face_skfr); FILE *file_skfr = NULL; if(cs_glob_rank_id < 1) { file_skfr = fopen("user_data_skfr.csv", "a"); if(nt==1){ for(int i=0;i<26;i++) { int id = lst_faces_skfr[i]; int proc = rank_faces_skfr[i]; fprintf(file_skfr,"face_id %d proc = %d, x = %g, y = %g, z = %g \n",id,proc,pos_faces_skfr[i*3+0],pos_faces_skfr[i*3+1],pos_faces_skfr[i*3+2]); } } if (file_skfr != NULL) fclose(file_skfr); } BFT_FREE(lst_elts_skfr); } // Usas faces... if((nt-199)%200==0){ const cs_real_3_t *surfbo = (const cs_real_3_t *)cs_glob_mesh_quantities->b_face_normal; const cs_real_t *surfbnval = cs_glob_mesh_quantities->b_face_surf; //printf("point 1 \n"); cs_real_t tot_stress_tang[78] = {0}; cs_field_t *b_forces = cs_field_by_name_try("boundary_forces"); //if (stress_skfr==NULL) printf("No lo reconoce"); //cs_lnum_t n_elts_skfr, *lst_elts_skfr; //cs_selector_get_b_face_list("gtsModel or gtsStrutsPads", &n_elts_skfr, lst_elts_skfr); //printf("point 2 \n"); /* cs_real_3_t total_b_stresses = {0., 0., 0.}; cs_lnum_t n_elts_str, *lst_elts_str; BFT_MALLOC(lst_elts_str, n_b_faces, cs_lnum_t); cs_selector_get_b_face_list("gtsModel or gtsStrutsPads", &n_elts_str, lst_elts_str); for (cs_lnum_t i_elt = 0; i_elt < n_elts_str; i_elt++) { cs_lnum_t face_id = lst_elts_str[i_elt]; for (int ii = 0; ii < 3; ii++) total_b_stresses[ii] += stress_skfr->val[face_id*3 + ii]; } BFT_FREE(lst_elts_str); cs_parall_sum(3, CS_DOUBLE, total_b_stresses); */ //cs_real_3_t *stresses; //BFT_MALLOC(stresses, 26, cs_real_3_t); //cs_post_stress_tangential(26, lst_faces_skfr, stresses); for(int i=0; i<26; i++) { if(cs_glob_rank_id == rank_faces_skfr[i]){ int face_id = lst_faces_skfr[i]; //printf("face_id: %d, rank: %d",face_id,cs_glob_rank_id); cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; //srfbn = sqrt(surfbo[face_id][0]*surfbo[face_id][0] + surfbo[face_id][1]*surfbo[face_id][1] + surfbo[face_id][2]*surfbo[face_id][2]); area_faces_skfr[i*3+0] = surfbo[face_id][0]; area_faces_skfr[i*3+1] = surfbo[face_id][1]; area_faces_skfr[i*3+2] = surfbo[face_id][2]; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; /// printf("point 3 \n"); //if (cs_glob_rank_id == rank_faces_skfr[i]) { tot_stress_tang[i*3+0] = (b_forces->val[face_id*3+0] - fornor*srfnor[0]) / srfbn; tot_stress_tang[i*3+1] = (b_forces->val[face_id*3+1] - fornor*srfnor[1]) / srfbn; tot_stress_tang[i*3+2] = (b_forces->val[face_id*3+2] - fornor*srfnor[2]) / srfbn; //tot_stress_tang[i*3+0] = b_forces->val[face_id*3+0]; //tot_stress_tang[i*3+1] = b_forces->val[face_id*3+1]; //tot_stress_tang[i*3+2] = b_forces->val[face_id*3+2]; //} } cs_parall_bcast(rank_faces_skfr[i],78,CS_DOUBLE,tot_stress_tang); cs_parall_bcast(rank_faces_skfr[i],78,CS_DOUBLE,area_faces_skfr); } //cs_parall_sum(78,CS_DOUBLE,tot_stress_tang); //printf("point 5 \n"); FILE *file_stress = NULL; if(cs_glob_rank_id < 1) { // printf("point 6 \n"); file_stress = fopen("user_data_stress.csv", "a"); for(int i=0;i<26;i++) { cs_real_t stress_tang_norm_act = sqrt(tot_stress_tang[i*3+0]*tot_stress_tang[i*3+0] + tot_stress_tang[i*3+1]*tot_stress_tang[i*3+1] + tot_stress_tang[i*3+2]*tot_stress_tang[i*3+2]); cs_real_t skfr_coeff_act = stress_tang_norm_act/q_inf; //cs_real_t skfr_coeff_act = stress_tang_norm_act; //fprintf(file_stress,"sx = %g, sy = %g, sz = %g \n",tot_stress_tang[i*3+0]),tot_stress_tang[i*3+1],tot_stress_tang[i*3+2]; if (i<25) {fprintf(file_stress,"%g, ",skfr_coeff_act);} else {fprintf(file_stress,"%g \n",skfr_coeff_act);} } //fprintf(file_stress,"Areas: \n"); //for(int i=0;i<26;i++) { //if (i<25) {fprintf(file_stress,"%g %g %g, ",area_faces_skfr[i*3+0],area_faces_skfr[i*3+1],area_faces_skfr[i*3+2]);} //else {fprintf(file_stress,"%g %g %g\n",area_faces_skfr[i*3+0],area_faces_skfr[i*3+1],area_faces_skfr[i*3+2]);} //} if (file_stress != NULL) fclose(file_stress); } //cs_real_3_t *stresses; //BFT_MALLOC(stresses, 26, cs_real_3_t); //cs_post_stress_tangential(26, lst_faces_skfr, stresses); // Number and array of boundary faces (to use only certain parts, // containing boundary faces that meet a certain criteria) cs_lnum_t n_elts_skfr_top, *lst_elts_skfr_top; cs_lnum_t n_elts_skfr_bottom, *lst_elts_skfr_bottom; cs_lnum_t n_elts_skfr_side, *lst_elts_skfr_side; cs_lnum_t n_elts_skfr_front, *lst_elts_skfr_front; cs_lnum_t n_elts_skfr_base, *lst_elts_skfr_base; cs_lnum_t n_elts_skfr_all, *lst_elts_skfr_all; // Allocate memory for n_b_faces elements of type cs_lnum_t, // in lst_elts pointer (to fill a list of boundary faces) BFT_MALLOC(lst_elts_skfr_top, n_b_faces, cs_lnum_t); BFT_MALLOC(lst_elts_skfr_bottom, n_b_faces, cs_lnum_t); BFT_MALLOC(lst_elts_skfr_side, n_b_faces, cs_lnum_t); BFT_MALLOC(lst_elts_skfr_front, n_b_faces, cs_lnum_t); BFT_MALLOC(lst_elts_skfr_base, n_b_faces, cs_lnum_t); BFT_MALLOC(lst_elts_skfr_all, n_b_faces, cs_lnum_t); // Fill a list of boundary faces verifying a given selection criteria // For this case, using a group name of mesh entities, only faces on GTS // (name set at snappyHexMeshDict) // Creo que se pueden usar operadores lógicos cs_selector_get_b_face_list("gtsModel and box[-2.475, -0.1618, 0.4500, -0.001, 0.1618, 0.4520]", &n_elts_skfr_top, lst_elts_skfr_top); cs_selector_get_b_face_list("(gtsModel or gtsStrutsPads) and box[-2.475, -0.1618, -0.004, -0.001, 0.1618, 0.001]", &n_elts_skfr_bottom, lst_elts_skfr_bottom); cs_selector_get_b_face_list("gtsModel and box[-2.475, -0.1622, 0.001, -0.001, -0.1618, 0.449]", &n_elts_skfr_side, lst_elts_skfr_side); cs_selector_get_b_face_list("gtsModel and box[-2.478, -0.1618, 0.001, -1.82878, 0.1618, 0.449]", &n_elts_skfr_front, lst_elts_skfr_front); cs_selector_get_b_face_list("gtsModel and box[-0.004, -0.1618, 0.001, 0.004, 0.1618, 0.449]", &n_elts_skfr_base, lst_elts_skfr_base); cs_selector_get_b_face_list("gtsModel or gtsStrutsPads", &n_elts_skfr_all, lst_elts_skfr_all); //cs_selector_get_b_face_list("gtsModel or gtsStruts or gtsStrutsPads", &n_elts, lst_elts); // Another example: // cs_selector_get_b_face_list // ("normal[0, -1, 0, 0.1] and box[-1000, -1000, -1000, 1000, 0.01, 1000]", // &nlelt, lstelt); cs_real_t sigma_act_top[3]={0}; cs_real_t sigma_act_bottom[3]={0}; cs_real_t sigma_act_side[3]={0}; cs_real_t sigma_act_front[3]={0}; cs_real_t sigma_act_base[3]={0}; cs_real_t sigma_act_all[3]={0}; cs_real_t area_top=0, area_bottom=0, area_side=0, area_front=0, area_base=0, area_all=0; cs_real_t cf_top, cf_bottom, cf_side, cf_front, cf_base, cf_all; //if(cs_glob_rank_id < 1) printf("point 1\n"); // Loop over GTS boundary faces (top) for (cs_lnum_t i_elt = 0; i_elt < n_elts_skfr_top; i_elt++) { int face_id = lst_elts_skfr_top[i_elt]; cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; area_top += srfbn; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; sigma_act_top[0] += (b_forces->val[face_id*3+0] - fornor*srfnor[0]); sigma_act_top[1] += (b_forces->val[face_id*3+1] - fornor*srfnor[1]); sigma_act_top[2] += (b_forces->val[face_id*3+2] - fornor*srfnor[2]); } //if(cs_glob_rank_id < 1) printf("point 2\n"); cs_parall_sum(3, CS_DOUBLE, sigma_act_top); //if(cs_glob_rank_id < 1) printf("point 2.1\n"); cs_parall_sum(1, CS_DOUBLE, &area_top); //if(cs_glob_rank_id < 1) printf("point 2.2\n"); if(cs_glob_rank_id < 1) { for(int i=0;i<3;i++) sigma_act_top[i] = sigma_act_top[i]/proj_area; //if(cs_glob_rank_id < 1) printf("point 2.3\n"); cs_real_t norm_sigma_top = sqrt(sigma_act_top[0]*sigma_act_top[+0] + sigma_act_top[1]*sigma_act_top[1] + sigma_act_top[2]*sigma_act_top[2]); //if(cs_glob_rank_id < 1) printf("point 2.4\n"); cf_top = norm_sigma_top/q_inf; } //if(cs_glob_rank_id < 1) printf("point 3\n"); // Loop over GTS boundary faces (bottom) for (cs_lnum_t i_elt = 0; i_elt < n_elts_skfr_bottom; i_elt++) { int face_id = lst_elts_skfr_bottom[i_elt]; cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; area_bottom += srfbn; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; sigma_act_bottom[0] += (b_forces->val[face_id*3+0] - fornor*srfnor[0]); sigma_act_bottom[1] += (b_forces->val[face_id*3+1] - fornor*srfnor[1]); sigma_act_bottom[2] += (b_forces->val[face_id*3+2] - fornor*srfnor[2]); } cs_parall_sum(3, CS_DOUBLE, sigma_act_bottom); cs_parall_sum(1, CS_DOUBLE, &area_bottom); if(cs_glob_rank_id < 1) { for(int i=0;i<3;i++) sigma_act_bottom[i] = sigma_act_bottom[i]/proj_area; cs_real_t norm_sigma_bottom = sqrt(sigma_act_bottom[0]*sigma_act_bottom[+0] + sigma_act_bottom[1]*sigma_act_bottom[1] + sigma_act_bottom[2]*sigma_act_bottom[2]); cf_bottom = norm_sigma_bottom/q_inf; } // Loop over GTS boundary faces (side) for (cs_lnum_t i_elt = 0; i_elt < n_elts_skfr_side; i_elt++) { int face_id = lst_elts_skfr_side[i_elt]; cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; area_side += srfbn; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; sigma_act_side[0] += (b_forces->val[face_id*3+0] - fornor*srfnor[0]); sigma_act_side[1] += (b_forces->val[face_id*3+1] - fornor*srfnor[1]); sigma_act_side[2] += (b_forces->val[face_id*3+2] - fornor*srfnor[2]); } cs_parall_sum(3, CS_DOUBLE, sigma_act_side); cs_parall_sum(1, CS_DOUBLE, &area_side); if(cs_glob_rank_id < 1) { for(int i=0;i<3;i++) sigma_act_side[i] = sigma_act_side[i]/proj_area; cs_real_t norm_sigma_side = sqrt(sigma_act_side[0]*sigma_act_side[+0] + sigma_act_side[1]*sigma_act_side[1] + sigma_act_side[2]*sigma_act_side[2]); cf_side = norm_sigma_side/q_inf; } // Loop over GTS boundary faces (front) for (cs_lnum_t i_elt = 0; i_elt < n_elts_skfr_front; i_elt++) { int face_id = lst_elts_skfr_front[i_elt]; cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; area_front += srfbn; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; sigma_act_front[0] += (b_forces->val[face_id*3+0] - fornor*srfnor[0]); sigma_act_front[1] += (b_forces->val[face_id*3+1] - fornor*srfnor[1]); sigma_act_front[2] += (b_forces->val[face_id*3+2] - fornor*srfnor[2]); } cs_parall_sum(3, CS_DOUBLE, sigma_act_front); cs_parall_sum(1, CS_DOUBLE, &area_front); if(cs_glob_rank_id < 1) { for(int i=0;i<3;i++) sigma_act_front[i] = sigma_act_front[i]/proj_area; cs_real_t norm_sigma_front = sqrt(sigma_act_front[0]*sigma_act_front[+0] + sigma_act_front[1]*sigma_act_front[1] + sigma_act_front[2]*sigma_act_front[2]); cf_front = norm_sigma_front/q_inf; } // Loop over GTS boundary faces (base) for (cs_lnum_t i_elt = 0; i_elt < n_elts_skfr_base; i_elt++) { int face_id = lst_elts_skfr_base[i_elt]; cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; area_base += srfbn; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; sigma_act_base[0] += (b_forces->val[face_id*3+0] - fornor*srfnor[0]); sigma_act_base[1] += (b_forces->val[face_id*3+1] - fornor*srfnor[1]); sigma_act_base[2] += (b_forces->val[face_id*3+2] - fornor*srfnor[2]); } cs_parall_sum(3, CS_DOUBLE, sigma_act_base); cs_parall_sum(1, CS_DOUBLE, &area_base); if(cs_glob_rank_id < 1) { for(int i=0;i<3;i++) sigma_act_base[i] = sigma_act_base[i]/proj_area; cs_real_t norm_sigma_base = sqrt(sigma_act_base[0]*sigma_act_base[+0] + sigma_act_base[1]*sigma_act_base[1] + sigma_act_base[2]*sigma_act_base[2]); cf_base = norm_sigma_base/q_inf; } for (cs_lnum_t i_elt = 0; i_elt < n_elts_skfr_all; i_elt++) { int face_id = lst_elts_skfr_all[i_elt]; cs_real_t srfbn,srfnor[3], fornor; srfbn = surfbnval[face_id]; area_all += srfbn; srfnor[0] = surfbo[face_id][0] / srfbn; srfnor[1] = surfbo[face_id][1] / srfbn; srfnor[2] = surfbo[face_id][2] / srfbn; fornor = b_forces->val[face_id*3+0]*srfnor[0] + b_forces->val[face_id*3+1]*srfnor[1] + b_forces->val[face_id*3+2]*srfnor[2]; sigma_act_all[0] += (b_forces->val[face_id*3+0] - fornor*srfnor[0]); sigma_act_all[1] += (b_forces->val[face_id*3+1] - fornor*srfnor[1]); sigma_act_all[2] += (b_forces->val[face_id*3+2] - fornor*srfnor[2]); } cs_parall_sum(3, CS_DOUBLE, sigma_act_all); cs_parall_sum(1, CS_DOUBLE, &area_all); if(cs_glob_rank_id < 1) { for(int i=0;i<3;i++) sigma_act_all[i] = sigma_act_all[i]/proj_area; cs_real_t norm_sigma_all = sqrt(sigma_act_all[0]*sigma_act_all[+0] + sigma_act_all[1]*sigma_act_all[1] + sigma_act_all[2]*sigma_act_all[2]); cf_all = norm_sigma_all/q_inf; } // Save to txt //if(cs_glob_rank_id < 1) printf("point 4\n"); FILE *file_cfparts = NULL; if(cs_glob_rank_id < 1) { file_cfparts = fopen("user_data_cfparts.csv", "a"); // cftop cfbottom cfside cffront cfbase fprintf(file_cfparts,"%g %g %g %g %g %g %g\n",cf_all,cf_top+cf_bottom+cf_side*2.0+cf_front+cf_base,cf_top,cf_bottom,cf_side,cf_front,cf_base); if (file_cfparts != NULL) fclose(file_cfparts); } //if(cs_glob_rank_id < 1) printf("point 5\n"); BFT_FREE(lst_elts_skfr_top); BFT_FREE(lst_elts_skfr_bottom); BFT_FREE(lst_elts_skfr_side); BFT_FREE(lst_elts_skfr_front); BFT_FREE(lst_elts_skfr_base); } } END_C_DECLS