IBComputeNormalGradient.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <arrayfunctions.h>
#include <immersedboundaries.h>
#include <mpivars.h>
#include <hypar.h>

int IBComputeNormalGradient(void*               s,       
                            void*               m,       
                            const double* const var,     
                            int                 nvars,   
                            double** const      grad_var 
                          )
{
  HyPar             *solver  = (HyPar*)          s;
  MPIVariables      *mpi     = (MPIVariables*)   m; 
  ImmersedBoundary  *IB      = (ImmersedBoundary*) solver->ib;

  if (!solver->flag_ib) return(0);

  if ((*grad_var) != NULL) {
    fprintf(stderr,"Error in IBComputeNormalGradient()\n");
    fprintf(stderr," grad_var is not NULL on rank %d\n",
            mpi->rank );
    return 1;
  }

  int nfacets_local = IB->nfacets_local;
  FacetMap *fmap = IB->fmap;

  if (nfacets_local > 0) {
    (*grad_var) = (double*) calloc (nvars*nfacets_local, sizeof(double));

    for (int n = 0; n < nfacets_local; n++) {
  
      double *alpha;
      int    *nodes, j, k;
  
      double v_c[nvars];
      alpha = &(fmap[n].interp_coeffs[0]);
      nodes = &(fmap[n].interp_nodes[0]);
      _ArraySetValue_(v_c,nvars,0.0);
      for (j=0; j<_IB_NNODES_; j++) {
        for (k=0; k<nvars; k++) {
          v_c[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
        }
      }
  
      double v_ns[nvars];
      alpha = &(fmap[n].interp_coeffs_ns[0]);
      nodes = &(fmap[n].interp_nodes_ns[0]);
      _ArraySetValue_(v_ns,nvars,0.0);
      for (j=0; j<_IB_NNODES_; j++) {
        for (k=0; k<nvars; k++) {
          v_ns[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
        }
      }
  
     double nx = fmap[n].facet->nx;
     double ny = fmap[n].facet->ny;
     double nz = fmap[n].facet->nz;
  
      for (k=0; k<nvars; k++) {
        double dv_dx = (v_ns[k] - v_c[k]) / fmap[n].dx;
        double dv_dy = (v_ns[k] - v_c[k]) / fmap[n].dy;
        double dv_dz = (v_ns[k] - v_c[k]) / fmap[n].dz;
        (*grad_var)[n*nvars+k] = dv_dx*nx + dv_dy*ny + dv_dz*nz;
      }
  
    }

  }

  return(0);
}