iso_axi.c

C 开发的有限元软件

/*
    This file is part of the FElt finite element analysis package.
    Copyright (C) 1993-1995 Jason I. Gobat and Darren C. Atkinson

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/*****************************************************************************
*
* File:		iso_axi.c
*
* Description:	contains the element definition routines for isoparametric
*		plane stress / plane strain elements with only four nodes
*		(faster, simpler, etc.)
*
******************************************************************************/

# include <stdio.h>
# include <math.h>
# include "allocate.h"
# include "fe.h"
# include "error.h"
# include "misc.h"

# define FORCE	 1
# define NOFORCE 0

void     AQuadLumpedMassMatrix ( );
void     AQuadConsistentMassMatrix ( );
unsigned LocalAQuadShapeFunctions   ( );
Vector   GlobalAQuadShapeFunctions  ( );
Matrix   IsoAQuadLocalB ( );
Vector	 IsoAQuadEquivNodalForces ( );
int      AQuadElementSetup ( );
int  	 AQuadElementStress    ( );

int quad_AxisymEltSetup ( ), quad_AxisymEltStress ( );

struct definition quad_AxisymDefinition = {
   "quad_Axisym", quad_AxisymEltSetup, quad_AxisymEltStress,
   Planar, 4, 4, 6, 2, {0, 1, 2, 0, 0, 0, 0}, 0
};


int quad_AxisymEltSetup (element, mass_mode)
   Element	element;
   char		mass_mode;
{
   return AQuadElementSetup (element, mass_mode);
}

int quad_AxisymEltStress (element)
   Element	element;
{
   return AQuadElementStress (element);
}


int AQuadElementSetup (element, mass_mode)
   Element	element;
   char		mass_mode;
{
   double xc1,xc3,rr,maxx,minx,sg;
   unsigned		numnodes;
   unsigned		i,j;
   int			ninteg;
   Matrix		B;
   Matrix		D;
   Vector		jac;
   Vector		equiv;
   int			count;
   static Vector	weights,xpoints;
   static Matrix	tempK;
   static Matrix	N, dNdxi, dNde,
                        dNdx, dNdy = NullMatrix;
   static Matrix	Bt, temp;

   if (dNdy == NullMatrix) {
  
      N     = CreateMatrix (4,4);
      dNdxi = CreateMatrix (4,4);
      dNde  = CreateMatrix (4,4);
      dNdx  = CreateMatrix (4,4);
      dNdy  = CreateMatrix (4,4);
      weights = CreateVector (4);
      xpoints = CreateVector (8);
      tempK = CreateMatrix (8,8);
      Bt    = CreateMatrix (8,4);
      temp  = CreateMatrix (8,4);
   }

   if (element -> material -> E == 0) {
      error ("isoparametric element %d has 0.0 for Young's modulus (E)",element -> number);
      return 1;
   }
   if (element -> material -> nu == 0) {
      error ("isoparametric element %d has 0.0 for Poisson's ratio (nu)",element -> number);
      return 1;
   }

   ninteg = 4;	/* 2 x 2 quadrature */
  
   numnodes = LocalAQuadShapeFunctions (element, ninteg, N, 
                                       dNdxi, dNde, weights,xpoints, NOFORCE);  

   jac = GlobalAQuadShapeFunctions (element,dNdxi,dNde,dNdx,dNdy,
                                   ninteg,numnodes);

   D = AxisymD (element);
/*	  PrintMatrix(D,stdout);   */

   if (D == NullMatrix)
      return 1;
   
   for (i = 1 ; i <= ninteg ; i++) {
      if (VectorData (jac) [i] <= 0.0) {
         error ("det |J| for elt %d is <= 0, check elt distortion",element -> number);
         return 1;
      }
   } 

   element -> K = CreateMatrix (8,8);

   if (numnodes == 3) {
      MatrixRows (element -> K) = 6;
      MatrixCols (element -> K) = 6;
   }

   ZeroMatrix (element -> K);

	/*
	 * set-up so that multiplications work right
	 */

   MatrixRows (tempK) = 2*numnodes;
   MatrixCols (tempK) = 2*numnodes;

      xc1 = element -> node[1] -> x;
      xc3 = element -> node[3] -> x;

   for (i = 1 ; i <= ninteg ; i++) {
      B = IsoAQuadLocalB (element, numnodes,N, dNdx, dNdy, i,xpoints);
/*	  PrintMatrix(B,stdout);   */
      if (B == NullMatrix)
         return 1;

      MatrixRows (Bt) = MatrixRows (temp) = MatrixCols (B);
      TransposeMatrix (Bt, B);

      MultiplyMatrices (temp, Bt, D);
      MultiplyMatrices (tempK, temp, B);

      ScaleMatrix (tempK, tempK, VectorData (weights) [i]*VectorData (jac) [i], 0.0);
      sg  = VectorData (xpoints) [i];
	  maxx=xc3;
	  minx=xc1;
      rr=(minx+maxx)/2+sg*(maxx-minx)/2; 
      ScaleMatrix (tempK, tempK, rr*4*acos(0), 0.0);
      AddMatrices (element -> K, element -> K, tempK);
   }

	/*
	 * clean out the 7 & 8th rows and columns if this was a triangular
	 * element so nothing extra gets assembled into the global
	 * stiffness routines
	 */

   if (numnodes == 3) {
      for (i = 1 ; i <= 8 ; i++) 
         for (j = 7 ; j <= 8 ; j++) 
            MatrixData (element -> K) [i][j] = 0.0;

      for (i = 7 ; i <= 8 ; i++) 
         for (j = 1 ; j <= 6 ; j++) 
            MatrixData (element -> K) [i][j] = 0.0;
   }

   if (mass_mode) {
      element -> M = CreateMatrix (8, 8);
   
      if (mass_mode == 'l') 
         AQuadLumpedMassMatrix (element, numnodes);
      else if (mass_mode == 'c')
         AQuadConsistentMassMatrix (element, numnodes);
   }

   if (element -> numdistributed > 0) {
      equiv = IsoAQuadEquivNodalForces (element, &count);
      if (equiv == NullMatrix)
         return count;

       for (i = 1; i <= numnodes ; i++) {
          element -> node[i] -> eq_force[1] += VectorData (equiv) [2*i - 1];
          element -> node[i] -> eq_force[2] += VectorData (equiv) [2*i];
       }
   }

   return 0;
} 

