sandbox/huet/src/lagrangian_caps/volume-conservation-ft.h

    Enforcing volume conservation

    In this file, we follow the method of Sigüenza et al. to enforce exact volume conservation of the capsule. The details can be found in appendix A of [1].

    We will need to find the smallest real root of a cubic equation, which we implement in a separate file.

    #include "smallest_root_cubic.h"

    The function below computes the quantity \bm{\alpha_m} as specified in Eq.,(41) of 1.

    trace
coord compute_alpha_m(lagMesh* mesh, int m) {
  coord alpha = {0., 0., 0.};

  for(int i=0; i<mesh->nodes[m].nb_triangles; i++) {
    int tid = mesh->nodes[m].triangle_ids[i];
    int local_id = 0;
    while (mesh->triangles[tid].node_ids[local_id] != m && local_id < 3) 
      local_id++;
    int nids[2];  // `nids` for "neighbor ids"
    nids[0] = mesh->triangles[tid].node_ids[(local_id + 1)%3];
    nids[1] = mesh->triangles[tid].node_ids[(local_id + 2)%3];

    foreach_dimension()
      alpha.x += -periodic_friendly_cross_product_x(
        mesh->nodes[nids[0]].pos, mesh->nodes[nids[1]].pos, 
        mesh->centroid)/12;
  }

  return alpha;
}
    The function below solves an optimization problem to enforce the conservation of volume of a capsule. The cost function to minimize is the following:

    where \bm{X} = \left[\bm{X}_1 \cdots \bm{X}_n \right] is a tensor storing the coordinates \bm{X}_i of the n capsule nodes, \bm{\delta X} = \left[\bm{\delta X}_1 \cdots \bm{\delta X}_n \right] is a tensor storing their displacements, V(\bm{X}) is the capsule volume in configuration \bm{X}, V_0 is the initial volume and \Lambda is a Lagrange multiplier. Sigüenza et al. show that \Lambda is a (real) root of a third order polynomial which coefficients are given by Eqs.,(43-46) in 1. An assumption for their derivation is \| \bm{\delta X} \| \ll \| \bm{X} \|, so the solution to the optimization will not quite be exact, but this approximation is asymptotically true as \Delta t tends to zero.

    trace
void enforce_optimal_volume_conservation(lagMesh* mesh) {

    First, the \bm{\alpha_m} are computed and stored in an array of size the number of Lagrangian nodes.

      coord* alpha = malloc(mesh->nln*sizeof(coord));
  for(int m=0; m<mesh->nln; m++)
    alpha[m] = compute_alpha_m(mesh, m);
    Once all the \bm{\alpha_m} are known, the coefficients of the polynomial in \Lambda can be computed: with and

    where \star_{i,1} is the quantity \star at vertex 1 of triangle i, V is the capsule volume and V_0 is the initial (conserved) volume. In practice, we normalize the polynomial so that its largest coefficient is 1.

      double coeff_polynomial[4];
  for(int k=0; k<4; k++) coeff_polynomial[k] = 0;
  for(int j=0; j<mesh->nlt; j++) {
    int tn[3]; // `tn` for "triangle nodes"
    for(int k=0; k<3; k++) 
      tn[k] = mesh->triangles[j].node_ids[k];
    for(int k=0; k<3; k++) {
      coord cp0; coord cp2; // `cp` for "cross-product"
      foreach_dimension() {
        cp0.x = alpha[tn[(k+1)%3]].y*alpha[tn[(k+2)%3]].z
          - alpha[tn[(k+1)%3]].z*alpha[tn[(k+2)%3]].y;
        cp2.x = periodic_friendly_cross_product_x(
          mesh->nodes[tn[(k+1)%3]].pos, mesh->nodes[tn[(k+2)%3]].pos,
          mesh->centroid);
      }
      coeff_polynomial[3] += cdot(alpha[tn[k]], cp0);
      coeff_polynomial[2] += periodic_friendly_dot_product(
      mesh->nodes[tn[k]].pos, cp0, mesh->centroid);
      coeff_polynomial[1] += cdot(alpha[tn[k]], cp2);
    }
  }
  coeff_polynomial[3] /= 18;
  coeff_polynomial[2] /= 6;
  coeff_polynomial[1] /= 6;
  coeff_polynomial[0] = mesh->volume - mesh->initial_volume;
  double normalize_factor = 
    max(max(fabs(coeff_polynomial[3]), coeff_polynomial[2]), 
    max(fabs(coeff_polynomial[1]), coeff_polynomial[0]));
  coeff_polynomial[3] /= normalize_factor;
  coeff_polynomial[2] /= normalize_factor;
  coeff_polynomial[1] /= normalize_factor;
  coeff_polynomial[0] /= normalize_factor;

    Once the coefficients are found, we find its smallest real root (in absolute value) thanks to a separate routine.

      double lambda = find_smallest_real_root(coeff_polynomial);
    Finally, the correction applied to each node i is 
      for(int i=0; i<mesh->nln; i++) {
    foreach_dimension() 
      mesh->nodes[i].pos.x += lambda*alpha[i].x;
  }
  free(alpha);
  correct_lag_pos(mesh);
}

    References

    [siguenza2016validation]

    Julien Sigüenza, Simon Mendez, Dominique Ambard, Frédéric Dubois, Franck Jourdan, Rémy Mozul, and Franck Nicoud. Validation of an immersed thick boundary method for simulating fluid–structure interactions of deformable membranes. Journal of Computational Physics, 322:723–746, 2016.