constsurface.c

Go to the documentation of this file.

/*
 * Global variables and symbols:
 */
#include <stdlib.h>
#include <math.h>
#include "constrain.h"
#include "rio.h"

#include "geo.h"
#include "pfs.h"

typedef struct _surfacedata 
{
    void   *pfsSurface;
} sCSurfaceData;

/*
** ------------------------------------------------------------------------
** Local variables and symbols:
*/

/*
%K Tolerance definition
*/
static double Tolerance = 1.0e+10;

/*
%K Zero definition
*/
#ifndef ZERO
#define ZERO(v)     ((fabs(v)<Tolerance)?1:0) 
#endif

/*
%K Minimum definition
*/
#ifndef MIN
#define MIN(a,b)    ((a<b)?a:b)
#endif

/*
%K Maximum definition
*/
#ifndef MAX
#define MAX(a,b)    ((a>b)?a:b)
#endif


/*
** ------------------------------------------------------------------------
** Local functions:
*/

static void Inverse( double m[3][3], double d[3][3] )
{
 double det = m[0][0] * ((m[1][1] * m[2][2]) - (m[1][2] * m[2][1])) -
              m[0][1] * ((m[1][0] * m[2][2]) - (m[1][2] * m[2][0])) +
              m[0][2] * ((m[1][0] * m[2][1]) - (m[1][1] * m[2][0]));

 d[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]) / det;
 d[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]) / det;
 d[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]) / det;
 d[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]) / det;
 d[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]) / det;
 d[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]) / det;
 d[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]) / det;
 d[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]) / det;
 d[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]) / det;

} /* End of Inverse */

static double VecLen( sCoord *vec )
{
  return(sqrt((vec->x * vec->x) + (vec->y * vec->y) + (vec->z * vec->z)));
}

static double Dist( sCoord *in0, sCoord *in1 )
{
  sCoord tmp;

  tmp.x = in1->x - in0->x;
  tmp.y = in1->y - in0->y;
  tmp.y = in1->z - in0->z;

  return(VecLen(&tmp));
}

static void EvaluateTolerance( sConstrain *c )
{
  int inc[3];
  sCSurfaceData *data = NULL;
  double dist1 = 0.0;
  double dist2 = 0.0;
  double dist3 = 0.0;
  double auxTol1 = 0.0;
  double auxTol2 = 0.0;
  double auxTol3 = 0.0;
  int i;

  data = (sCSurfaceData *)c->data;

  for(i = 0; i < c->numElems; i++)
  {
    ElmConnect(c->elems[i],inc);
    dist1 = Dist(&c->nodes[inc[1]-1].coord, &c->nodes[inc[0]-1].coord);
    dist2 = Dist(&c->nodes[inc[2]-1].coord, &c->nodes[inc[1]-1].coord);
    dist3 = Dist(&c->nodes[inc[0]-1].coord, &c->nodes[inc[2]-1].coord);
    auxTol1 = dist1 * 0.1;
    auxTol2 = dist2 * 0.1;
    auxTol3 = dist3 * 0.1;
    if(auxTol1 < Tolerance)
      Tolerance = auxTol1;
    if(auxTol2 < Tolerance)
      Tolerance = auxTol2;
    if(auxTol3 < Tolerance)
      Tolerance = auxTol3;
  }

} /* End of EvaluateTolerance */

/* =====================  _surfaceCurrElemTransfMtx =================== */
static int _surfaceCurrElemTransfMtx( double transf[3][3] )
{
  int         i;
  int         nnode = 3; /* it works only for T3 elements */
  PFSGeoPoint pn[3], n[3];
 
  for(i = 0; i < nnode; i++)
  {
    int id;
    pfsRItrElemNode( &id );
    pfsRItrNodeCoords(id, &(pn[i].x), &(pn[i].y), &(pn[i].z));  
  }

  /* First direction, oriented vector from node 1 to node 2 of the element.
  */
  PFSGeoDiffVec(&(pn[1]), &(pn[0]), &(n[0]));

  /* Third direction, normal vector of the element surface.
   */
  pfsRItrElemNormal(&(n[2].x), &(n[2].y), &(n[2].z));

  /* Second direction.
   */
  PFSGeoCrossProd(&n[2], &n[0], &n[1]);

  /* Normalize the vectors.
   */
  PFSGeoVecNormalize(&n[0], &n[0]);
  PFSGeoVecNormalize(&n[1], &n[1]);
  PFSGeoVecNormalize(&n[2], &n[2]);

  /* Return the transformation matrix.
   */
  for(i=0; i<nnode; ++i)
  {
    transf[i][0] = n[i].x;
    transf[i][1] = n[i].y;
    transf[i][2] = n[i].z;
  }

  return 1;
}

