Removing mesh volume elements - Evolving (??) mesh

[Pisolino]

Hi all,
i'm wondering if it is possible, in order to perform shape evolution, in CS to remove some elements (tetra) at the boundary from the mesh and to project the same boundary condition to the new-generated boundary faces ? This would be done every N simulations, then the simulation re-starts again.
Thanks in advance

Andrea

[Yvan Fournier]

Hello,

you could easily use the postprocessing features to export the "mesh - removed cells" to MED format, and run a new computation from that, but we currently do not have true restart (with previous values) options for that (it would require restarting with interpolation, which is on the wish list, but not implemented yet).

Regards,

Yvan

[Pisolino]

Thank you for the reply, for my purposes it's ok to restart the simulation from initial values (maybe it's a little slow but it is ok), i'm very happy that it seems possible to perform such simulation. I'm taking a look at cs_user_postprocessing functions. My objective is to automatically:
- after N iteration in steady state check every face of a given boundary
- remove the cells associated to some boundary faces that do not satisfy a given criteria (value lower than..)
- transfer the boundary property to the new boundary faces
- save it as a mesh
- start another N iterations simulation using the new mesh

I'm wondering if there is any example about transferring the boundary face property to the new generated boundary, maybe i can just set the old boundary ID of the removed faces to the new created ?

Again, is there any example about automatically restarting (from initial values, not previous values) using the new mesh ?

Thanks for every support you can give me

Andrea

[Yvan Fournier]

Hello,

For postprocessing output of a volume mesh, check cs_post_define_volume_mesh_by_func in the Doxygen examples (adapt a function which selects all cells except the ones you want to remove).
Fo this use, you need a time-varying connectivity (so the extracted mesh is defined at the current and not initial time step) and MED output format (to keep volume and boundary groups).

To transfer boundary face groups, things are a bit more trick, as the postprocessing output does not include hooks to add those.

