1 /*============================================================================
2 * Management of the GUI parameters file: mobile mesh
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <stdio.h>
34 #include <stdlib.h>
35 #include <math.h>
36 #include <string.h>
37 #include <fcntl.h>
38 #include <unistd.h>
39 #include <assert.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_error.h"
47 #include "bft_printf.h"
48
49 #include "cs_ale.h"
50 #include "cs_base.h"
51 #include "cs_boundary.h"
52 #include "cs_boundary_zone.h"
53 #include "cs_convection_diffusion.h"
54 #include "cs_field_pointer.h"
55 #include "cs_gui.h"
56 #include "cs_gui_util.h"
57 #include "cs_gui_boundary_conditions.h"
58 #include "cs_mesh.h"
59 #include "cs_prototypes.h"
60 #include "cs_timer.h"
61 #include "cs_time_step.h"
62 #include "cs_volume_zone.h"
63
64 /*----------------------------------------------------------------------------
65 * Header for the current file
66 *----------------------------------------------------------------------------*/
67
68 #include "cs_gui_mobile_mesh.h"
69
70 /*----------------------------------------------------------------------------*/
71
72 BEGIN_C_DECLS
73
74 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
75
76 /*=============================================================================
77 * Local Macro Definitions
78 *============================================================================*/
79
80 /* debugging switch */
81 #define _XML_DEBUG_ 0
82
83 /*============================================================================
84 * Static variables
85 *============================================================================*/
86
87 /*-----------------------------------------------------------------------------
88 * Possible values for boundary nature
89 *----------------------------------------------------------------------------*/
90
91 enum ale_boundary_nature
92 {
93 ale_boundary_nature_none,
94 ale_boundary_nature_fixed_wall,
95 ale_boundary_nature_sliding_wall,
96 ale_boundary_nature_internal_coupling,
97 ale_boundary_nature_external_coupling,
98 ale_boundary_nature_fixed_velocity,
99 ale_boundary_nature_fixed_displacement,
100 ale_boundary_nature_free_surface
101 };
102
103 /*============================================================================
104 * Private function definitions
105 *============================================================================*/
106
107 /*-----------------------------------------------------------------------------
108 * Return value for ALE method
109 *
110 * parameters:
111 * tn_ale <-- tree node associated with ALE
112 *
113 * return:
114 * 0 for isotropic, 1 for orthotropic
115 *----------------------------------------------------------------------------*/
116
117 static int
_ale_visc_type(cs_tree_node_t * tn_ale)118 _ale_visc_type(cs_tree_node_t *tn_ale)
119 {
120 int mvisc_type = 0;
121
122 cs_tree_node_t *tn_mv = cs_tree_get_node(tn_ale, "mesh_viscosity");
123
124 const char *type = cs_tree_node_get_tag(tn_mv, "type");
125 if (type != NULL) {
126 if (strcmp(type, "isotrop") != 0) {
127 if (strcmp(type, "orthotrop") == 0)
128 mvisc_type = 1;
129 else
130 bft_error(__FILE__, __LINE__, 0,
131 "invalid mesh viscosity type: %s", type);
132 }
133 }
134
135 return mvisc_type;
136 }
137
138 /*-----------------------------------------------------------------------------
139 * Get the ale boundary formula
140 *
141 * parameters:
142 * tn_w <-- pointer to tree node for a given BC
143 * choice <-- nature: "fixed_velocity" or "fixed_displacement"
144 *----------------------------------------------------------------------------*/
145
146 static const char*
_get_ale_boundary_formula(cs_tree_node_t * tn_w,const char * choice)147 _get_ale_boundary_formula(cs_tree_node_t *tn_w,
148 const char *choice)
149 {
150 cs_tree_node_t *tn = cs_tree_get_node(tn_w, "ale");
151 tn = cs_tree_node_get_sibling_with_tag(tn, "choice", choice);
152
153 const char *formula = cs_tree_node_get_child_value_str(tn, "formula");
154
155 return formula;
156 }
157
158 /*-----------------------------------------------------------------------------
159 * Get uialcl data for fixed displacement
160 *
161 * parameters:
162 * tn_w <-- pointer to tree node for a given mobile_boundary BC
163 * z <-- selected boundary zone
164 * impale --> impale
165 * disale --> disale
166 *----------------------------------------------------------------------------*/
167
168 static void
_uialcl_fixed_displacement(cs_tree_node_t * tn_w,const cs_zone_t * z,int * impale,cs_real_t disale[][3])169 _uialcl_fixed_displacement(cs_tree_node_t *tn_w,
170 const cs_zone_t *z,
171 int *impale,
172 cs_real_t disale[][3])
173 {
174
175 const cs_mesh_t *m = cs_glob_mesh;
176
177 /* Get formula */
178 const char *formula = _get_ale_boundary_formula(tn_w, "fixed_displacement");
179
180 if (!formula)
181 bft_error(__FILE__, __LINE__, 0,
182 _("Boundary nature formula is null for %s."),
183 cs_gui_node_get_tag(tn_w, "label"));
184
185 /* Evaluate formula using meg */
186 cs_real_t *bc_vals = cs_meg_boundary_function(z,
187 "mesh_velocity",
188 "fixed_displacement");
189
190 /* Loop over boundary faces */
191 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
192 const cs_lnum_t face_id = z->elt_ids[elt_id];
193
194 /* ALE BC on vertices */
195 const cs_lnum_t s = m->b_face_vtx_idx[face_id];
196 const cs_lnum_t e = m->b_face_vtx_idx[face_id+1];
197
198 /* Compute the portion of surface associated to v_id_1 */
199 for (cs_lnum_t k = s; k < e; k++) {
200
201 const cs_lnum_t v_id = m->b_face_vtx_lst[k];
202 cs_real_t *_val = disale[v_id];
203
204 impale[v_id] = 1;
205 //FIXME prorata
206 for (int d = 0; d < 3; d++)
207 _val[d] = bc_vals[elt_id + d * z->n_elts];
208
209 }
210 }
211 }
212
213 /*-----------------------------------------------------------------------------
214 * Get uialcl data for fixed velocity
215 *
216 * parameters:
217 * tn_w <-- pointer to tree node for a given mobile_boundary BC
218 * iuma <-- iuma
219 * ivma <-- ivma
220 * iwma <-- iwma
221 * nfabor <-- Number of boundary faces
222 * z <-- selected boundary zone
223 * rcodcl --> rcodcl
224 *----------------------------------------------------------------------------*/
225
226 static void
_uialcl_fixed_velocity(cs_tree_node_t * tn_w,int iuma,int ivma,int iwma,int ivimpo,cs_lnum_t nfabor,const cs_zone_t * z,int ialtyb[],double * rcodcl)227 _uialcl_fixed_velocity(cs_tree_node_t *tn_w,
228 int iuma,
229 int ivma,
230 int iwma,
231 int ivimpo,
232 cs_lnum_t nfabor,
233 const cs_zone_t *z,
234 int ialtyb[],
235 double *rcodcl)
236 {
237 /* Get formula */
238 const char *formula = _get_ale_boundary_formula(tn_w, "fixed_velocity");
239
240 if (!formula)
241 bft_error(__FILE__, __LINE__, 0,
242 _("Boundary nature formula is null for %s."),
243 cs_gui_node_get_tag(tn_w, "label"));
244
245 /* Evaluate formula using meg */
246 cs_real_t *bc_vals = cs_meg_boundary_function(z,
247 "mesh_velocity",
248 "fixed_velocity");
249
250 /* Loop over boundary faces */
251 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
252 const cs_lnum_t face_id = z->elt_ids[elt_id];
253
254
255 /* Fill rcodcl */
256 rcodcl[(iuma-1) * nfabor + face_id] = bc_vals[elt_id + 0 * z->n_elts];
257 rcodcl[(ivma-1) * nfabor + face_id] = bc_vals[elt_id + 1 * z->n_elts];
258 rcodcl[(iwma-1) * nfabor + face_id] = bc_vals[elt_id + 2 * z->n_elts];
259
260 ialtyb[face_id] = ivimpo;
261
262 }
263
264 /* Free memory */
265 BFT_FREE(bc_vals);
266
267 }
268
269 /*-----------------------------------------------------------------------------
270 * Return the ale boundary nature of a given boundary condition
271 *
272 * parameters:
273 * tn <-> pointer to tree node for a given mobile_boundary BC
274 *
275 * return:
276 * associated nature
277 *----------------------------------------------------------------------------*/
278
279 static enum ale_boundary_nature
_get_ale_boundary_nature(cs_tree_node_t * tn)280 _get_ale_boundary_nature(cs_tree_node_t *tn)
281 {
282 const char *nat_bndy = cs_tree_node_get_tag(tn, "nature");
283
284 if (cs_gui_strcmp(nat_bndy, "free_surface"))
285 return ale_boundary_nature_free_surface;
286
287 else {
288
289 const char *label = cs_tree_node_get_tag(tn, "label");
290
291 /* get the matching BC node */
292 tn = cs_tree_node_get_child(tn->parent, nat_bndy);
293
294 /* Now searh from siblings */
295 tn = cs_tree_node_get_sibling_with_tag(tn, "label", label);
296
297 /* Finaly get child node ALE */
298 cs_tree_node_t *tn_ale = cs_tree_get_node(tn, "ale/choice");
299 const char *nat_ale = cs_tree_node_get_value_str(tn_ale);
300
301 if (cs_gui_strcmp(nat_ale, "fixed_boundary"))
302 return ale_boundary_nature_fixed_wall;
303 if (cs_gui_strcmp(nat_ale, "sliding_boundary"))
304 return ale_boundary_nature_sliding_wall;
305 else if (cs_gui_strcmp(nat_ale, "internal_coupling"))
306 return ale_boundary_nature_internal_coupling;
307 else if (cs_gui_strcmp(nat_ale, "external_coupling"))
308 return ale_boundary_nature_external_coupling;
309 else if (cs_gui_strcmp(nat_ale, "fixed_velocity"))
310 return ale_boundary_nature_fixed_velocity;
311 else if (cs_gui_strcmp(nat_ale, "fixed_displacement"))
312 return ale_boundary_nature_fixed_displacement;
313 else
314 return ale_boundary_nature_none;
315 }
316 }
317
318 /*-----------------------------------------------------------------------------
319 * Return the ale boundary type of a given boundary condition
320 *
321 * parameters:
322 * tn_bndy <-- pointer to tree node for a given BC
323 *
324 * return:
325 * associated boundary type
326 *----------------------------------------------------------------------------*/
327
328 static cs_boundary_type_t
_get_ale_boundary_type(cs_tree_node_t * tn_bndy)329 _get_ale_boundary_type(cs_tree_node_t *tn_bndy)
330 {
331 const char *nat_bndy = cs_tree_node_get_tag(tn_bndy, "nature");
332
333 if (cs_gui_strcmp(nat_bndy, "free_surface"))
334 return CS_BOUNDARY_ALE_FREE_SURFACE;
335
336 else {
337
338 const char *label = cs_tree_node_get_tag(tn_bndy, "label");
339
340 /* get the matching BC node */
341 cs_tree_node_t *tn = cs_tree_node_get_child(tn_bndy->parent, nat_bndy);
342
343 /* Now search from siblings */
344 tn = cs_tree_node_get_sibling_with_tag(tn, "label", label);
345
346 /* Finally get child node ALE */
347 tn = cs_tree_get_node(tn, "ale/choice");
348 const char *nat = cs_tree_node_get_value_str(tn);
349
350 if (cs_gui_strcmp(nat, "fixed_boundary"))
351 return CS_BOUNDARY_ALE_FIXED;
352 else if (cs_gui_strcmp(nat, "sliding_boundary"))
353 return CS_BOUNDARY_ALE_SLIDING;
354 else if (cs_gui_strcmp(nat, "internal_coupling"))
355 return CS_BOUNDARY_ALE_INTERNAL_COUPLING;
356 else if (cs_gui_strcmp(nat, "external_coupling"))
357 return CS_BOUNDARY_ALE_EXTERNAL_COUPLING;
358 else if (cs_gui_strcmp(nat, "fixed_velocity"))
359 return CS_BOUNDARY_ALE_IMPOSED_VEL;
360 else if (cs_gui_strcmp(nat, "fixed_displacement"))
361 return CS_BOUNDARY_ALE_IMPOSED_DISP;
362 else if (cs_gui_strcmp(nat, "free_surface"))
363 return CS_BOUNDARY_ALE_FREE_SURFACE;
364 else
365 return CS_BOUNDARY_UNDEFINED;
366 }
367 }
368
369 /*-----------------------------------------------------------------------------
370 * Retrieve internal coupling x, y and z values
371 *
372 * parameters:
373 * tn_ic <-- tree node for a given BC's internal coupling definitions
374 * name <-- node name
375 * xyz --> result matrix
376 *----------------------------------------------------------------------------*/
377
378 static void
_get_internal_coupling_xyz_values(cs_tree_node_t * tn_ic,const char * name,double xyz[3])379 _get_internal_coupling_xyz_values(cs_tree_node_t *tn_ic,
380 const char *name,
381 double xyz[3])
382 {
383 cs_tree_node_t *tn = cs_tree_node_get_child(tn_ic, name);
384
385 cs_gui_node_get_child_real(tn, "X", xyz);
386 cs_gui_node_get_child_real(tn, "Y", xyz+1);
387 cs_gui_node_get_child_real(tn, "Z", xyz+2);
388 }
389
390 /*-----------------------------------------------------------------------------
391 * Retrieve data for internal coupling for a specific boundary
392 *
393 * parameters:
394 * label <-- boundary label
395 * xmstru --> Mass matrix
396 * xcstr --> Damping matrix
397 * xkstru --> Stiffness matrix
398 * forstr --> Fluid force matrix
399 * istruc <-- internal coupling boundary index
400 *----------------------------------------------------------------------------*/
401
402 static void
_get_uistr2_data(const char * label,double * xmstru,double * xcstru,double * xkstru,double * forstr,int istruc)403 _get_uistr2_data(const char *label,
404 double *xmstru,
405 double *xcstru,
406 double *xkstru,
407 double *forstr,
408 int istruc)
409 {
410 /* Get mass matrix, damping matrix and stiffness matrix */
411
412 cs_meg_fsi_struct("mass_matrix", label, NULL,
413 &xmstru[istruc * 9]);
414 cs_meg_fsi_struct("damping_matrix", label, NULL,
415 &xcstru[istruc * 9]);
416 cs_meg_fsi_struct("stiffness_matrix", label, NULL,
417 &xkstru[istruc * 9]);
418
419 /* Set variable for fluid force vector */
420 const cs_real_t fluid_f[3] = {forstr[istruc*3],
421 forstr[istruc*3 + 1],
422 forstr[istruc*3 + 2]};
423
424 /* Get fluid force matrix */
425 cs_meg_fsi_struct("fluid_force", label, fluid_f,
426 &forstr[istruc * 3]);
427 }
428
429 /*-----------------------------------------------------------------------------
430 * Return the external coupling dof ("DDL") value
431 *
432 * <boundary_conditions>
433 * <wall label=label_argument">
434 * <ale choice="external_coupling">
435 * <node_name_argument choice="off"/>
436 *
437 * parameters:
438 * label <-- boundary label
439 * node_name <-- Node name: DDLX, DDLY or DDLZ.
440 *----------------------------------------------------------------------------*/
441
442 static int
_get_external_coupling_dof(cs_tree_node_t * tn_ec,const char * name)443 _get_external_coupling_dof(cs_tree_node_t *tn_ec,
444 const char *name)
445 {
446 cs_tree_node_t *tn = cs_tree_node_get_child(tn_ec, name);
447 const char *choice = cs_tree_node_get_child_value_str(tn, "choice");
448
449 int is_on = cs_gui_strcmp(choice, "on");
450
451 return is_on;
452 }
453
454 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
455
456 /*============================================================================
457 * Public Fortran function definitions
458 *============================================================================*/
459
460 /*----------------------------------------------------------------------------
461 * ALE related keywords
462 *
463 * Fortran Interface:
464 *
465 * SUBROUTINE UIALIN
466 * *****************
467 *
468 * nalinf <-> number of sub iteration of initialization
469 * of fluid
470 * nalimx <-> max number of iterations of implicitation of
471 * the displacement of the structures
472 * epalim <-> realtive precision of implicitation of
473 * the displacement of the structures
474 *
475 *----------------------------------------------------------------------------*/
476
CS_PROCF(uialin,UIALIN)477 void CS_PROCF (uialin, UIALIN) (int *nalinf,
478 int *nalimx,
479 double *epalim)
480 {
481 cs_tree_node_t *tn
482 = cs_tree_get_node(cs_glob_tree, "thermophysical_models/ale_method");
483
484 cs_gui_node_get_status_int(tn, &cs_glob_ale);
485
486 if (cs_glob_ale) {
487 cs_gui_node_get_child_int(tn, "fluid_initialization_sub_iterations",
488 nalinf);
489 cs_gui_node_get_child_int(tn, "max_iterations_implicitation",
490 nalimx);
491 cs_gui_node_get_child_real(tn, "implicitation_precision",
492 epalim);
493 }
494
495 #if _XML_DEBUG_
496 bft_printf("==> %s\n", __func__);
497 bft_printf("--cs_glob_ale = %i\n", cs_glob_ale);
498 if (cs_glob_ale > 0) {
499 bft_printf("--nalinf = %i\n", *nalinf);
500 bft_printf("--nalimx = %i\n", *nalimx);
501 bft_printf("--epalim = %g\n", *epalim);
502 }
503 #endif
504 }
505
506 /*----------------------------------------------------------------------------
507 * ALE diffusion type
508 *
509 * Fortran Interface:
510 *
511 * SUBROUTINE UIALVM
512 * *****************
513 *----------------------------------------------------------------------------*/
514
CS_PROCF(uialvm,UIALVM)515 void CS_PROCF (uialvm, UIALVM) ()
516 {
517 cs_tree_node_t *tn
518 = cs_tree_get_node(cs_glob_tree, "thermophysical_models/ale_method");
519
520 int iortvm = _ale_visc_type(tn);
521
522 cs_var_cal_opt_t vcopt;
523 int key_cal_opt_id = cs_field_key_id("var_cal_opt");
524 cs_field_t *f_mesh_u = cs_field_by_name("mesh_velocity");
525 cs_field_get_key_struct(f_mesh_u, key_cal_opt_id, &vcopt);
526
527 if (iortvm == 1) { /* orthotropic viscosity */
528 vcopt.idften = CS_ANISOTROPIC_LEFT_DIFFUSION;
529 } else { /* isotropic viscosity */
530 vcopt.idften = CS_ISOTROPIC_DIFFUSION;
531 }
532 cs_field_set_key_struct(f_mesh_u, key_cal_opt_id, &vcopt);
533
534 #if _XML_DEBUG_
535 bft_printf("==> %s\n", __func__);
536 bft_printf("--iortvm = %i\n", iortvm);
537 }
538 #endif
539 }
540
541 /*-----------------------------------------------------------------------------
542 * uialcl
543 *
544 * parameters:
545 * ialtyb <-- ialtyb
546 * impale <-- uialcl_fixed_displacement
547 * disale <-- See uialcl_fixed_displacement
548 * iuma <-- See uialcl_fixed_velocity
549 * ivma <-- See uialcl_fixed_velocity
550 * iwma <-- See uialcl_fixed_velocity
551 * rcodcl <-- See uialcl_fixed_velocity
552 *----------------------------------------------------------------------------*/
553
554 void CS_PROCF (uialcl, UIALCL) (const int *const ibfixe,
555 const int *const igliss,
556 const int *const ivimpo,
557 const int *const ifresf,
558 int *const ialtyb,
559 int *const impale,
560 cs_real_3_t *disale,
561 const int *const iuma,
562 const int *const ivma,
563 const int *const iwma,
564 double *const rcodcl)
565 {
566 const cs_mesh_t *m = cs_glob_mesh;
567
568 cs_tree_node_t *tn_b0 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
569
570 /* Loop on boundary zones */
571
572 for (cs_tree_node_t *tn_bndy = cs_tree_node_get_child(tn_b0, "boundary");
573 tn_bndy != NULL;
574 tn_bndy = cs_tree_node_get_next_of_name(tn_bndy)) {
575
576 const char *label = cs_tree_node_get_tag(tn_bndy, "label");
577
578 const cs_zone_t *z = cs_boundary_zone_by_name_try(label);
579 if (z == NULL) /* possible in case of old XML file with "dead" nodes */
580 continue;
581
582 cs_lnum_t n_faces = z->n_elts;
583 const cs_lnum_t *faces_list = z->elt_ids;
584
585 /* Get ALE nature and get sibling tn */
586 enum ale_boundary_nature nature = _get_ale_boundary_nature(tn_bndy);
587
588 if (nature == ale_boundary_nature_none)
589 continue;
590
591 /* get the matching BC node */
592 const char *nat_bndy = cs_tree_node_get_tag(tn_bndy, "nature");
593 cs_tree_node_t *tn_bc = cs_tree_node_get_child(tn_bndy->parent, nat_bndy);
594 tn_bc = cs_tree_node_get_sibling_with_tag(tn_bc, "label", label);
595
596 if (nature == ale_boundary_nature_fixed_wall) {
597 for (cs_lnum_t ifac = 0; ifac < n_faces; ifac++) {
598 cs_lnum_t ifbr = faces_list[ifac];
599 ialtyb[ifbr] = *ibfixe;
600 }
601 }
602 else if (nature == ale_boundary_nature_sliding_wall) {
603 for (cs_lnum_t ifac = 0; ifac < n_faces; ifac++) {
604 cs_lnum_t ifbr = faces_list[ifac];
605 ialtyb[ifbr] = *igliss;
606 }
607 }
608 else if (nature == ale_boundary_nature_free_surface) {
609 for (cs_lnum_t ifac = 0; ifac < n_faces; ifac++) {
610 cs_lnum_t ifbr = faces_list[ifac];
611 ialtyb[ifbr] = *ifresf;
612 }
613 }
614 else if (nature == ale_boundary_nature_fixed_displacement) {
615 _uialcl_fixed_displacement(tn_bc,
616 z,
617 impale,
618 disale);
619 }
620 else if (nature == ale_boundary_nature_fixed_velocity) {
621 _uialcl_fixed_velocity(tn_bc, *iuma, *ivma, *iwma, *ivimpo,
622 m->n_b_faces,
623 z,
624 ialtyb, rcodcl);
625 }
626 }
627 }
628
629 /*-----------------------------------------------------------------------------
630 * Retrieve data for internal coupling. Called once at initialization
631 *
632 * parameters:
633 * idfstr --> Structure definition
634 * mbstru <-- number of previous structures (-999 or by restart)
635 * aexxst --> Displacement prediction alpha
636 * bexxst --> Displacement prediction beta
637 * cfopre --> Stress prediction alpha
638 * ihistr --> Monitor point synchronisation
639 * xstr0 <-> Values of the initial displacement
640 * xstreq <-> Values of the equilibrium displacement
641 * vstr0 <-> Values of the initial velocity
642 *----------------------------------------------------------------------------*/
643
644 void CS_PROCF (uistr1, UISTR1) (cs_lnum_t *idfstr,
645 const int *mbstru,
646 double *aexxst,
647 double *bexxst,
648 double *cfopre,
649 int *ihistr,
650 double *xstr0,
651 double *xstreq,
652 double *vstr0)
653 {
654 int istruct = 0;
655
656 cs_tree_node_t *tn0
657 = cs_tree_get_node(cs_glob_tree, "thermophysical_models/ale_method");
658
659 /* Get advanced data */
660 cs_gui_node_get_child_real(tn0, "displacement_prediction_alpha", aexxst);
661 cs_gui_node_get_child_real(tn0, "displacement_prediction_beta", bexxst);
662 cs_gui_node_get_child_real(tn0, "stress_prediction_alpha", cfopre);
663 cs_gui_node_get_child_status_int(tn0, "monitor_point_synchronisation", ihistr);
664
665 cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
666
667 cs_tree_node_t *tn_b0 = cs_tree_node_get_child(tn, "boundary");
668 cs_tree_node_t *tn_w0 = cs_tree_node_get_child(tn, "boundary");//FIXME
669
670 /* At each time-step, loop on boundary faces */
671
672 int izone = 0;
673
674 for (tn = tn_b0;
675 tn != NULL;
676 tn = cs_tree_node_get_next_of_name(tn), izone++) {
677
678 const char *label = cs_tree_node_get_tag(tn, "label");
679
680 cs_tree_node_t *tn_w
681 = cs_tree_node_get_sibling_with_tag(tn_w0, "label", label);
682
683 /* Keep only internal coupling */
684 if ( _get_ale_boundary_nature(tn_w)
685 == ale_boundary_nature_internal_coupling) {
686
687 if (istruct+1 > *mbstru) { /* Do not overwrite restart data */
688 /* Read initial_displacement, equilibrium_displacement
689 and initial_velocity */
690 cs_tree_node_t *tn_ic = cs_tree_get_node(tn_w, "ale");
691 tn_ic = cs_tree_node_get_sibling_with_tag(tn_ic,
692 "choice",
693 "internal_coupling");
694 _get_internal_coupling_xyz_values(tn_ic, "initial_displacement",
695 &xstr0[3 * istruct]);
696 _get_internal_coupling_xyz_values(tn_ic, "equilibrium_displacement",
697 &xstreq[3 * istruct]);
698 _get_internal_coupling_xyz_values(tn_ic, "initial_velocity",
699 &vstr0[3 * istruct]);
700 }
701
702 const cs_zone_t *z = cs_boundary_zone_by_name_try(label);
703 if (z == NULL) /* possible in case of old XML file with "dead" nodes */
704 continue;
705
706 cs_lnum_t n_faces = z->n_elts;
707 const cs_lnum_t *faces_list = z->elt_ids;
708
709 /* Set idfstr to positive index starting at 1 */
710 for (cs_lnum_t ifac = 0; ifac < n_faces; ifac++) {
711 cs_lnum_t ifbr = faces_list[ifac];
712 idfstr[ifbr] = istruct + 1;
713 }
714 ++istruct;
715 }
716 }
717 }
718
719 /*-----------------------------------------------------------------------------
720 * Retrieve data for internal coupling. Called at each step
721 *
722 * parameters:
723 * xmstru --> Mass matrix
724 * xcstr --> Damping matrix
725 * xkstru --> Stiffness matrix
726 * forstr --> Fluid force matrix
727 *----------------------------------------------------------------------------*/
728
729 void CS_PROCF (uistr2, UISTR2) (double *const xmstru,
730 double *const xcstru,
731 double *const xkstru,
732 double *const forstr)
733 {
734 int istru = 0;
735
736 cs_tree_node_t *tn_b0 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
737
738 /* At each time-step, loop on boundary faces */
739
740 for (cs_tree_node_t *tn_bndy = cs_tree_node_get_child(tn_b0, "boundary");
741 tn_bndy != NULL;
742 tn_bndy = cs_tree_node_get_next_of_name(tn_bndy)) {
743
744 const char *label = cs_tree_node_get_tag(tn_bndy, "label");
745
746 enum ale_boundary_nature nature =_get_ale_boundary_nature(tn_bndy);
747
748 /* Keep only internal coupling */
749 if (nature == ale_boundary_nature_internal_coupling) {
750
751 /* get the matching BC node */
752 const char *nat_bndy = cs_tree_node_get_tag(tn_bndy, "nature");
753 cs_tree_node_t *tn_bc = cs_tree_node_get_child(tn_bndy->parent, nat_bndy);
754 tn_bc = cs_tree_node_get_sibling_with_tag(tn_bc, "label", label);
755
756 cs_tree_node_t *tn_ic = cs_tree_get_node(tn_bc, "ale");
757 tn_ic = cs_tree_node_get_sibling_with_tag(tn_ic,
758 "choice",
759 "internal_coupling");
760
761 _get_uistr2_data(label,
762 xmstru,
763 xcstru,
764 xkstru,
765 forstr,
766 istru);
767 ++istru;
768 }
769 }
770 }
771
772 /*-----------------------------------------------------------------------------
773 * Retrieve data for external coupling
774 *
775 * parameters:
776 * idfstr <-- Structure definition
777 * asddlf --> Block of the DDL forces
778 *----------------------------------------------------------------------------*/
779
780 void
781 CS_PROCF(uiaste, UIASTE)(int *idfstr,
782 int *asddlf)
783 {
784 int istruct = 0;
785
786 cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
787
788 cs_tree_node_t *tn_b0 = cs_tree_node_get_child(tn, "boundary");
789 cs_tree_node_t *tn_w0 = cs_tree_node_get_child(tn, "boundary");//FIXME
790
791 /* Loop on boundary faces */
792
793 for (tn = tn_b0; tn != NULL; tn = cs_tree_node_get_next_of_name(tn)) {
794
795 const char *label = cs_tree_node_get_tag(tn, "label");
796
797 cs_tree_node_t *tn_w
798 = cs_tree_node_get_sibling_with_tag(tn_w0, "label", label);
799
800 enum ale_boundary_nature nature =_get_ale_boundary_nature(tn_w);
801
802 if (nature == ale_boundary_nature_external_coupling) {
803
804 const cs_zone_t *z = cs_boundary_zone_by_name_try(label);
805 if (z == NULL) /* possible in case of old XML file with "dead" nodes */
806 continue;
807
808 cs_lnum_t n_faces = z->n_elts;
809 const cs_lnum_t *faces_list = z->elt_ids;
810
811 cs_tree_node_t *tn_ec = cs_tree_get_node(tn_w, "ale");
812 tn_ec = cs_tree_node_get_sibling_with_tag(tn_ec,
813 "choice",
814 "external_coupling");
815
816 /* Get DDLX, DDLY and DDLZ values */
817 asddlf[istruct*3 + 0] = _get_external_coupling_dof(tn_ec, "DDLX") ? 0 : 1;
818 asddlf[istruct*3 + 1] = _get_external_coupling_dof(tn_ec, "DDLY") ? 0 : 1;
819 asddlf[istruct*3 + 2] = _get_external_coupling_dof(tn_ec, "DDLZ") ? 0 : 1;
820
821 /* Set idfstr with negative value starting from -1 */
822 for (cs_lnum_t ifac = 0; ifac < n_faces; ifac++) {
823 cs_lnum_t ifbr = faces_list[ifac];
824 idfstr[ifbr] = -istruct - 1;
825 }
826 istruct++;
827 }
828
829 }
830 }
831
832 /*============================================================================
833 * Public function definitions
834 *============================================================================*/
835
836 /*-----------------------------------------------------------------------------
837 * Return the viscosity's type of ALE method
838 *
839 * parameters:
840 * type --> type of viscosity's type
841 *----------------------------------------------------------------------------*/
842
843 void
844 cs_gui_get_ale_viscosity_type(int *type)
845 {
846 cs_tree_node_t *tn
847 = cs_tree_get_node(cs_glob_tree,
848 "thermophysical_models/ale_method");
849
850 *type = _ale_visc_type(tn);
851 }
852
853 /*----------------------------------------------------------------------------
854 * Mesh viscosity setting.
855 *----------------------------------------------------------------------------*/
856
857 void
858 cs_gui_mesh_viscosity(void)
859 {
860 cs_tree_node_t *tn0
861 = cs_tree_get_node(cs_glob_tree, "thermophysical_models/ale_method");
862
863 /* Get formula */
864 const char *mvisc_expr = cs_tree_node_get_child_value_str(tn0, "formula");
865
866 if (mvisc_expr == NULL)
867 return;
868
869 /* Compute mesh viscosity using the MEG function.
870 * The function accepts a field of arbitrary dimension, hence no need to
871 * check here. The formula is generated by the code.
872 */
873 const cs_zone_t *z_all = cs_volume_zone_by_name("all_cells");
874 cs_field_t *fmeg[1] = {CS_F_(vism)};
875 cs_meg_volume_function(z_all, fmeg);
876
877 }
878
879 /*----------------------------------------------------------------------------*/
880 /*!
881 * \brief Translate the user settings for the domain boundaries into a
882 * structure storing the ALE boundaries (New mechanism used in CDO)
883 *
884 * \param[in, out] domain pointer to a \ref cs_domain_t structure
885 */
886 /*----------------------------------------------------------------------------*/
887
888 void
889 cs_gui_mobile_mesh_get_boundaries(cs_domain_t *domain)
890 {
891 assert(domain != NULL);
892
893 cs_tree_node_t *tn_b0 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
894
895 /* Loop on boundary zones */
896
897 for (cs_tree_node_t *tn_bndy = cs_tree_node_get_child(tn_b0, "boundary");
898 tn_bndy != NULL;
899 tn_bndy = cs_tree_node_get_next_of_name(tn_bndy)) {
900
901 const char *label = cs_tree_node_get_tag(tn_bndy, "label");
902 const cs_zone_t *z = cs_boundary_zone_by_name_try(label);
903 if (z == NULL) /* possible in case of old XML file with "dead" nodes */
904 continue;
905
906 cs_boundary_type_t ale_bdy = _get_ale_boundary_type(tn_bndy);
907
908 if (ale_bdy == CS_BOUNDARY_UNDEFINED)
909 continue;
910
911 cs_boundary_add(domain->ale_boundaries,
912 ale_bdy,
913 z->name);
914
915 } /* Loop on mobile_boundary zones */
916
917 /* TODO */
918 /* else if (nature == ale_boundary_nature_fixed_displacement) { */
919 /* _uialcl_fixed_displacement(tn_bndy, z, */
920 /* impale, disale); */
921 /* } */
922 /* else if (nature == ale_boundary_nature_fixed_velocity) { */
923 /* _uialcl_fixed_velocity(tn_bndy, *iuma, *ivma, *iwma, *ivimpo, */
924 /* m->n_b_faces, z, */
925 /* ialtyb, rcodcl); */
926 /* } */
927
928 }
929
930 /*----------------------------------------------------------------------------*/
931 /*!
932 * \brief Return the fixed velocity for a boundary
933 *
934 * \param[in] label boundary condition label
935 *
936 * \return a pointer to an array of cs_real_t values
937 */
938 /*----------------------------------------------------------------------------*/
939
940 cs_real_t *
941 cs_gui_mobile_mesh_get_fixed_velocity(const char *label)
942 {
943 cs_tree_node_t *tn_b0 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
944
945 /* Loop on boundary zones */
946
947 for (cs_tree_node_t *tn_bndy = cs_tree_node_get_child(tn_b0, "boundary");
948 tn_bndy != NULL;
949 tn_bndy = cs_tree_node_get_next_of_name(tn_bndy)) {
950
951 const char *nat_bndy = cs_tree_node_get_tag(tn_bndy, "nature");
952 const char *label_bndy = cs_tree_node_get_tag(tn_bndy, "label");
953
954 /* get the matching BC node */
955 cs_tree_node_t *tn = cs_tree_node_get_child(tn_bndy->parent, nat_bndy);
956
957 /* Now searh from siblings */
958 tn = cs_tree_node_get_sibling_with_tag(tn, "label", label_bndy);
959
960 if (strcmp(label_bndy, label) == 0) {
961
962 /* Get formula */
963 const char *formula = _get_ale_boundary_formula(tn, "fixed_velocity");
964
965 if (!formula)
966 bft_error(__FILE__, __LINE__, 0,
967 _("Boundary nature formula is null for %s."),
968 cs_gui_node_get_tag(tn, "label"));
969
970 const cs_zone_t *bz = cs_boundary_zone_by_name(label);
971
972 /* Evaluate formula using meg */
973 return cs_meg_boundary_function(bz,
974 "mesh_velocity",
975 "fixed_velocity");
976
977 }
978 }
979
980 return NULL; /* avoid a compilation warning */
981 }
982
983 /*----------------------------------------------------------------------------*/
984
985 END_C_DECLS
986