static int _surfaceToLocal( sConstrain *c, sCoord *p, 
                            void **elm, double transf[3][3] )
{
  int found = 0;
  sCSurfaceData *data = (sCSurfaceData *)c->data;

  pfsRActivateSurf(data->pfsSurface);

  if( !(*elm) )
  {
    /* Traverses Pfs's elements for findind in which one the given point is */
    found = pfsRGetElem( p->x, p->y, p->z, elm );

    if(!found)
    {
      /* o pto não está sobre nenhum elemento, então procura o elemento
       * da borda mais próximo */
      found = pfsRGetLoopElem( p->x, p->y, p->z, elm );
      if( !found )
      {
        fprintf(stderr, "Could not find element(ToLocal).\n\n");
        *elm = NULL; /* not necessary */
        return 0;
      }
    }
  }

  /* sets current pfs element */
  pfsRSetCurrElem( *elm );

  /* gets local element transformation matrix */
  _surfaceCurrElemTransfMtx( transf ); /* not using the relax jacobian */

  return 1;
}
static void _invertMtx( double m[3][3], double i[3][3] )
{
  double det = m[0][0] * ((m[1][1] * m[2][2]) - (m[1][2] * m[2][1])) -
               m[0][1] * ((m[1][0] * m[2][2]) - (m[1][2] * m[2][0])) +
               m[0][2] * ((m[1][0] * m[2][1]) - (m[1][1] * m[2][0]));

  i[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]) / det;
  i[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]) / det;
  i[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]) / det;
  i[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]) / det;
  i[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]) / det;
  i[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]) / det;
  i[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]) / det;
  i[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]) / det;
  i[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]) / det;

}

static int _surfaceToGlobal( sConstrain *c, sCoord *p,
                             void **elm, double itransf[3][3] )
{
  int    found = 0;
  double transf[3][3];
  sCSurfaceData *data = (sCSurfaceData *)c->data;

  pfsRActivateSurf(data->pfsSurface);

  if( !(*elm) )
  {
    /* Traverses Pfs's elements for findind in which one the given point is */
    found = pfsRGetElem( p->x, p->y, p->z, elm );

    if(!found)
    {
      /* o pto não está sobre nenhum elemento, então procura o elemento
       * da borda mais próximo */
      found = pfsRGetLoopElem( p->x, p->y, p->z, elm );
      if( !found )
      {
        fprintf(stderr, "Could not find element(ToGlobal).\n\n");
        *elm = NULL; /* not necessary */
        return 0;
      }
    }
  }

  /* sets current pfs element */
  pfsRSetCurrElem( *elm );

  /* gets local element transformation matrix */
  _surfaceCurrElemTransfMtx( transf ); /* not using the relax jacobian */

  /* inverts local element transformation matrix */
  _invertMtx( transf, itransf );

  return 1;
} /* End of _surfaceToGlobal */

/* Evaluate the coordinates defined in the local system of the current 
 * element. 
 */