But before removing those cells, in a prior stage (such as cs_user_extra_operations.f90, you can add groups to interior faces belonging to those cells (using the cs_mesh_group_i_faces_add, which is new in 5.0-beta). There are no examples yet for this, though the Doxygen documentation is up-to-date.

Regards,

Yvan

[Pisolino]

Hi,
again thank you for the support. Just to be sure i understood everything and as a future starting point for other users:
in the example cs_user_postprocess.c file :

the example function (_t_lt_21_select) select the desired cells with temperature == 21.

then a writer is defined, i'll define the mesh as med, about time_dep i suppose i have to choose FVM_WRITER_FIXED_MESH since i use a fixed mesh during the calculation.

then cs_post_define_volume_mesh_by_func is called and i have to set time-varying to true so at the desired timestep the mesh will be extracted with the calculated (and not initial) values.

Question #1 about the 5.0 cs_mesh_group_i_faces_add : am i forced to use the all new beta version or i can just import the function in my 4.0 environment ?

Question #2 from what i understood the modified mesh is saved after N timestep, but then it is not so clear to me (sorry) how to automatically start another calculation with the new mesh. Since it is an optimization loop every simulation should count something like 20/30 restarts and it is not feasible to perform it manually.

Thanks again

[Yvan Fournier]

Hello,

You actually need to set the writer to FVM_WRITER_TRANSIENT_CONNECT, otherwise it will force the mesh to "fixed" and force the output to the initial connectivity.

Actually, this might be an issue, because the current MED options do not allow transient connectivity. If you only do one output before restarting, removing the check (simply replacing FVM_WRITER_FIXED_MESH with FVM_WRITER_TRANSIENT_CONNECT around line 138-140 in src/fvm/fvm_writer.c and reinstalling) should probably work.

Adding the new function from 5.0-beta to 4.0 is quite complex, as the code change was significant (with quite a bit of refactorization).

To use 4.0, you could also simply edit cs_glob_mesh->i_face_family[] and cs_glob_mesh->b_face_family[], copying the value from the removed boundary face (b_face_family[face_id]) to the neighboring interior faces (i_face_family[face_id]). To do this, the simplest solution is probably to use a temporary integer array allocated at the is of cs_glob_mesh->n_cells_with_ghosts, initialize all values to 0, then loop on boundary faces to mark adjacent cells with the family number of boundary faces where those are marked for removal. Then do a second loop on interior faces and copy the value on cells to adjacent interior faces.

Regards,

Yvan

[Pisolino]

Sorry to bother you again, after a full day in solving this i'm almost losing my mind. I'm continuing to read from fortran and C subroutines and i don't have all the commands to perform something with a language or another one. Despite the logic of what i would do is quite clear to me i don't understand how to perform it, i'll try to be as clear as i can:

-1 select all the faces belonging to a group called "boun1". In f90 this is done with:
subroutine getfbr (fstr, facnb, faces)
Build the list of boundary faces matching a criteria string
in C i can call all of them with cs_glob_mesh->b_face_cells but then i don't know how to check if they belong to "boun1" ? it is the value of the field family (cs_glob_mesh->b_face_family) ?

-2 scroll the boundary faces list and for every face belonging to "boun1" check the property value (yplus), depending on the condition perform an action. Importing the value is described both for f90 and C. For C:
cs_real_t *bpro_yplus = cs_field_by_name("yplus")->val

-3 if condition is reached: find the cell c-id that owns the boundary face b-id
don't know how to do, neither in f90 or C

-4 once the cell is found, select all the interior faces of the cell and set em belonging to "boun1"
don't know how to do, neither in f90 or C

Any hint will be very appreciated

Andrea

[Yvan Fournier]

Hello,

You may use something like this:

Code: Select all

const cs_mesh_t *mesh = cs_glob_mesh;
cs_lnum_t n_elts;
cs_lnum_t *elt_ids;

const int n_cells_ext = mesh->n_cells_with_ghosts;

/* Initialize cell marker */
int *cell_marker;
BFT_MALLOC(cell_marker, n_cells_ext, int);
for (cs_lnum_t i = 0; i < n_cells; i++)
  cell_marker[i] = 0;

/* Select boundary faces */
BFT_MALLOC(elt_ids, n_b_faces, cs_lnum_t);
cs_selector_get_b_face_list("boun1", &n_elts, elt_ids);

/* Mark adjacent cells (with boundary face family) */
for (cs_lnum_t i = 0; i < n_elts; i++) {
  cs_lnum_t face_id = elt_list[i];
  cell_marker[mesh->b_face_cells[face_id] = mesh->b_face_family[face_id];
}

/* Now propagate to interior faces */
for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++) {
  for (cs_lnum_t j = 0; j < 2; j++) {
    cs_lnum_t cell_id = mesh->i_face_cells[i][j];
    if (cell_marker[cell_id] != 0)
       mesh->i_face_family[face_id] = cell_marker[cell_id];
  }
}

...

BFT_FREE(cell_marker);
BFT_FREE(elt_ids);

The above handles the propagation of boundary face groups to interior faces.

You can actually simultaneously do this and then use the cell marker (before freeing it) in a selection function for cells (see the Doxygen example for cs_post_define_volume_mesh_by_func, adapt it by selecting all cells for which cell_marker is still 0).

Regards,

Yvan

[Pisolino]

Hi,
thank you again for the big hint. i've started from your function

Code: Select all

static void
bounfunc(void        *input,
                cs_lnum_t   *n_cells,
                cs_lnum_t  **cell_ids)
{
 cs_lnum_t i, face_id;

  cs_lnum_t _n_cells = 0;
  cs_lnum_t *_cell_ids = NULL;
  cs_lnum_t n_b_faces = 0;
  //cs_lnum_t *b_face_ids = NULL;
  
 /* import the mesh */
 const cs_mesh_t *mesh = cs_glob_mesh;
 /* define variable and pointers*/
cs_lnum_t n_elts;
cs_lnum_t *elt_ids;
/*import all cells, also with halo - ghost */
const int n_cells_ext = mesh->n_cells_with_ghosts;

/* Initialize cell marker */
int *cell_marker;
/* malloc in cell marker the size of the cells, ghost included*/
BFT_MALLOC(cell_marker, n_cells_ext, int);

/* mark all cells as 0 */
for (cs_lnum_t i = 0; i < n_cells; i++)
  cell_marker[i] = 0;

/* Select boundary faces belonging to "heat" and save 
*  their number and id  */
BFT_MALLOC(elt_ids, n_b_faces, cs_lnum_t);
cs_selector_get_b_face_list("heat", &n_elts, elt_ids);

/* Mark adjacent cells (with boundary face family) */
/* scroll all the selected element list   */
for (cs_lnum_t i = 0; i < n_elts; i++) {
  cs_lnum_t face_id = elt_ids[i];
/* marking cells cell_marker[cell_id]= Boundary face family  */
  cell_marker[mesh->b_face_cells[face_id]] = mesh->b_face_family[face_id];
}

/* Now propagate to interior faces */
for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++) { //for all the i_faces
  for (cs_lnum_t j = 0; j < 2; j++) {
    cs_lnum_t cell_id = mesh->i_face_cells[i][j]; //pick the cell id of i_face
    if (cell_marker[cell_id] != 0) //if it is the cell to be removed
       mesh->i_face_family[face_id] = cell_marker[cell_id]; //
  }
}

/* now pick only the all - undesired cells to a new mesh*/
BFT_MALLOC(_cell_ids, mesh->n_cells, cs_lnum_t); /* Allocate selection list */
    for (i = 0; i < mesh->n_cells; i++) {
      if (cell_marker[i] == 0) {
        _cell_ids[_n_cells] = i;
        _n_cells += 1;
      }
    }

    BFT_REALLOC(_cell_ids, _n_cells, cs_lnum_t); /* Adjust size (good practice,
                                                    but not required) */
 
 /* Set return values */
  *n_cells = _n_cells;
  *cell_ids = _cell_ids;
 
BFT_FREE(cell_marker);
BFT_FREE(elt_ids);
} 

