/****************************************************************************** UDF for defining the anisotropic conductivity matrix for a cylindrical shell This UDF is from Ansys Help document and Dr Du Mengjie modified and tested! ******************************************************************************/
#include "udf.h"
#define CELL_NUM 2    //number of cellsreal origin[CELL_NUM][3];  //coordinates of every single cell axisstatic const real axis[3]  = {0, 1, 0};int id[CELL_NUM] = {144,147};  //zone id of cells, a function of assignment can be adapted if there are too many cellsstatic const real cond[3] = {1, 100, 100};  //radial, tangential and axial directions
DEFINE_ON_DEMAND(on_demand_func){  int i;  cell_t c;  Domain *d = Get_Domain(1);  Thread *t;  real xmin,xmax,ymin,ymax,zmin,zmax;  real xcent[ND_ND];  for ( i = 0; i < CELL_NUM; i++)   {    xmin = 1e30;xmax = -1e30;    ymin = 1e30;ymax = -1e30;    zmin = 1e30;zmax = -1e30;    t = Lookup_Thread(d,id[i]); //obtain cell[i] thread #if !RP_HOST  begin_c_loop(c,t) //obtain cell[i] domain x/y/z min/max coordinates  {    C_CENTROID(xcent,c,t);    if ( xcent[0] < xmin ) xmin = xcent[0];    if ( xcent[0] > xmax ) xmax = xcent[0];    if ( xcent[1] < ymin ) ymin = xcent[1];    if ( xcent[1] > ymax ) ymax = xcent[1];    if ( xcent[2] < zmin ) zmin = xcent[2];    if ( xcent[2] > zmax ) zmax = xcent[2];  }  end_c_loop(c,t)#if RP_NODE  //parallel reduction  xmin = PRF_GRLOW1(xmin);  xmax = PRF_GRHIGH1(xmax);  ymin = PRF_GRLOW1(ymin);  ymax = PRF_GRHIGH1(ymax);  zmin = PRF_GRLOW1(zmin);  zmax = PRF_GRHIGH1(zmax);#endif#endif  node_to_host_real_6(xmin,xmax,ymin,ymax,zmin,zmax);  //obtain cell[i] center coordinates value and assign them to global array origin  origin[i][0] = (xmax + xmin) / 2;  origin[i][1] = (ymax + ymin) / 2;  origin[i][2] = (zmax + zmin) / 2;  }#if !RP_NODE  for ( i = 0 ; i < CELL_NUM; i++)    Message("In zone %d, center is %lf, %lf, %lf.\n",id[i], origin[i][0],origin[i][1], origin[i][2]);#endif}
DEFINE_ANISOTROPIC_CONDUCTIVITY(cyl_ortho_cond,c,t,dmatrix){   real x[3][3]; /* principal direction matrix for cell in cartesian coords. */   real xcent[ND_ND];   real R;   int i;   for ( i = 0; i < CELL_NUM; i++)     if ( id[i] == THREAD_ID(t))       break;   C_CENTROID(xcent,c,t);   //Message("In zone %d, center is %lf, %lf, %lf.\n",id[i], origin[i][0],origin[i][1], origin[i][2]);   NV_VV(x[0],=,xcent,-,origin[i]);  #if RP_3D     NV_V(x[2],=,axis);  #endif  #if RP_3D     R = NV_DOT(x[0],x[2]);     NV_VS(x[0],-=,x[2],*,R);  #endif   R = NV_MAG(x[0]);   if (R > 0.0)      NV_S(x[0],/=,R);  #if RP_3D   N3V_CROSS(x[1],x[2],x[0]);  #else     x[1][0] = -x[0][1];     x[1][1] = x[0][0];  #endif   /* dmatrix is computed as xT*cond*x */   dmatrix[0][0] = cond[0]*x[0][0]*x[0][0]      + cond[1]*x[1][0]*x[1][0]  #if RP_3D    + cond[2]*x[2][0]*x[2][0]  #endif  ;   dmatrix[1][1] = cond[0]*x[0][1]*x[0][1]    + cond[1]*x[1][1]*x[1][1]  #if RP_3D    + cond[2]*x[2][1]*x[2][1]  #endif    ;   dmatrix[1][0] = cond[0]*x[0][1]*x[0][0]    + cond[1]*x[1][1]*x[1][0]  #if RP_3D    + cond[2]*x[2][1]*x[2][0]  #endif  ;   dmatrix[0][1] = dmatrix[1][0];
  #if RP_3D   dmatrix[2][2] = cond[0]*x[0][2]*x[0][2]    + cond[1]*x[1][2]*x[1][2]    + cond[2]*x[2][2]*x[2][2]    ;   dmatrix[0][2] = cond[0]*x[0][0]*x[0][2]    + cond[1]*x[1][0]*x[1][2]    + cond[2]*x[2][0]*x[2][2]    ;   dmatrix[2][0] = dmatrix[0][2];   dmatrix[1][2] = cond[0]*x[0][1]*x[0][2]    + cond[1]*x[1][1]*x[1][2]    + cond[2]*x[2][1]*x[2][2]     ;   dmatrix[2][1] = dmatrix[1][2];  #endif
}