pltgz4.c

C 开发的有限元软件

/*
Fri Mar 25 10:56:32 MET 1994

This file is for the FElt finite element analysis package.
Copyright (C) 1993,1994 Jason I. Gobat and Darren C. Atkinson

The PLTGZ4 element has to be checked by some benchmarks in the x-y plane.
It is the simplest possible triangular plate element I can provide.
It is integrated analytically. I did not check for errors after codeing.
It is provided for discussion only WITHOUT ANY WARRANTY.
 
Let me know if there is any interest for a more developed formulation.
(arbitary Poisson ratio, stress evaluation, arbitray location in 3D,
loads, anisitropic material, ...)

You should not remove the following headers in every further use related to
the PLTGZ4 element formulation.

Best wishes
Gerhard Zirwas

----------------------------------------------------------------------------
Lehrstuhl fuer Baumechanik, TU-Muenchen, Arcisstr. 21, D-80333 Muenchen

Dipl.-Ing. Gerhard Zirwas

priv to    Schwojerstr. 68c, D-82140 Olching
mail to    [e-mail]
           or [e-mail]
phone to   [phone] or [phone] (private)
FAX to     [phone]
-------------------------------------- Confitemini domino quoniam bonus. ---
*/

/************************************************************************
 * File:	pltgz4.c						*
 *									*
 * Description: This file contains the definition structure and the	*
 *		stiffness function for an PLTGZ4			*
 *									*
 *		PLate elemet with 					*
 *		Triangular shape, where energy due to the shearmodulus	*
 *		G is 							*	
 *		Zeroed,	and the poisson ratio nu is simplified to 	*
 *		1/4							*
 *									*
 *		The energy is zero obviously only for the shears	*
 *  		involving the direction normal to the plate,		*
 * 		by locking the normals to the plate as normals 		*
 * 		after deformation. (Kirchhoff plate theory).		*
 *		The poisson ratio is simplified to nu=1/4 by		*
 *		setting the Lame coefficiant lambda=mu. (mu=G)		*
 ************************************************************************/

/* To add the element to Felt look into the Users Guide and Reference Manual
at 9.5 "Putting it all together"
The procedure is something like this:

copy this file to ../../FElt-2.0/lib/Elements/pltgz4.c

Add  PLTGZ4 to    ../../FElt-2.0/lib/Elements/elements.db

Do an make element.h in ../../FElt-2.0/lib/Elements

Add pltgz4.o to the Makefile in ../../FElt-2.0/lib/Elements which will then
look like:

OBJS	= misc.o beam.o beam3d.o cst.o truss.o iso_2d.o iso_quad.o\
          timoshenko.o pltgz4.o 

Go to the ../../FElt-2.0 and do a

make clean

and

make all

Now you should be able to run Felt with the plate element.
*/

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

void    PLTGZ4LumpedMassMatrix ( );
Matrix  PLTGZ4LocalK 	    ( );
int	PLTGZ4ElementSetup ( );
int	PLTGZ4ElementStress    ( );
void	PLTGZ4mpl ( );

struct definition PLTGZ4Definition = {
    "PLTGZ4",
    PLTGZ4ElementSetup,
    PLTGZ4ElementStress,
    Planar, 			/* The shape of this element		   	*/
    3, 				/* 3 nodes per element			   	*/
    3, 				/* 3 nodes define the shape (triangl)	   	*/
    0, 				/* 0,6 magnitudes in each stress structure	*/
    3, 				/* 3 global DOF / node			   	*/
   {0, 3, 4, 5, 0, 0, 0},      	/* DOF 1 is Tz, DOF 2 is Rx DOF 3 is Ry .. 	*/
    0				/* retain stiffness after assembling	   	*/
};