and i've set an ensight writer ( id = 1) just to test the function behaviour :

Code: Select all

cs_post_define_writer(1,                            /* writer_id */
                          "user",                       /* writer name */
                          "postprocessing",             /* directory name */
                          "EnSight Gold",               /* format name */
                          "",
                          FVM_WRITER_TRANSIENT_CONNECT,
                          1,
                          50,
                          -1.0);

and then this is the final part:

Code: Select all

void
cs_user_postprocess_meshes(void)
{
      cs_post_define_volume_mesh_by_func( -1,               /* mesh id */
                                       "opti_cfd",
                                       bounfunc,
                                       NULL,            /*  */
                                       true,            /* time varying */
                                       true,           /* add_groups */
                                       true,           /* auto_variables */
                                       1,
                                       1);
}

No matter what, if i turn on the function cs_post_define_volume_mesh i have a SIGSEGV error

Code: Select all

 Warning:
 --------
    Mesh quality issue has been detected

    The mesh should be re-considered using the listed criteria.

    The calculation will run but the solution quality may be degraded...

 Computing geometric quantities (0.094 s)
SIGSEGV signal (forbidden memory area access) intercepted!

i also tryied with the standard function "_t_lt_21_select" of the tutorial, the same error occurs, the only way to remove it is to turn off the function cs_post_define_volume_mesh. I've already changed the mesh id, switch another writer (-1 the standard one ) with no changes at all in the final error.
Do you see any big mistake ?

Andrea

[Yvan Fournier]

Hello,

I did not notice anything obvious, but to be safe:

- test with a user mesh and user writer (1 and A for example), as you already did in some variants
- test this on a debug build

With a debug build, you can either run under a debugger, or at least get a backtrace which, combined with the addr2line utility, will tell you where the SIGSEGV happens. As the function is not so trivial, using tools to see where the error happens and working from there should help.

Regards,

Yvan