static void _surfaceLocalT3Coord( int inc[3], sCoord coord[3], double jac[3][3] )
{
  const int    nnode = 3;
        int    i;
  PFSGeoPoint  tmp;
  PFSGeoPoint  pn[3], n[3];

  for(i = 0; i < nnode; i++)
  {
    pfsRItrNodeCoords(inc[i], &(pn[i].x), &(pn[i].y), &(pn[i].z));
  }

  /* First direction, oriented vector from node 1 to node 2 of the element.
  */
  PFSGeoDiffVec(&(pn[1]), &(pn[0]), &(n[0]));

  /* Third direction, normal vector of the element surface.
   */
  pfsRItrElemNormal(&(n[2].x), &(n[2].y), &(n[2].z));

  /* Second direction.
   */
  PFSGeoCrossProd(&n[2], &n[0], &n[1]);

  /* Normalize the vectors.
   */
  for(i=0;i<nnode;++i)
  {
    PFSGeoVecNormalize(&n[i], &n[i]);
  }

  /* node 1 */
  coord[0].x = 0.0;
  coord[0].y = 0.0;
  coord[0].z = 0.0;

  /* node 2 */
  PFSGeoDiffVec(&(pn[1]), &(pn[0]), &tmp);
  coord[1].x = PFSGeoVecLen(&tmp);
  coord[1].y = 0.0;
  coord[1].z = 0.0;

  /* node 3 */
  PFSGeoDiffVec(&(pn[2]), &(pn[0]), &tmp);
  coord[2].x = PFSGeoDotProd(&n[0], &tmp);
  coord[2].y = PFSGeoDotProd(&n[1], &tmp);
  coord[2].z = 0.0;

  /* Preenche a matriz q transforma do sistema global para o sistema
   * local do elemento 
   */
  for( i=0; i<3; i++ )
  {
    jac[i][0] = n[i].x;
    jac[i][1] = n[i].y;
    jac[i][2] = n[i].z;
  }
}

static int _surfaceLocalTri2PFS( sCoord coord[3], 
                                 double u1, double v1,
                                 double u2, double v2, 
                                 double u3, double v3, double jacinv[2][2] )
{
  double x1 = coord[0].x; double y1 = coord[0].y;
  double x2 = coord[1].x; double y2 = coord[1].y;
  double x3 = coord[2].x; double y3 = coord[2].y;

  double dt;
  double dudx;
  double dudy;
  double dvdx;
  double dvdy;
 
  dt = (x1*y2 - y1*x2) + (x2*y3 - y2*x3) + (x3*y1 - y3*x1);

  if( dt < 1e-8 )
    return( 0 ); 
   
  dudx = (y2 - y3)*u1 + (y3 - y1)*u2 + (y1 - y2)*u3;
  dudy = (x3 - x2)*u1 + (x1 - x3)*u2 + (x2 - x1)*u3;
  dvdx = (y2 - y3)*v1 + (y3 - y1)*v2 + (y1 - y2)*v3;
  dvdy = (x3 - x2)*v1 + (x1 - x3)*v2 + (x2 - x1)*v3;
  
  jacinv[0][0] = dudx / dt;
  jacinv[0][1] = dudy / dt;
  jacinv[1][0] = dvdx / dt;
  jacinv[1][1] = dvdy / dt;
  
  return( 1 );
}
static int _surfaceSlide( sConstrain *c, sCoord *p, void **elm, sCoord *d,
                          double ijac[3][3], sCoord *pl )
{
  int    i;
  int    nnode = 3;
  int    nodeid[3];
  double du_pfs, dv_pfs;
  double u[3], v[3];
  double jac[3][3];       /* transforma do sistema global para o 
                           * local do elemento */ 
  double local2pfs[2][2]; /* transforma do sistema local para o
                           * sistema da pfs (u,v) */
  double path_len;
  sCoord localT3Coord[3];
  int    status;

  sCSurfaceData *data = (sCSurfaceData *)c->data;

  if( !(*elm) )
  {
    /* Traverses Pfs's elements to find in which one the given point is */
    int found = pfsRGetElem( p->x, p->y, p->z, elm );

    if( !found )
    {
      /* o pto não está sobre nenhum elemento, então procura o elemento
       * da borda mais próximo */
      found = pfsRGetLoopElem( p->x, p->y, p->z, elm );
      if( !found )
      {
        fprintf(stderr, "Could not find element(Slide).\n\n");
        *elm = NULL; /* not necessary */
        return 0;
      }
    }
  }

  /* Sets the given surface mesh as the current active one.
   */
  pfsRActivateSurf(data->pfsSurface);

  /* Sets the PFS current element. */
  pfsRSetCurrElem(*elm);

  /* Each iteration:
   * gets PFS node indices of current element;
   * gets PFS parametric values of nodes of current element.
   */
  for( i=0; i< nnode; ++i )
  {
    pfsRItrElemNode(&nodeid[i]);
    pfsRItrNodeParVals(nodeid[i], &u[i], &v[i]);
  }

  /* Computes local T3 element coordinate nodes */
  _surfaceLocalT3Coord( nodeid, localT3Coord, jac );
  
  Inverse( jac, ijac );
 
  /* Computes local2pfs */
  _surfaceLocalTri2PFS( localT3Coord, u[0], v[0], 
                                      u[1], v[1],
                                      u[2], v[2], local2pfs );

  du_pfs = local2pfs[0][0]*d->x + local2pfs[0][1]*d->y;
  dv_pfs = local2pfs[1][0]*d->x + local2pfs[1][1]*d->y;

  path_len = sqrt( d->x*d->x + d->y*d->y );

  status = pfsRShot( *elm, p->x, p->y, p->z,
                     du_pfs, dv_pfs, path_len,
                     elm, &pl->x, &pl->y, &pl->z );

  if( !status )
  {
   pl->x = p->x;
   pl->y = p->y;
   pl->z = p->z;
  }
  return 1;
} /* End of SurfaceToSlide */

