/*
This file is part of the FElt finite element analysis package.
Copyright (C) 1993-2000 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_2d.c
*
* Description: contains the element definition routines for isoparametric
* plane stress / plane strain elements
*
******************************************************************************/
# include <stdio.h>
# include <math.h>
# include "fe.h"
# include "error.h"
# include "misc.h"
# define PLANESTRESS 1
# define PLANESTRAIN 2
extern int iso2d_PlaneStressEltSetup ( ), iso2d_PlaneStressEltStress ( );
extern int iso2d_PlaneStrainEltSetup ( ), iso2d_PlaneStrainEltStress ( );
struct definition iso2d_PlaneStressDefinition = {
"iso2d_PlaneStress",
iso2d_PlaneStressEltSetup, iso2d_PlaneStressEltStress,
Planar, 9, 4, 6, 2, {0, 1, 2, 0, 0, 0, 0}, 0
};
struct definition iso2d_PlaneStrainDefinition = {
"iso2d_PlaneStrain",
iso2d_PlaneStrainEltSetup, iso2d_PlaneStrainEltStress,
Planar, 9, 4, 6, 2, {0, 1, 2, 0, 0, 0, 0}, 0
};
unsigned LocalIsoShapeFunctions ( );
Vector GlobalIsoShapeFunctions ( );
Matrix Iso2dLocalB ( );
int iso2dElementSetup ( );
int iso2dElementStress ( );
int iso2d_PlaneStressEltSetup (element, mass_mode, tangent)
Element element;
char mass_mode;
int tangent;
{
return iso2dElementSetup (element, mass_mode, tangent, PLANESTRESS);
}
int iso2d_PlaneStrainEltSetup (element, mass_mode, tangent)
Element element;
char mass_mode;
int tangent;
{
return iso2dElementSetup (element, mass_mode, tangent, PLANESTRAIN);
}
int iso2d_PlaneStressEltStress (element)
Element element;
{
return iso2dElementStress (element, PLANESTRESS);
}
int iso2d_PlaneStrainEltStress (element)
Element element;
{
return iso2dElementStress (element, PLANESTRAIN);
}
int iso2dElementSetup (element, mass_mode, tangent, type)
Element element;
char mass_mode;
int tangent;
unsigned type;
{
unsigned numnodes;
unsigned i,j;
int ninteg;
Matrix B;
Matrix D;
Vector jac;
static Vector weights;
static Matrix tempK;
static Matrix N, dNdxi, dNde,
dNdx, dNdy = NullMatrix;
static Matrix Bt, temp;
if (dNdy == NullMatrix) {
N = CreateMatrix (9,9);
dNdxi = CreateMatrix (9,9);
dNde = CreateMatrix (9,9);
dNdx = CreateMatrix (9,9);
dNdy = CreateMatrix (9,9);
weights = CreateVector (9);
tempK = CreateMatrix (18,18);
Bt = CreateMatrix (18,3);
temp = CreateMatrix (18,3);
}
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;
}
if (element -> material -> t == 0) {
error ("isoparametric element %d has 0.0 for thickness (t)",element -> number);
return 1;
}
numnodes = LocalIsoShapeFunctions (element, N, dNdxi, dNde, weights);
if (element -> node[3] -> number == element -> node[4] -> number) {
if (numnodes != 4) {
error ("triangular isoparametric elt %d has more than 3 unique nodes",
element -> number);
return 1;
}
else
numnodes = 3;
}
if (numnodes > 4) {
ninteg = 9;
numnodes = 9;
} else
ninteg = 4;
ninteg = 4;
jac = GlobalIsoShapeFunctions (element,N,dNdxi,dNde,dNdx,dNdy,
ninteg,numnodes);
if (type == PLANESTRESS)
D = PlaneStressD (element);
else if (type == PLANESTRAIN)
D = PlaneStrainD (element);
else
D = NullMatrix; /* gcc -Wall */
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;
}
}
if (element -> K == NullMatrix) {
if (numnodes == 3) {
element -> K = CreateMatrix (8,8);
MatrixRows (element -> K) = 6;
MatrixCols (element -> K) = 6;
}
else
element -> K = CreateMatrix (2*numnodes, 2*numnodes);
}
ZeroMatrix (element -> K);
/*
* We need to fool the matrix routines since really we can have
* a variable size stiffness array ... there's no need to operate
* on all those zeros if we don't have to
*/
MatrixRows (tempK) = 2*numnodes;
MatrixCols (tempK) = 2*numnodes;
for (i = 1 ; i <= ninteg ; i++) {
B = Iso2dLocalB (element, numnodes, dNdx, dNdy, i);
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);
AddMatrices (element -> K, element -> K, tempK);
}
ScaleMatrix (element -> K, element -> K, element -> material -> t, 0.0);
/*
* 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;
}
return 0;
}
int iso2dElementStress (element, type)
Element element;
unsigned type;
{
element -> ninteg = 0;
return 0;
}
Matrix Iso2dLocalB (element, numnodes, dNdx, dNdy, point)
Element element;
unsigned numnodes;
Matrix dNdx, dNdy;
unsigned point;
{
unsigned i;
static Matrix B = NullMatrix;
if (B == NullMatrix)
B = CreateMatrix (3,18);
for (i = 1 ; i <= numnodes ; i++) {
MatrixData (B) [1][2*i - 1] = MatrixData (dNdx) [i][point];
MatrixData (B) [1][2*i] = 0.0;
MatrixData (B) [2][2*i - 1] = 0.0;
MatrixData (B) [2][2*i] = MatrixData (dNdy) [i][point];
MatrixData (B) [3][2*i - 1] = MatrixData (dNdy) [i][point];
MatrixData (B) [3][2*i] = MatrixData (dNdx) [i][point];
}
MatrixCols (B) = numnodes*2;
return B;
}
Vector GlobalIsoShapeFunctions (element, N, dNdxi, dNde, dNdx, dNdy,
ninteg,nodes)
Element element;
int ninteg;
Matrix N,
dNdxi, dNde,
dNdx, dNdy;
unsigned nodes;
{
unsigned i,j;
static Vector jac,
dxdxi, dxde,
dydxi, dyde = NullMatrix;
if (dyde == NullMatrix) {
dxdxi = CreateVector (9);
dxde = CreateVector (9);
dydxi = CreateVector (9);
dyde = CreateVector (9);
jac = CreateVector (9);
}
for (i = 1 ; i <= 9 ; i++) {
VectorData (dxdxi) [i] = 0.0;
VectorData (dxde) [i] = 0.0;
VectorData (dydxi) [i] = 0.0;
VectorData (dyde) [i] = 0.0;
VectorData (jac) [i] = 0.0;
}
for (i = 1 ; i <= ninteg ; i++) {
for (j = 1 ; j <= nodes ; j++) {
if (element -> node[j] != NULL) {
VectorData (dxdxi) [i] += MatrixData (dNdxi) [j][i]*
(element -> node[j] -> x);
VectorData (dxde) [i] += MatrixData (dNde) [j][i]*
(element -> node[j] -> x);
VectorData (dydxi) [i] += MatrixData (dNdxi) [j][i]*
(element -> node[j] -> y);
VectorData (dyde) [i] += MatrixData (dNde) [j][i]*
(element -> node[j] -> y);
}
}
VectorData (jac) [i] = VectorData (dxdxi)[i] * VectorData (dyde)[i] -
VectorData (dxde)[i] * VectorData (dydxi)[i];
for (j = 1 ; j <= nodes ; j++) {
if (element -> node[j] != NULL) {
MatrixData (dNdx) [j][i] = (MatrixData (dNdxi) [j][i]*VectorData (dyde)[i] -
MatrixData (dNde)[j][i]*VectorData (dydxi)[i])/
VectorData (jac)[i];
MatrixData (dNdy) [j][i] = -(MatrixData (dNdxi)[j][i]*VectorData (dxde)[i] -
MatrixData (dNde)[j][i]*VectorData (dxdxi)[i])/
VectorData (jac)[i];
}
else {
MatrixData (dNdx) [j][i] = 0.0;
MatrixData (dNdy) [j][i] = 0.0;
}
}
}
return jac;
}
/*****************************************************************************
*
* Function: LocalIsoShapeFunctions
*
* Description: calculates the shape functions and the derivatives (w/ respect
* to xi, eta coordinates) of the shape function for a four to
* nine node plane stress / plane strain element
*
* Note: The approach looks rather brutish, but it seems much clearer
* to me this way and we have to do each individual computation
* anyways, whether we put ourselves in a loop and use index
* notation or not ...
*
******************************************************************************/
unsigned LocalIsoShapeFunctions (element, N, dNdx, dNde, weights)
Element element;
Matrix N,
dNdx,
dNde;
Vector weights;
{
unsigned i,j,k;
int ninteg;
double eta,xi;
double Nt[10];
double de[10],dx[10];
double *gauss_points;
double *gauss_wts;
static unsigned numnodes;
static int points [10];
static int prev_points [10] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
unsigned same_flag;
same_flag = 1;
numnodes = 4;
for (i = 5 ; i <= 9 ; i++) {
points [i] = 0;
if (element -> node[i] != NULL) {
numnodes ++;
points [i] = 1;
}
if (points [i] != prev_points [i])
same_flag = 0;
}
if (same_flag)
return numnodes;
if (numnodes == 4)
ninteg = 2;
else
ninteg = 3;
for (i = 5 ; i <= 9 ; i++)
prev_points [i] = points [i];
ninteg = 2;
GaussPoints (ninteg, &gauss_points, &gauss_wts);
for (i = 0 ; i < ninteg ; i++) {
xi = gauss_points [i];
for (j = 0 ; j < ninteg ; j++) {
eta = gauss_points [j];
Nt [1] = 0.25*(1 - eta)*(1 - xi);
dx [1] = 0.25*(-1 + eta);
de [1] = 0.25*(-1 + xi);
Nt [2] = 0.25*(1 - eta)*(1 + xi);
dx [2] = 0.25*(1 - eta);
de [2] = 0.25*(-1 - xi);
Nt [3] = 0.25*(1 + eta)*(1 + xi);
dx [3] = 0.25*(1 + eta);
de [3] = 0.25*(1 + xi);
Nt [4] = 0.25*(1 + eta)*(1 - xi);
dx [4] = 0.25*(-1 - eta);
de [4] = 0.25*(1 - xi);
if (numnodes == 4 &&
element -> node[3] -> number == element -> node[4] -> number) {
Nt [3] += Nt [4];
dx [3] += dx [4];
de [3] += de [4];
/*
printf ("%g %g %g %g %g\n", eta, xi, Nt [3], dx [3], de [3]);
*/
}
if (numnodes > 4) {
if (points [9]) {
Nt [9] = (1 - eta*eta)*(1 - xi*xi);
dx [9] = 2*xi*(eta*eta - 1);
de [9] = 2*eta*(xi*xi - 1);
for (k = 1 ; k <= 4 ; k++) {
Nt [k] -= 0.25*Nt [9];
dx [k] -= 0.25*dx [9];
de [k] -= 0.25*de [9];
}
}
else
Nt [9] = dx [9] = de [9] = 0.0;
if (points [5]) {
Nt [5] = 0.5*((1 - eta)*(1 - xi*xi) - Nt[9]);
dx [5] = -xi + xi*eta - 0.5*dx[9];
de [5] = 0.5*(-1 + xi*xi - de[9]);
}
else
Nt [5] = dx [5] = de [5] = 0.0;
if (points [6]) {
Nt [6] = 0.5*((1 - eta*eta)*(1 + xi) - Nt[9]);
dx [6] = 0.5*(1 - eta*eta - dx[9]);
de [6] = -eta - eta*xi - 0.5*de[9];
}
else
Nt [6] = dx [6] = de [6] = 0.0;
if (points [7]) {
Nt [7] = 0.5*((1 + eta)*(1 - xi*xi) - Nt[9]);
dx [7] = -xi - eta*xi - 0.5*dx[9];
de [7] = 0.5*(1 - xi*xi - de[9]);
}
else
Nt [7] = dx [7] = de [7] = 0.0;
if (points [8]) {
Nt [8] = 0.5*((1 - eta*eta)*(1 - xi) - Nt[9]);
dx [8] = 0.5*(-1 + eta*eta - dx[9]);
de [8] = -eta + xi*eta - 0.5*de[9];
}
else
Nt [8] = dx [8] = de [8] = 0.0;
Nt [1] -= 0.5*(Nt[5] + Nt[8]);
dx [1] -= 0.5*(dx[5] + dx[8]);
de [1] -= 0.5*(de[5] + de[8]);
Nt [2] -= 0.5*(Nt[5] + Nt[6]);
dx [2] -= 0.5*(dx[5] + dx[6]);
de [2] -= 0.5*(de[5] + de[6]);
Nt [3] -= 0.5*(Nt[6] + Nt[7]);
dx [3] -= 0.5*(dx[6] + dx[7]);
de [3] -= 0.5*(de[6] + de[7]);
Nt [4] -= 0.5*(Nt[7] + Nt[8]);
dx [4] -= 0.5*(dx[7] + dx[8]);
de [4] -= 0.5*(de[7] + de[8]);
}
for (k = 1 ; k <= 9 ; k++) {
MatrixData (N) [k][i*ninteg + j+1] = Nt [k];
MatrixData (dNdx) [k][i*ninteg + j+1] = dx [k];
MatrixData (dNde) [k][i*ninteg + j+1] = de [k];
}
VectorData (weights) [i*ninteg + j+1] = gauss_wts [j];
}
}
return numnodes;
}
syntax highlighted by Code2HTML, v. 0.9.1