load.c

Go to the documentation of this file.

/*
%M This modules contains the load methods and definitions
%a Joao Luiz Elias Campos.
%d September 2nd, 1998.
%r $Id: load.c,v 1.1 2004/06/22 05:29:59 joaoluiz Exp $
%w (C) COPYRIGHT 1995-1996, Eduardo Nobre Lages.
   (C) COPYRIGHT 1997-1999, Joao Luiz Elias Campos.
   All Rights Reserved
   Duplication of this program or any part thereof without the express
   written consent of the author is prohibited.
*/

/*
** ------------------------------------------------------------------------
** Global variables and symbols:
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "load.h"
#include "elm.h"
#include "node.h"
#include "alg.h"

/*
%V Load case and its numbers
*/
int        CurrLoadCase = 0;
int        NumLoadCase  = 0;
sLoadCase *LoadCase     = 0L;


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

#ifndef ZERO
#define ZERO 0.000001
#endif

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


/*
** ------------------------------------------------------------------------
** Public functions:
*/

/*
%F This function apply the prescribed values.
*/
void PrescribedValues( void )
{
 int        i, j, node, key;
 int        noi, noj, nok, nol, id, nnode;
 int        inc[8];
 double     vx, vy, vz;
 double     q1x, q1y, q1z;
 double     q2x, q2y, q2z;
 double     q3x, q3y, q3z;
 double     q4x, q4y, q4z;
 double     intforce[24];
 eLoadType  type;
 sElement  *elm = 0L;
 

/* Initialize prescribed values
 */
 NodeInitPrescValues( );

/* Update nodal support
 */
 for( i = 0; i < NumNodalSupp; i++ )
 {
  id = NodalSupp[i].node - 1;
  if( NodalSupp[i].tx == 0 )
   NodeVector[id].dof.x = FORCE;
  if( NodalSupp[i].ty == 0 ) 
   NodeVector[id].dof.y = FORCE;
  if( NodalSupp[i].tz == 0 ) 
   NodeVector[id].dof.z = FORCE;
 }

/* Apply nodal forces
 */
 for( i = 0; i < NumNodalForce; i++ )
 {
  node = NodalForce[i].node;
  vx   = NodalForce[i].fx;
  vy   = NodalForce[i].fy;
  vz   = NodalForce[i].fz;
  NodeApplyForce( node, vx, vy, vz );
 }

/* Apply nodal displacement
 */
 for( i = 0; i < NumNodalDisp; i++ )
 {
  node = NodalDisp[i].node;
  vx   = NodalDisp[i].dx;
  vy   = NodalDisp[i].dy;
  vz   = NodalDisp[i].dz;
  NodeApplyDisplacement( node, vx, vy, vz );
 }

/* Apply nodal velocity
 */
 for( i = 0; i < NumNodalVel; i++ )
 {
  node = NodalVel[i].node;
  vx   = NodalVel[i].vx;
  vy   = NodalVel[i].vy;
  vz   = NodalVel[i].vz;
  NodeApplyVelocity( node, vx, vy, vz );
 }

/* Apply line distributed loads
 */
 if( NumLineForce )
  ElementBuildAdjacence( );

 for( i = 0; i < NumLineForce; i++ )
 {
#if 0
 /* Check for the element id
  */
  if( LineForce[i].elem )
  {
   if( !LineForce[i].adj )
    continue;
  }

 /* Check for the adjacent element id
  */
  if( LineForce[i].adj )
  {
   if( !LineForce[i].elem )
    continue;
  }
#endif

 /* Get load information
  */
  q1x  = LineForce[i].q1x;
  q1y  = LineForce[i].q1y;
  q1z  = LineForce[i].q1z;
  q2x  = LineForce[i].q2x;
  q2y  = LineForce[i].q2y;
  q2z  = LineForce[i].q2z;
  noi  = LineForce[i].noi;
  noj  = LineForce[i].noj;
  nok  = LineForce[i].nok;
  key  = LineForce[i].key;
  type = LineForce[i].type;

 /* Get element address
  */
  elm = ElmList[LineForce[i].elem-1];

 /* Compute the element contribution
  */
  ElmLoad(elm,type,key,noi,noj,nok,0L,&q1x,&q1y,&q1z,
          &q2x,&q2y,&q2z,0L,0L,0L,0L,0L,0L);
 /* Add to the prescribed values vector
  */
  if( NodeVector[noi-1].dof.x == FORCE )
   NodeVector[noi-1].dof.vpx += q1x;
  if( NodeVector[noi-1].dof.y == FORCE )
   NodeVector[noi-1].dof.vpy += q1y;
  if( NodeVector[noi-1].dof.z == FORCE )
   NodeVector[noi-1].dof.vpz += q1z;

  if( NodeVector[noj-1].dof.x == FORCE )
   NodeVector[noj-1].dof.vpx += q2x;
  if( NodeVector[noj-1].dof.y == FORCE )
   NodeVector[noj-1].dof.vpy += q2y;
  if( NodeVector[noj-1].dof.z == FORCE )
   NodeVector[noj-1].dof.vpz += q2z;

  if( (NDof == 3) && (NodeVector[nok-1].dof.x == FORCE) )
   NodeVector[nok-1].dof.vpx += q1x;
  if( (NDof == 3) && (NodeVector[nok-1].dof.y == FORCE) )
   NodeVector[nok-1].dof.vpy += q1y;
  if( (NDof == 3) && (NodeVector[nok-1].dof.z == FORCE) )
   NodeVector[nok-1].dof.vpz += q1z;
 }
 
 /* Apply area uniform loads
 */
 if( NumAreaForceUniform )
  ElementBuildAdjacence( );

 for( i = 0; i < NumAreaForceUniform; i++ )
 {
 /* Get load information
  */
  q1x   = AreaForceUniform[i].qx;
  q1y   = AreaForceUniform[i].qy;
  q1z   = AreaForceUniform[i].qz;
  noi  = AreaForceUniform[i].noi;
  noj  = AreaForceUniform[i].noj;
  nok  = AreaForceUniform[i].nok;
  key  = AreaForceUniform[i].key;
  type = AreaForceUniform[i].type;

 /* Get element address
  */
  elm = ElmList[AreaForceUniform[i].elem-1];

 /* Compute the element contribution
  */
  ElmLoad(elm,type,key,noi,noj,nok,&nol,&q1x,&q1y,&q1z,
          &q2x,&q2y,&q2z,&q3x,&q3y,&q3z,&q4x,&q4y,&q4z);

 /* Add to the prescribed values vector
  */
  if( NodeVector[noi-1].dof.x == FORCE )
   NodeVector[noi-1].dof.vpx += q1x;
  if( NodeVector[noi-1].dof.y == FORCE )
   NodeVector[noi-1].dof.vpy += q1y;
  if( NodeVector[noi-1].dof.z == FORCE )
   NodeVector[noi-1].dof.vpz += q1z;

  if( NodeVector[noj-1].dof.x == FORCE )
   NodeVector[noj-1].dof.vpx += q2x;
  if( NodeVector[noj-1].dof.y == FORCE )
   NodeVector[noj-1].dof.vpy += q2y;
  if( NodeVector[noj-1].dof.z == FORCE )
   NodeVector[noj-1].dof.vpz += q2z;

  if( NodeVector[nok-1].dof.x == FORCE )
   NodeVector[nok-1].dof.vpx += q3x;
  if( NodeVector[nok-1].dof.y == FORCE )
   NodeVector[nok-1].dof.vpy += q3y;
  if( NodeVector[nok-1].dof.z == FORCE )
   NodeVector[nok-1].dof.vpz += q3z;
   
  if( nol != (-1)){
   AreaForceUniform[NumAreaForceUniform].nol  = nol;
   if( NodeVector[nol-1].dof.x != DISPLACEMENT ){
    NodeVector[nol-1].dof.x = FORCE;
    NodeVector[nol-1].dof.vpx += q4x;
   }
   if( NodeVector[nol-1].dof.y != DISPLACEMENT ){
    NodeVector[nol-1].dof.y = FORCE;
    NodeVector[nol-1].dof.vpy += q4y;
   }
   if( NodeVector[nol-1].dof.z != DISPLACEMENT ){
    NodeVector[nol-1].dof.z = FORCE;
    NodeVector[nol-1].dof.vpz += q4z;
   }
  }
  
 }
 
 /* Apply area pressure loads
 */
 if( NumAreaForcePressure )
  ElementBuildAdjacence( );

 for( i = 0; i < NumAreaForcePressure; i++ )
 {
 /* Get load information
  */
  q1x   = AreaForcePressure[i].qx;
  q1y   = AreaForcePressure[i].qy;
  q1z   = AreaForcePressure[i].qz;
  noi  = AreaForcePressure[i].noi;
  noj  = AreaForcePressure[i].noj;
  nok  = AreaForcePressure[i].nok;
  key  = AreaForcePressure[i].key;
  type = AreaForcePressure[i].type;

 /* Get element address
  */
  elm = ElmList[AreaForcePressure[i].elem-1];

 /* Compute the element contribution
  */
  ElmLoad(elm,type,key,noi,noj,nok,&nol,&q1x,&q1y,&q1z,
          &q2x,&q2y,&q2z,&q3x,&q3y,&q3z,&q4x,&q4y,&q4z);

 /* Add to the prescribed values vector
  */
  if( NodeVector[noi-1].dof.x == FORCE )
   NodeVector[noi-1].dof.vpx += q1x;
  if( NodeVector[noi-1].dof.y == FORCE )
   NodeVector[noi-1].dof.vpy += q1y;
  if( NodeVector[noi-1].dof.z == FORCE )
   NodeVector[noi-1].dof.vpz += q1z;

  if( NodeVector[noj-1].dof.x == FORCE )
   NodeVector[noj-1].dof.vpx += q2x;
  if( NodeVector[noj-1].dof.y == FORCE )
   NodeVector[noj-1].dof.vpy += q2y;
  if( NodeVector[noj-1].dof.z == FORCE )
   NodeVector[noj-1].dof.vpz += q2z;

  if( NodeVector[nok-1].dof.x == FORCE )
   NodeVector[nok-1].dof.vpx += q3x;
  if( NodeVector[nok-1].dof.y == FORCE )
   NodeVector[nok-1].dof.vpy += q3y;
  if( NodeVector[nok-1].dof.z == FORCE )
   NodeVector[nok-1].dof.vpz += q3z;
   
  if( nol != (-1)){
   AreaForcePressure[NumAreaForcePressure].nol  = nol;
   if( NodeVector[nol-1].dof.x != DISPLACEMENT ){
    NodeVector[nol-1].dof.x = FORCE;
    NodeVector[nol-1].dof.vpx += q4x;
   }
   if( NodeVector[nol-1].dof.y != DISPLACEMENT ){
    NodeVector[nol-1].dof.y = FORCE;
    NodeVector[nol-1].dof.vpy += q4y;
   }
   if( NodeVector[nol-1].dof.z != DISPLACEMENT ){
    NodeVector[nol-1].dof.z = FORCE;
    NodeVector[nol-1].dof.vpz += q4z;
   }
  }
  
 }

/* Gravitational load
 */
 if( (fabs(GravForce.gx) > ZERO) || 
     (fabs(GravForce.gy) > ZERO) || 
     (fabs(GravForce.gz) > ZERO) )
 {
  for( i = 0; i < NumElements; i++ )
  {
   if( ElmList[i]->rezone != NONE )
    continue;

   vx = vy = vz = 0.0;

   ElmGravity( ElmList[i], &vx, &vy, &vz );
   ElmConnect( ElmList[i], inc );
   ElmNumNodes( ElmList[i], &nnode );

   for( j = 0; j < nnode; j++ )
   {
    if( NodeVector[inc[j]-1].dof.x == FORCE ) 
     NodeVector[inc[j]-1].dof.vpx += vx;
    if( NodeVector[inc[j]-1].dof.y == FORCE ) 
     NodeVector[inc[j]-1].dof.vpy += vy;
    if( (NDof == 3) && (NodeVector[inc[j]-1].dof.z == FORCE) ) 
     NodeVector[inc[j]-1].dof.vpz += vz;
   }
  }
 }

/* Initial stress
 */
 for( i = 0; i < NumElements; i++ )
 {
  if( ElmList[i]->rezone != NONE )
   continue;

  for( j = 0; j < 24; j++ )
   intforce[j] = 0.0;

  ElmInterForce( ElmList[i], 0L, intforce );
  ElmConnect( ElmList[i], inc );
  ElmNumNodes( ElmList[i], &nnode );

  for( j = 0; j < nnode; j++ )
  {
   if(NodeVector[inc[j]-1].dof.x != DISPLACEMENT){
    NodeVector[inc[j]-1].dof.x = FORCE;
    NodeVector[inc[j]-1].dof.vpx += intforce[NDof*j];
   }
   if(NodeVector[inc[j]-1].dof.y != DISPLACEMENT){
    NodeVector[inc[j]-1].dof.y = FORCE;
    NodeVector[inc[j]-1].dof.vpy += intforce[NDof*j+1];
   }
   if( (NDof == 3) && (NodeVector[inc[j]-1].dof.z != DISPLACEMENT) ){
    NodeVector[inc[j]-1].dof.z = FORCE;
    NodeVector[inc[j]-1].dof.vpz += intforce[NDof*j+2];
  }
 }
 }

/* Element internal pressure
 */
 for( i = 0; i < NumElemPoten; i++ )
 {
  if( ElmList[i]->rezone != NONE )
   continue;

  elm = ElmList[ElemPoten[i].elem-1];

  ElmSetPressure( elm, ElemPoten[i].phi );
 }

} /* End of PrescribedValues */


