constshot.c

Go to the documentation of this file.

/*
 *   
 * @file constshot.c
 * @author Juan Pablo Ibanez.
 * @date September 14, 2006
 *
 * Module created to support the use of SHOT constrain
 * (stolen from constcurve.c)
 *
 */

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

typedef struct _curvedata 
{
  double jac[3][3];
  double ijac[3][3];
} sCCurveData;

sCoord* shotVector = NULL;

/*
** ------------------------------------------------------------------------
** 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     DiffVec( sCoord *in0, sCoord *in1, sCoord *out )    
{
  out->x = in0->x - in1->x;
  out->y = in0->y - in1->y;
  out->z = in0->z - in1->z;
} 

static double ProdEscalar( sCoord *in0, sCoord* in1 )
{
  double prod;
  prod = (in0->x)*(in1->x) + (in0->y)*(in1->y) + (in0->z)*(in1->z); 
  
  return (prod);
}

static void CrossProd( sCoord *in0, sCoord *in1, sCoord *out )
{
  out->x = (in0->y * in1->z) - (in0->z * in1->y); 
  out->y = (in0->z * in1->x) - (in0->x * in1->z); 
  out->z = (in0->x * in1->y) - (in0->y * in1->x); 
}  

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;

  DiffVec(in1, in0, &tmp);
  
  return(VecLen(&tmp));
}

static void EvaluateTolerance( sConstrain *c )
{
  int inc[3];
  double dist = 0.0;
  double auxTol = 0.0;
  int i;

  for(i = 0; i < c->numElems; i++)
  {
    ElmConnect(c->elems[i],inc);
    dist = Dist(&c->nodes[inc[1]-1].coord, &c->nodes[inc[0]-1].coord);
    auxTol = dist * 0.1;
    if(auxTol < Tolerance)
      Tolerance = auxTol;
  }

} /* End of EvaluateTolerance */

static void TranslateShot(sConstrain *c, sCoord* p)
{
  int i;
  int sign;
  sCoord* vec = shotVector;

  /* aplly tranlation to shot constrain */
  /* p = current node */
  sign = -1;
  for(i = 0; i < c->numNodes; i++)
  {
    c->nodes[i].coord.x = p->x + sign*vec->x;
    c->nodes[i].coord.y = p->y + sign*vec->y;
    c->nodes[i].coord.z = p->z + sign*vec->z;
    sign += 2; 
  }

} /* end of TranslateShot */

/*
 * Subclass methods:
 */
static void ConstShotNew   (int, int, int, sConstrain **, sConstrain **);
static void ConstShotFree  (sConstrain *);
static void ConstShotRead  (sConstrain *);
static void ConstShotBuild (sConstrain *);

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

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

 (*c) = (sConstrain *)malloc(sizeof(sConstrain));
 (*c)->type = SHOT;
 (*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[ContNumConstrainRead] = (*c);
 ContNumConstrainRead++;

} /* End of ConstShotNew */