int AQuadElementStress (element)
   Element	element;
{
   static Vector	stress = NullMatrix,
			d;
   static Matrix	temp;
   static Vector	weights,xpoints;
   static Matrix	N, dNdxi, dNde,
                        dNdx, dNdy = NullMatrix;
   unsigned		numnodes;
   int			ninteg;
   Matrix		D,
			B;
   double		diameter;
   double		sigma1,sigma2,
			theta,
			sigma_x,
			sigma_y,
			tau_xy;
   Vector		jac;
   unsigned		i,j;
   double		x,y;

   /* NOT IMPLEMENTED IN A WORKING WAY !!*/
   printf("Warning. This part does not work!!!\n");

   if (dNdy == NullMatrix) {
  
      N     = CreateMatrix (4,4);
      dNdxi = CreateMatrix (4,4);
      dNde  = CreateMatrix (4,4);
      dNdx  = CreateMatrix (4,4);
      dNdy  = CreateMatrix (4,4);
      weights = CreateVector (4);
      xpoints = CreateVector (8);
   }
   
   if (stress == NullMatrix) {
      stress = CreateVector (3);
      d = CreateVector (8);
      temp = CreateMatrix (3,8);
   }

   ninteg = 4;

      D = AxisymD (element);

   if (D == NullMatrix)
      return 1;

   if (element -> number == 1)
      numnodes = LocalAQuadShapeFunctions (element, ninteg, N, 
                                          dNdxi, dNde, weights,xpoints, FORCE);  
   else
      numnodes = LocalAQuadShapeFunctions (element, ninteg, N, 
                                          dNdxi, dNde, weights,xpoints, NOFORCE);  

   ninteg = 4;
   
   jac = GlobalAQuadShapeFunctions (element,dNdxi,dNde,dNdx,dNdy,
                                   ninteg,numnodes);

   for (i = 1 ; i <= numnodes ; i++) {
      VectorData (d) [2*i - 1] = element -> node[i] -> dx[1];
      VectorData (d) [2*i] = element -> node[i] -> dx[2];
   }

   MatrixRows (d) = numnodes*2;
   MatrixCols (temp) = numnodes*2;

   element -> ninteg = ninteg;
   SetupStressMemory (element);

   for (i = 1 ; i <= ninteg ; i++) {
      B = IsoAQuadLocalB (element, numnodes, N, dNdx, dNdy, i,xpoints);
      if (B == NullMatrix)
         return 1;

      x = y = 0.0;
      for (j = 1 ; j <= numnodes ; j++) {
         x += MatrixData (N)[j][i]*element -> node[j] -> x;
         y += MatrixData (N)[j][i]*element -> node[j] -> y;
      }
   
      MultiplyMatrices (temp, D, B);  

      MultiplyMatrices (stress, temp, d);
    
      sigma_x = VectorData (stress) [1];
      sigma_y = VectorData (stress) [2];
      tau_xy = VectorData (stress) [3];

      diameter = sqrt((sigma_x-sigma_y)*(sigma_x-sigma_y)/4 + tau_xy*tau_xy); 

      if (sigma_x - sigma_y != 0) {
         theta = 0.5*atan2(2*tau_xy,(sigma_x - sigma_y));
         if (sigma_x - sigma_y < 0) {
            if (tau_xy < 0)
               theta += M_PI_2;
            else
               theta -= M_PI_2;
         }
      }
      else
         theta = 0;

      sigma1 = (sigma_x + sigma_y)/2 + diameter;
      sigma2 = (sigma_x + sigma_y)/2 - diameter;

      element -> stress [i] -> x = x;
      element -> stress [i] -> y = y; 
      element -> stress [i] -> values [1] = sigma_x;
      element -> stress [i] -> values [2] = sigma_y;
      element -> stress [i] -> values [3] = tau_xy;
      element -> stress [i] -> values [4] = sigma1;
      element -> stress [i] -> values [5] = sigma2;
      element -> stress [i] -> values [6] = theta*180/M_PI;
   }

   return 0;
} 