/*
%F
*/
void UpdatePrescribedValues( double *FVector )
{
 int        i, j, key;
 int        noi, noj, nok, nnode;
 int        inc[8];
 double     vx, vy, vz;
 double     q1x, q1y, q1z;
 double     q2x, q2y, q2z;
 eLoadType  type;
 sElement  *elm = 0L;

/* Update rezoned element contribution to the load vector.
 */
 for( i = 0; i < NumLineForce; i++ )
 {
 /* Check for the element id
  */
  if( (ElmList[LineForce[i].elem-1]->rezone == NONE) ||
      ((LineForce[i].adj > 0) && 
       (ElmList[LineForce[i].adj-1]->rezone == NONE)) ) continue;

 /* Get load information
  */
  q1x  = LineForce[i].q1x;
  q1y  = LineForce[i].q1y;
  q1z  = LineForce[i].q1z;
  q2x  = LineForce[i].q2x;
  q2y  = LineForce[i].q2y;
  q2z  = LineForce[i].q2z;
  noi  = LineForce[i].noi;
  noj  = LineForce[i].noj;
  nok  = LineForce[i].nok;
  key  = LineForce[i].key;
  type = LineForce[i].type;

 /* Get element address
  */
  elm = ElmList[LineForce[i].elem-1];

 /* Compute the element contribution
  */
  ElmLoad(elm,type,key,noi,noj,nok,0L,
          &q1x,&q1y,&q1z,&q2x,&q2y,&q2z,
          0L,0L,0L,0L,0L,0L);

 /* Add to the prescribed values vector
  */
  if( NodeVector[noi-1].dof.x == FORCE )
  {
   NodeVector[noi-1].dof.vpx -= q1x;
   FVector[NDof*(noi-1)]     -= q1x;
  }
  if( NodeVector[noi-1].dof.y == FORCE )
  {
   NodeVector[noi-1].dof.vpy -= q1y;
   FVector[NDof*(noi-1)+1]   -= q1y;
  }
  if( NodeVector[noi-1].dof.z == FORCE )
  {
   NodeVector[noi-1].dof.vpz -= q1z;
   FVector[NDof*(noi-1)+2]   -= q1z;
  }

  if( NodeVector[noj-1].dof.x == FORCE )
  {
   NodeVector[noj-1].dof.vpx -= q2x;
   FVector[NDof*(noj-1)]     -= q2x;
  }
  if( NodeVector[noj-1].dof.y == FORCE )
  {
   NodeVector[noj-1].dof.vpy -= q2y;
   FVector[NDof*(noj-1)+1]   -= q2y;
  }
  if( NodeVector[noj-1].dof.z == FORCE )
  {
   NodeVector[noj-1].dof.vpz -= q2z;
   FVector[NDof*(noj-1)+2]   -= q2z;
  }

  if( (NDof == 3) && (NodeVector[nok-1].dof.x == FORCE) )
  {
   NodeVector[nok-1].dof.vpx -= q1x;
   FVector[NDof*(nok-1)]     -= q1x;
  }
  if( (NDof == 3) && (NodeVector[nok-1].dof.y == FORCE) )
  {
   NodeVector[nok-1].dof.vpy -= q1y;
   FVector[NDof*(nok-1)+1]   -= q1y;
  }
  if( (NDof == 3) && (NodeVector[nok-1].dof.z == FORCE) )
  {
   NodeVector[nok-1].dof.vpz -= q1z;
   FVector[NDof*(nok-1)+2]   -= q1z;
  }
 }

/* Gravitational load
 */
 if( (fabs(GravForce.gz) > ZERO) || 
     (fabs(GravForce.gy) > ZERO) || 
     (fabs(GravForce.gz) > ZERO) )
 {
  for( i = 0; i < NumElements; i++ )
  {
   if( ElmList[i]->rezone != TODO )
    continue;

   vx = vy = vz = 0.0;

   ElmGravity( ElmList[i], &vx, &vy, &vz );
   ElmConnect( ElmList[i], inc );
   ElmNumNodes( ElmList[i], &nnode );

   for( j = 0; j < nnode; j++ )
   {
    if( NodeVector[inc[j]-1].dof.x == FORCE ) 
    {
     NodeVector[inc[j]-1].dof.vpx -= vx;
     FVector[NDof*(inc[j]-1)]     -= vx;
    }
    if( NodeVector[inc[j]-1].dof.y == FORCE ) 
    {
     NodeVector[inc[j]-1].dof.vpy -= vy;
     FVector[NDof*(inc[j]-1)+1]   -= vy;
    }
    if( (NDof == 3) && (NodeVector[inc[j]-1].dof.z == FORCE) ) 
    {
     NodeVector[inc[j]-1].dof.vpz -= vz;
     FVector[NDof*(inc[j]-1)+2]   -= vz;
    }
   }
  }
 }

} /* End of UpdatePrescribedValues */


/* =======================================================  End of File  */

---