static void ConstShotFree(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 ConstShotFree */

static void ConstShotRead(sConstrain *c)
{
 int i;
 int sign;
 sElement *elm = NULL;
 sCoord   *vec = NULL;    /* vetor unitario */
 int       inc[2];        /* Vetor auxiliar de incidencia dos segmentos. */

  /* Le as coordenadas do vetor unitario.
   */
   vec = (sCoord*)malloc(sizeof(sCoord));
   fscanf(nf, "%lf %lf %lf", &(vec->x), &(vec->y), &(vec->z));
   shotVector = vec;
    
  /* inicia informacao comum aos dois nos */
   for(i = 0; i < c->numNodes; i++)
   {  
     c->nodes[i].id = i+1;
     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;
   }

  /* seta coordenadas dos dois nos */
   sign = -1;
   for(i = 0; i < c->numNodes; i++)
   {
     c->nodes[i].coord.x = sign*vec->x;
     c->nodes[i].coord.y = sign*vec->y;
     c->nodes[i].coord.z = sign*vec->z;
     sign += 2;
   }

  /* Cria os segmentos.
   */
  for(i = 0; i < c->numElems; i++)
  {
    ElmNew(LINE2, i+1, 1, 1, 1, &elm, c->elems, c->nodes);
    inc[0] = i+1;
    inc[1] = i+2;
    ElmSetConnect(elm, inc);
  }

 EvaluateTolerance(c);

 /*for(i = 0; i < c->numNodes; i++) {
   printf("node = %d\tx = %f\ty = %f\n",
          i,c->nodes[i].coord.x,c->nodes[i].coord.y);
 }*/

} /* End of ConstShotRead */

static void ConstShotBuild(sConstrain *c)
{
 int i, j, n, k, l;
 int inc[3];
 double jac[3][3], ijac[3][3];
 int *count = NULL;
 sCCurveData *data = NULL;

 count = (int *)calloc(c->numNodes, sizeof(int));
 data = (sCCurveData *)c->data;

 for( i = 0; i < c->numElems; i++ ) {
  ElmNumNodes(c->elems[i],&n);
  ElmConnect(c->elems[i],inc);
  ElmJacobian(c->elems[i],jac,ijac);
  for( j = 0; j < n; j++ ) {
   for( k = 0; k < 3; k++ ) {
    for( l = 0; l < 3; l++ ) {
     data[inc[j]-1].jac[k][l] += jac[k][l];
     data[inc[j]-1].ijac[k][l] += ijac[k][l];
    }
   }
   count[inc[j]-1]++;
  }
 }

 for( i = 0; i < c->numNodes; i++ ) {
  for( k = 0; k < 3; k++ ) {
   for( l = 0; l < 3; l++ ) {
    data[i].jac[k][l] /= count[i];
    data[i].ijac[k][l] /= count[i];
   }
  }
 }
 
 free(count);

} /* End of ConstShotBuild */

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

  /* update shot constrain */
  TranslateShot(c,p);

  ElmJacobian(c->elems[0],jac,ijac);

  *elm = c->elems[0];

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

  return 1;

} /* End of ConstShotLocal */

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

  ElmJacobian(c->elems[0],jac,ijac);

  *elm = c->elems[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 ConstShotGlobal */

static int ConstShotSlide( sConstrain *c, sCoord *p, void **elm, sCoord *d,
                           double *l, double *g, sCoord *pl )
{
  int      i, j;
  double   jace[3][3], ijace[3][3];
  double   d_gbl[3];
  double   d_local[3];

  d_local[0] = d->x;
  d_local[1] = 0.0;
  d_local[2] = 0.0;
      
  *elm = c->elems[0];
  ElmJacobian( ((sElement*)(*elm)), jace, ijace );

  /* Computes displacements in global coordinates */
  for(i = 0; i < NDof; i++)
  {
    d_gbl[i] = 0.0;
    for(j = 0; j < NDof; j++)
      d_gbl[i] += ijace[i][j] * d_local[j];
  }

  /* Computes new node coordinates in shot constrain. */
  pl->x = p->x + d_gbl[0]; 
  pl->y = p->y + d_gbl[1]; 
  if( NDof == 3 )
    pl->z = p->z + d_gbl[2]; 

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

  return 1;
}

/*
 * Public functions:
 */

void ConstShotInit(void);
void ConstShotInit(void)
{
  ConstClass[SHOT].new    = ConstShotNew;
  ConstClass[SHOT].free   = ConstShotFree;
  ConstClass[SHOT].read   = ConstShotRead;
  ConstClass[SHOT].build  = ConstShotBuild;
  ConstClass[SHOT].local  = ConstShotLocal;
  ConstClass[SHOT].global = ConstShotGlobal;
  ConstClass[SHOT].slide  = ConstShotSlide;

  Tolerance = 1.0e+10;

} /* End of ConstShotInit */


/*
 * End of File.
 */

---