/*
 * Subclass methods:
 */
static void ConstSurfaceNew(int, int, int, sConstrain **, sConstrain **);
static void ConstSurfaceFree(sConstrain *);
static void ConstSurfaceRead(sConstrain *);
static void ConstSurfaceBuild(sConstrain *);

static void ConstSurfaceNew( int label, int num_nodes, int num_elems,
                             sConstrain **c, sConstrain **clst)
{
  sCSurfaceData *data = NULL;

  data = (sCSurfaceData *)calloc(num_nodes, sizeof(sCSurfaceData));

  (*c) = (sConstrain *)calloc(1, sizeof(sConstrain));
  (*c)->type = SURFACE;
  (*c)->label = label;
  (*c)->numNodes = num_nodes;
  (*c)->numElems = num_elems;
  (*c)->nodes = (sNode *)calloc(num_nodes, sizeof(sNode));
  (*c)->elems = (sElement **)calloc(num_elems, sizeof(sElement *));
  (*c)->data = data;

  //clst[label-1] = (*c);
  //clst[i] = (*c);

  clst[ContNumConstrainRead] = (*c);
  ContNumConstrainRead++;

} /* End of ConstSurfaceNew */

static void ConstSurfaceFree( sConstrain *c )
{
  int i;

  if(c->data) 
    free(c->data);

  if(c->numNodes != 0)
    free(c->nodes);

  if(c->numElems != 0) 
  {
    for(i = 0; i < c->numElems; i++) 
    {
      ElmFree(c->elems[i]);
    }
    free(c->elems);
  }
  c->numNodes = 0;
  c->numElems = 0;
  c->nodes = NULL;
  c->elems = NULL;
  c->data = NULL;

} /* End of ConstSurfaceFree */

static void ConstSurfaceRead( sConstrain *c )
{
  int i, id;
  double x, y, z;
  int matid, tckid, ordid;
  sElement *elm = NULL;

  /* Lê os nós da superfície de restrição. 
   */
  for(i = 0; i < c->numNodes; i++) 
  {
    fscanf(nf, "%d %lf %lf %lf", &id, &x, &y, &z);
    c->nodes[i].id = id;
    c->nodes[i].coord.x = x;
    c->nodes[i].coord.y = y;
    c->nodes[i].coord.z = z;
    c->nodes[i].dof.x = FORCE;
    c->nodes[i].dof.y = FORCE;
    c->nodes[i].dof.z = FORCE;
    c->nodes[i].dof.psi = 0.0;
    c->nodes[i].rezone = NONE;
    c->nodes[i].curve = 0;
    c->nodes[i].dspIteration = DSP_NO;
    c->nodes[i].dspTime = DSP_NO;
  }

  /* Lê os elementos da superfície de restrição.
   */
  for(i = 0; i < c->numElems; i++) 
  {
    fscanf(nf, "%d%d%d%d", &id, &matid, &tckid, &ordid);
    ElmNew(T3, id, matid, ordid, tckid, &elm, c->elems, c->nodes);
    ElmRead(elm);
  }

  EvaluateTolerance(c);

} /* End of ConstSurfaceRead */