void AQuadLumpedMassMatrix (element, numnodes)
   Element	element;
   unsigned	numnodes;
{
   double	area, mass;
   unsigned	i;

   area = ElementArea (element, numnodes);
   mass = area * element -> material -> rho * 
          element -> material -> t / (double) numnodes;

   ZeroMatrix (element -> M);

   for (i = 1 ; i <= numnodes ; i++) {
      MatrixData (element -> M) [2*i - 1][2*i - 1] = mass;
      MatrixData (element -> M) [2*i][2*i] = mass;
   }

   return;
}

void AQuadConsistentMassMatrix (element, numnodes)
   Element  element;
   unsigned numnodes;
{
   double xc1,xc3,rr,maxx,minx,sg;
   unsigned     i,j;
   int          ninteg;
   Vector       jac;
   Vector       equiv;
   int          count;
   Vector   weights,xpoints;
   Matrix   tempM;
   Matrix   N, dNdxi, dNde,
                        dNdx, dNdy;
   Matrix   Nt, temp,NN;



      N     = CreateMatrix (4,4);
      dNdxi = CreateMatrix (4,4);
      dNde  = CreateMatrix (4,4);
      dNdx  = CreateMatrix (4,4);
      dNdy  = CreateMatrix (4,4);
      weights = CreateVector (4);
      xpoints = CreateVector (8);
      tempM = CreateMatrix (8,8);
      Nt    = CreateMatrix (8,2);
      NN    = CreateMatrix (2,8);
      temp  = CreateMatrix (8,2);




   ninteg = 4;  /* 2 x 2 quadrature */

   numnodes = LocalAQuadShapeFunctions (element, ninteg, N,
                                       dNdxi, dNde, weights,xpoints, FORCE);

   jac = GlobalAQuadShapeFunctions (element,dNdxi,dNde,dNdx,dNdy,
                                   ninteg,numnodes);

   element -> M = CreateMatrix (8,8);
   ZeroMatrix (element -> M);

   MatrixRows (tempM) = 2*numnodes;
   MatrixCols (tempM) = 2*numnodes;

   xc1 = element -> node[1] -> x;
   xc3 = element -> node[3] -> x;
   maxx=xc3;
   minx=xc1;