int PLTGZ4ElementSetup (element, mass_mode)
   Element	element;
   char		mass_mode;
{
   unsigned		i;
   Vector		equiv;
   int			count;
   Matrix		Kb;
   double		factor;
   double		area;

   if (element -> material -> nu != 0.25) {
      error ("PLTGZ4 element %d has Poisson's ratio (nu=1/4=0.25)", element -> number);
      return 1;
   }
   if (element -> material -> E == 0) {
      error ("PLTGZ4 element %d has 0.0 for material modulus E  ", element -> number);
      return 1;
   }
   if (element -> material -> t == 0) {
      error ("PLTGZ4 element %d has 0.0 for thickness (t)", element -> number);
      return 1;
   }

   Kb = PLTGZ4LocalK (element, &area);
   if (Kb == NullMatrix)
      return 1;
   
   element -> K = CreateMatrix (9,9);
   if (element -> K == NullMatrix) 
      Fatal ("allocation error creating element %d stiffness matrix", element -> number);

   factor = 1.0;
   ScaleMatrix (element -> K, Kb, factor, 0.0);

   if (element -> K == NullMatrix)
      Fatal ("element %d K matrix is null after definition",element -> number);

	/*	
	 * form the element mass matrix if necessary (note that we only
	 * have a lumped formulation for now)
	 */

   if (mass_mode) {
      element -> M = CreateMatrix (9,9);
      if (element -> M == NullMatrix)
         AllocationError (element, "mass matrix");

      if (mass_mode == 'l')
         PLTGZ4LumpedMassMatrix (element, area);
      else
         PLTGZ4LumpedMassMatrix (element, area);
   }

   return 0;
}
/*--------------------------------------------------------------------------------*/
int PLTGZ4ElementStress (element)
   Element	element;
{
   static Vector	stress = NullMatrix,
			d;
   
   if (stress == NullMatrix) {
      stress = CreateVector (3);
      d = CreateVector (9);

      if (stress == NullMatrix || d == NullMatrix )
         Fatal ("allocation error in element %d stresses", element -> number);
   }

/* stress evaluation will be plugged in. */

   return 0;
} 
/*--------------------------------------------------------------------------------*/
void PLTGZ4LumpedMassMatrix (element, area)
   Element	element;
   double	area;
{
   double	factor;
   unsigned	i;

   ZeroMatrix (element -> M); 

   factor = (element -> material -> t * element -> material -> rho * area)/3.0;

      MatrixData (element -> M) [1][1] = factor;   
      MatrixData (element -> M) [4][4] = factor;   
      MatrixData (element -> M) [7][7] = factor;   
   
   return;
}
/*--------------------------------------------------------------------------------*/
Matrix PLTGZ4LocalK (element, area)
   Element	element;
   double	*area;
{
   static Matrix 	Kb = NullMatrix;
   double		xc1,yc1,
			xc2,yc2,
			xc3,yc3,
			A,
			factor;
   double		xb,yb,xc,yc,h,		/* shifted coordinates; thickness	*/
			Em,lambda,mu,n;		/* constants of lin. elast. material 	*/
   double  		K_mpl[10][10];		/* Matrix coded by MapleV 		*/
   unsigned		i1,i2;

   if (Kb == NullMatrix) {
      Kb = CreateMatrix (9,9);

      if (Kb == NullMatrix)
         Fatal ("allocation error creating element stiffness matrix");
   }

   ZeroMatrix (Kb);

   h = element -> material -> t;

   xc1 = element -> node[1] -> x;
   xc2 = element -> node[2] -> x;
   xc3 = element -> node[3] -> x;
   yc1 = element -> node[1] -> y;
   yc2 = element -> node[2] -> y;
   yc3 = element -> node[3] -> y;
 
   xb=xc2-xc1;
   yb=yc2-yc1;
   xc=xc3-xc1;
   yc=yc3-yc1;
 
   A = 0.5*(xb*yc-xc*yb);
   
   if (A < 0) {
      error("incorrect node ordering for element %d (must be ccw)",element -> number);
      return NullMatrix;
   }
   if (A == 0) {
      error ("area of element %d is zero, check node numbering",element -> number);
      return NullMatrix;
   }

   if (area != NULL)
      (*area) = A;
   
   PLTGZ4mpl(xb,yb,xc,yc,h,&K_mpl[0][0]);
    
   for (i1 = 1 ; i1 <= 9 ; i1++) {
     for (i2 = 1 ; i2 <= 9 ; i2++) {
       MatrixData (Kb) [i1][i2] = K_mpl[i1][i2];
     }
   }

   Em = element -> material -> E;
   n  = element -> material -> nu;