static void ConstSurfaceBuild( sConstrain *c )
{
  sCSurfaceData *data = NULL;
  int            inc[3];
  int            i;

  data = (sCSurfaceData *)c->data;
 
  /* Creates the pfs surface data structure.
   */
  data->pfsSurface = pfsRInitSurf(c->numNodes, c->numElems);

  /* Sets the given surface mesh as the current active one.
   */
  pfsRActivateSurf(data->pfsSurface);

  /* Inserts the nodes in the database of the current active surface.
   */
  for(i = 0; i < c->numNodes; i++) 
    pfsRAddNode(c->nodes[i].coord.x, c->nodes[i].coord.y, c->nodes[i].coord.z);

  /* Inserts the elements (faces) int the database of the current
   * active surface.
   */
  for(i = 0; i < c->numElems; i++) 
  {
    /* Gets the element connectivites.
     */
    ElmConnect(c->elems[i],inc);

    /* Inserts the current element. 
     */
    pfsRAddElem(inc[0]-1, inc[1]-1, inc[2]-1);
  }

  /* Finishes the creation of the data structure of the current active
   * surface mesh.
   */
  pfsRCompleteSurf();

  /* Computes by a Gauss-Seidel iterative scheme the surface parametric
   * values of all nodes of the mesh in current active surface.
   */
  pfsRSolveSurfPar(0);
  
  /* finishes up definition of surface parametric values of all nodes 
   * of the mesh, computes bounding box for each element and its respactivy
   * Jacobian matrix */
  pfsCompleteParSpace();

} /* End of ConstSurfaceBuild */

static int ConstSurfaceLocal( sConstrain *c, sCoord *p, void **elm,
                              double *g, double *l )
{
  int    i, j;
  double transf[3][3];

  if(!_surfaceToLocal(c, p, elm, transf))
  {
    return 0;
  }

  for(i = 0; i < NDof; i++) 
  {
    l[i] = 0.0;
    for(j = 0; j < NDof; j++) 
    {
      l[i] += transf[i][j] * g[j];
    }
  }

  return 1;

} /* End of ConstSurfaceLocal */

static int ConstSurfaceGlobal( sConstrain *c, sCoord *p, void **elm, 
                               double *l, double *g )
{
  int i, j;
  double itransf[3][3]; /* not using relax jacobian matrix */

  if(!_surfaceToGlobal(c, p, elm, itransf))
    return 0;

  for( i = 0; i < NDof; i++ ) 
  {
    g[i] = 0.0;
    for( j = 0; j < NDof; j++ ) 
    {
      g[i] += itransf[i][j] * l[j];
    }
  }

  return 1;

} /* End of ConstSurfaceGlobal */

static int ConstSurfaceSlide( sConstrain *c, sCoord *p, void **elm, sCoord *d,
                               double *l, double *g, sCoord *pl )
{
  int i, j;
  double ijac[3][3];

  if(!_surfaceSlide(c, p, elm, d, ijac, pl))
    return 0;

  for(i = 0; i < NDof; i++)
  {
    g[i] = 0.0;
    for(j = 0; j < NDof; j++)
    {
      g[i] += ijac[i][j] * l[j];
    }
  }

  return 1;

} /* End of ConstSurfaceSlide */

/*
 * Public functions:
 */

void ConstSurfaceInit(void);
void ConstSurfaceInit(void)
{
 ConstClass[SURFACE].new    = ConstSurfaceNew;
 ConstClass[SURFACE].free   = ConstSurfaceFree;
 ConstClass[SURFACE].read   = ConstSurfaceRead;
 ConstClass[SURFACE].build  = ConstSurfaceBuild;
 ConstClass[SURFACE].local  = ConstSurfaceLocal;
 ConstClass[SURFACE].global = ConstSurfaceGlobal;
 ConstClass[SURFACE].slide  = ConstSurfaceSlide;

 Tolerance = 1.0e+10;

} /* End of ConstSurfaceInit */

/*
 * End of File.
 */

---