/*
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: fefunc.c *
* *
* Description: This file contains the function definitions for the *
* finite element intrinsic functions. *
************************************************************************/
# include <stdio.h>
# include <string.h>
# include "felt.h"
# include "trap.h"
# include "apply.h"
# include "debug.h"
# include "error.h"
# include "coerce.h"
# include "fefunc.h"
# include "execute.h"
# include "problem.h"
# include "allocate.h"
# include "definition.h"
# include "our-stdlib.h"
# include "pathsearch.h"
static struct {
Definition definition;
Function set_up;
Function stress;
} burlap_defs [32];
static int num_defs;
static int is_static PROTO ((void));
static int check_analysis PROTO ((char *));
static int is_permutation PROTO ((char *, Array));
static int element_set_up PROTO ((Element, int));
static int element_stress PROTO ((Element));
static Matrix check_matrix PROTO ((char *, Matrix));
static unsigned *check_dofs PROTO ((descriptor *, char *, unsigned *));
/* These are only in misc.h. */
extern double ElementArea PROTO ((Element, unsigned));
extern double ElementLength PROTO ((Element, unsigned));
/************************************************************************
* Function: is_static *
* *
* Description: Returns nonzero if the problem is a static problem, and *
* zero otherwise. *
************************************************************************/
static int is_static ( )
{
return problem.mode == Static || problem.mode == StaticThermal;
}
/************************************************************************
* Function: check_analysis *
* *
* Description: Returns nonzero if the analysis parameters are properly *
* initialized, and zero otherwise. *
************************************************************************/
static int check_analysis (func)
char *func;
{
if (CheckAnalysisParameters (problem.mode))
return 0;
if (problem.mode == TransientThermal) {
analysis.numdofs = 1;
analysis.dofs [1] = Tx;
}
if (analysis.numnodes * analysis.numdofs == 0) {
rterror ("missing analysis parameters");
return 1;
}
return 1;
}
/************************************************************************
* Function: is_permutation *
* *
* Description: Returns nonzero if the array is a permutation vector *
* suitable for node renumbering, otherwise zero. *
************************************************************************/
static int is_permutation (func, array)
char *func;
Array array;
{
int i;
int j;
int k;
int n;
int *ptr;
if ((n = array -> length) != problem.num_nodes) {
rterror ("illegal array length (%d) in %s()", n, func);
return 1;
}
ptr = (int *) array -> ptr;
for (i = 1; i <= n; i ++) {
k = ptr [i];
if (k < 1 || k > n) {
rterror ("illegal node number (%d) in %s()", k, func);
return 0;
}
for (j = i + 1; j <= n; j ++)
if (ptr [j] == k) {
rterror ("duplicate node number (%d) in %s()", k, func);
return 0;
}
}
return 1;
}
/************************************************************************
* Function: element_set_up *
* *
* Description: Calls the burlap code for the set-up function for the *
* specified element. The FElt routines expect to call *
* a C function for the element set-up step. In this case *
* a new element definition has been added and has a *
* burlap function as its set-up function. We simply *
* look up the definition in the private array and call *
* the approriate function with the arguments transformed *
* from C objects into burlap descriptors. *
************************************************************************/
static int element_set_up (element, mass_mode)
Element element;
int mass_mode;
{
int i;
int status;
descriptor *arg1;
descriptor *arg2;
descriptor *function;
descriptor *result;
status = 1;
for (i = 0; i < num_defs; i ++)
if (burlap_defs [i].definition == element -> definition) {
function = push ( );
arg1 = push ( );
arg2 = push ( );
D_Type (function) = T_Function;
D_Temp (function) = F_False;
D_Trapped (function) = F_False;
D_Function (function) = burlap_defs [i].set_up;
D_Type (arg1) = T_Element;
D_Temp (arg1) = F_False;
D_Trapped (arg1) = F_False;
arg1->u.ptr = &element;
D_Type (arg2) = T_Int;
D_Temp (arg2) = F_False;
D_Trapped (arg2) = F_False;
D_Int (arg2) = &mass_mode;
if (!function_call (function, 2)) {
result = pop ( );
result = CoerceData (deref (result), T_Int);
status = D_Type (result) == T_Int ? *D_Int (result) : 1;
RecycleData (result);
}
}
return status;
}
/************************************************************************
* Function: element_stress *
* *
* Description: Calls the burlap code for the stress function for the *
* specified element. The FElt routines expect to call *
* a C function for the element stress computation. In *
* this case a new element definition has been added and *
* has a burlap function as its stress function. We *
* simply look up the definition in the private array and *
* call the approriate function with the arguments *
* transformed from C objects into burlap descriptors. *
************************************************************************/
static int element_stress (element)
Element element;
{
int i;
int status;
descriptor *arg1;
descriptor *function;
descriptor *result;
status = 1;
for (i = 0; i < num_defs; i ++)
if (burlap_defs [i].definition == element -> definition) {
function = push ( );
arg1 = push ( );
D_Type (function) = T_Function;
D_Temp (function) = F_False;
D_Trapped (function) = F_False;
D_Function (function) = burlap_defs [i].stress;
D_Type (arg1) = T_Element;
D_Temp (arg1) = F_False;
D_Trapped (arg1) = F_False;
arg1->u.ptr = &element;
if (!function_call (function, 1)) {
result = top ( );
result = CoerceData (deref (result), T_Int);
status = D_Type (result) == T_Int ? *D_Int (result) : 1;
RecycleData (result);
}
}
return status;
}
/************************************************************************
* Function: check_matrix *
* *
* Description: Verifies that the specified matrix is the specified *
* size and that it is compact, or can be compacted. *
************************************************************************/
static Matrix check_matrix (func, a)
char *func;
Matrix a;
{
Matrix b;
unsigned size;
size = problem.num_nodes * problem.num_dofs;
if (Mrows (a) != size || Mcols (a) != size) {
MatrixError (func, a, NULL, M_SIZEMISMATCH, F_True);
return NULL;
}
if (IsCompact (a))
return a;
if (!(b = MakeCompactFromFull (a))) {
MatrixError (func, a, NULL, M_NOTCOMPACT, F_True);
return NULL;
}
return b;
}
/************************************************************************
* Function: check_dofs *
* *
* Description: Verifies that the specified descriptor is a proper *
* array of the degrees of freedom. *
************************************************************************/
static unsigned *check_dofs (d, function, num_dofs)
descriptor *d;
char *function;
unsigned *num_dofs;
{
unsigned i;
unsigned *dofs;
Array array;
int last;
array = D_Array (d);
if (array -> length != 6) {
rterror ("improper array length in %s()", function);
return NULL;
} else if (D_Type (array) != T_Int) {
rterror ("improper array type in %s()", function);
return NULL;
} else {
last = 0;
*num_dofs = 0;
dofs = (unsigned *) array -> ptr;
for (i = 1; i <= 6; i ++)
if (dofs [i] > 6) {
rterror ("illegal DOF in %s()", function);
return NULL;
} else if (dofs [i] && dofs [i] != last + 1) {
rterror ("illegal DOF array in %s()", function);
return NULL;
} else if (dofs [i]) {
last = dofs [i];
(*num_dofs) ++;
}
if (*num_dofs == 0) {
rterror ("no active dofs in %s()", function);
return NULL;
}
return dofs;
}
}
/************************************************************************
* Function: area_func *
* *
* Description: Pops and computes the area of the descriptor on the *
* top of the stack and places the result on the stack. *
* This function is only defined for planar elements. *
************************************************************************/
int area_func (n)
int n;
{
Element e;
descriptor *arg;
descriptor *result;
descriptor temp;
int type_error;
int status;
result = top ( );
temp = *result;
arg = &temp;
status = 0;
type_error = F_False;
arg = deref (arg);
switch (D_Type (arg)) {
case T_Element:
e = *D_Element (arg);
if (e -> definition -> shape == Planar) {
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = ElementArea (e, e -> definition -> shapenodes);
} else {
rterror ("area of non-planar element is undefined");
status = 1;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("area", arg, NULL, NULL, F_True);
RecycleData (arg);
d_printf ("area ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: felt_func *
* *
* Description: Pops the descriptor on the top of the stack and treats *
* it as the name of a FElt file to load. The path named *
* by the FELT_PATH environment variable is searched for *
* the file. The result is the full path name of the *
* loaded file. The following types are legal: *
* *
* felt (string) -> string (FElt file load) *
* *
* An attempt is first made to coerce the argument to a *
* string value. *
************************************************************************/
int felt_func (n)
int n;
{
char *name;
static char *path;
descriptor *arg;
descriptor *result;
descriptor temp;
int type_error;
result = top ( );
temp = *result;
arg = &temp;
type_error = F_False;
arg = CoerceData (deref (arg), T_String);
if (D_Type (arg) == T_Null) {
read_felt (NULL);
CreateData (result, NULL, NULL, T_Null);
} else if (D_Type (arg) == T_String) {
if (!path)
path = getenv ("FELT_PATH");
name = pathsearch (path, *D_String (arg), ".flt", F_True);
CreateData (result, NULL, NULL, T_String, strlen (name) + 1);
strcpy (*D_String (result), name);
if (read_felt (name)) {
rterror ("unable to read '%s'", name);
**D_String (result) = 0;
}
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("felt", arg, NULL, NULL, F_True);
RecycleData (arg);
return type_error == F_True;
}
/************************************************************************
* Function: length_func *
* *
* Description: Pops and computes the length of the descriptor on the *
* top of the stack and places the result on the stack. *
* The following types are legal: *
* *
* length (string) -> scalar (number of characters) *
* length (scalar) -> scalar (always one) *
* length (matrix) -> scalar (number of elements) *
* length (array) -> scalar (number of elements) *
* length (element) -> scalar (length of linear element) *
* *
* An attempt is first made to coerce the argument to a *
* double value if the argument is not a string. *
************************************************************************/
int length_func (n)
int n;
{
Matrix a;
Element e;
descriptor *arg;
descriptor *result;
descriptor temp;
int type_error;
int status;
result = top ( );
temp = *result;
arg = &temp;
arg = deref (arg);
type_error = F_False;
status = 0;
if (D_Type (arg) != T_String)
arg = CoerceData (arg, T_Double);
switch (D_Type (arg)) {
case T_String:
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = strlen (*D_String (arg));
break;
case T_Double:
D_Type (result) = T_Double;
D_Temp (result) = F_False;
D_Trapped (result) = F_False;
D_Double (result) = dbllit (1);
break;
case T_Matrix:
a = D_Matrix (arg);
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = Mrows (a) * Mcols (a);
break;
case T_Array:
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = D_Array (arg) -> length;
break;
case T_Element:
e = *D_Element (arg);
if (e -> definition -> shape == Linear) {
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = ElementLength (e, 3);
} else {
rterror ("length of non-linear element is undefined");
status = 1;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("length", arg, NULL, NULL, F_True);
RecycleData (arg);
d_printf ("length ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: volume_func *
* *
* Description: Not available yet. *
************************************************************************/
int volume_func (n)
int n;
{
rterror ("volume() function is not available");
return 1;
}
/************************************************************************
* Function: add_definition_func *
* *
* Description: Pops the descriptors on the top of the stack and uses *
* them as structure members of an element definition *
* structure and places the result of adding the *
* definition on the stack. *
************************************************************************/
int add_definition_func (n)
int n;
{
descriptor *arg1;
descriptor *arg2;
descriptor *arg3;
descriptor *arg4;
descriptor *arg5;
descriptor *arg6;
descriptor *arg7;
descriptor *arg8;
descriptor *arg9;
descriptor *result;
descriptor temp;
Definition definition;
int type_error;
int success;
int status;
unsigned i;
unsigned num_dofs;
unsigned *dofs;
arg9 = pop ( );
arg8 = pop ( );
arg7 = pop ( );
arg6 = pop ( );
arg5 = pop ( );
arg4 = pop ( );
arg3 = pop ( );
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_String);
arg2 = CoerceData (deref (arg2), T_Function);
arg3 = CoerceData (deref (arg3), T_Function);
arg4 = CoerceData (deref (arg4), T_Int);
arg5 = CoerceData (deref (arg5), T_Int);
arg6 = CoerceData (deref (arg6), T_Int);
arg7 = CoerceData (deref (arg7), T_Int);
arg9 = CoerceData (deref (arg9), T_Int);
arg8 = deref (arg8);
arg8 = CollapseMatrix (arg8);
arg8 = CoerceToArray (arg8, T_Int);
if (D_Type (arg1) != T_String)
type_error = F_True;
else if (D_Type (arg2) != T_Function)
type_error = F_True;
else if (D_Type (arg3) != T_Function)
type_error = F_True;
else if (D_Type (arg4) != T_Int)
type_error = F_True;
else if (D_Type (arg5) != T_Int)
type_error = F_True;
else if (D_Type (arg6) != T_Int)
type_error = F_True;
else if (D_Type (arg7) != T_Int)
type_error = F_True;
else if (D_Type (arg8) != T_Array)
type_error = F_True;
else if (!(dofs = check_dofs (arg8, "add_definition", &num_dofs)))
status = 1;
else {
if (!(definition = LookupDefinition (*D_String (arg1)))) {
definition = New (struct definition);
definition -> name = *D_String (arg1);
definition -> setup = element_set_up;
definition -> stress = element_stress;
definition -> shape = *D_Int (arg4);
definition -> numnodes = *D_Int (arg5);
definition -> shapenodes = *D_Int (arg6);
definition -> numstresses = *D_Int (arg7);
definition -> numdofs = num_dofs;
definition -> retainK = D_Type (arg9) != T_Null;
for (i = 1; i <= 6; i ++)
definition -> dofs [i] = dofs [i];
burlap_defs [num_defs].definition = definition;
burlap_defs [num_defs].set_up = CopyFunction (D_Function (arg2));
burlap_defs [num_defs].stress = CopyFunction (D_Function (arg3));
num_defs ++;
AddDefinition (definition);
success = F_True;
} else
success = F_False;
D_Type (result) = T_Double;
D_Temp (result) = F_False;
D_Trapped (result) = F_False;
D_Double (result) = dbllit (success ? 0 : 1);
}
if (type_error == F_True)
TypeError ("add_definition", NULL, NULL, NULL, F_True);
RecycleData (arg1);
RecycleData (arg2);
RecycleData (arg3);
RecycleData (arg4);
RecycleData (arg5);
RecycleData (arg6);
RecycleData (arg7);
RecycleData (arg8);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: remove_definition_func *
* *
* Description: Pops the descriptor from the top of the stack and *
* removes the element definition named by the descriptor. *
* If the definition is successfully removed then a zero *
* descriptor is pushed on the stack. Otherwise, a one *
* descriptor is pushed on the stack. The following types *
* are legal: *
* *
* remove_definition (string) -> double *
* *
* An attempt is first made to coerce the argument to a *
* string value. *
************************************************************************/
int remove_definition_func (n)
int n;
{
descriptor *arg;
descriptor *result;
descriptor temp;
Definition definition;
int type_error;
int success;
unsigned i;
result = top ( );
temp = *result;
arg = &temp;
type_error = F_False;
arg = CoerceData (deref (arg), T_String);
switch (D_Type (arg)) {
case T_String:
success = F_True;
if ((definition = LookupDefinition (*D_String (arg)))) {
for (i = 1; i <= problem.num_elements; i ++)
if (problem.elements [i] -> definition == definition) {
printf ("remove_definition: definition still in use\n");
success = F_False;
break;
}
if (success == F_True) {
RemoveDefinition (definition);
for (i = 0; i < num_defs; i ++)
if (burlap_defs [i].definition == definition) {
DestroyFunction (burlap_defs [i].set_up);
DestroyFunction (burlap_defs [i].stress);
burlap_defs [i] = burlap_defs [-- num_defs];
Delete (definition);
break;
}
}
}
D_Type (result) = T_Double;
D_Temp (result) = F_False;
D_Trapped (result) = F_False;
D_Double (result) = dbllit (success ? 0 : 1);
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("remove_definition", arg, NULL, NULL, F_True);
RecycleData (arg);
return type_error == F_True;
}
/************************************************************************
* Function: assemble_func *
* *
* Description: Interface to ConstructStiffness() & ConstructDynamic(). *
* The stiffness matrix is returned. The damping and mass *
* matrices may be returned by specifying them as the *
* parameters. The active DOFs must have previously been *
* computed. *
************************************************************************/
int assemble_func (n)
int n;
{
Matrix K;
Matrix M;
Matrix C;
descriptor *arg1;
descriptor *arg2;
descriptor *result;
descriptor temp;
int type_error;
int status;
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
if (!assignable (arg1)) {
RecycleData (arg1);
arg1 = NULL;
} else
arg1 = deref (arg1);
if (!assignable (arg2)) {
RecycleData (arg2);
arg2 = NULL;
} else
arg2 = deref (arg2);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in assemble()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in assemble()");
status = 1;
} else {
if (is_static ( )) {
if ((K = ConstructStiffness (&status))) {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = K;
} else
CreateData (result, NULL, NULL, T_Null);
} else {
if (!ConstructDynamic (&K, &M, &C)) {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = K;
if (!arg1 || AssignObject (arg1, T_Matrix, F_True, M))
DestroyMatrix (M);
if (!arg2 || AssignObject (arg2, T_Matrix, F_True, C))
DestroyMatrix (C);
} else
CreateData (result, NULL, NULL, T_Null);
}
}
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: clear_nodes_func *
* *
* Description: Interface to ClearNodes(). The return value is always *
* a null descriptor. *
************************************************************************/
int clear_nodes_func (n)
int n;
{
descriptor *result;
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in clear_nodes()");
return 1;
}
result = push ( );
CreateData (result, NULL, NULL, T_Null);
ClearNodes ( );
return 0;
}
/************************************************************************
* Function: compute_modes_func *
* *
* Description: Interface to ComputeEigenModes(). The active DOFs must *
* have previously been computed. The return value is the *
* vector of eigenvalues. The matrix of eigenvectors can *
* of retrieved by specifying it as the remaining *
* parameter. The following types are legal: *
* *
* compute_modes (matrix, matrix, ...) -> matrix *
* *
* An attempt is first made to coerce the first two *
* arguments to matrices. *
************************************************************************/
int compute_modes_func (n)
int n;
{
Matrix l;
Matrix x;
Matrix K;
Matrix M;
descriptor *arg1;
descriptor *arg2;
descriptor *arg3;
descriptor *result;
descriptor temp;
int type_error;
int status;
arg3 = pop ( );
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_Matrix);
arg2 = CoerceData (deref (arg2), T_Matrix);
if (!assignable (arg3)) {
RecycleData (arg3);
arg3 = NULL;
} else
arg3 = deref (arg3);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in compute_modes()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in compute_modes()");
status = 1;
} else if (D_Type (arg1) == T_Matrix && D_Type (arg2) == T_Matrix) {
K = D_Matrix (arg1);
M = D_Matrix (arg2);
if ((status = ComputeEigenModes (K, M, &l, &x)))
MatrixError ("compute_modes", M, NULL, status, F_True);
else {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = l;
if (!arg3 || AssignObject (arg3, T_Matrix, F_True, x))
DestroyMatrix (x);
}
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("compute_modes", arg1, arg2, arg3, F_True);
RecycleData (arg1);
RecycleData (arg2);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: compute_stresses_func *
* *
* Description: Not available yet. *
************************************************************************/
int compute_stresses_func (n)
int n;
{
rterror ("compute_stresses() function is not available");
return 1;
}
/************************************************************************
* Function: construct_forces_func *
* *
* Description: Interface to ConstructForceVector() and *
* AssembleTransientForce(). The active DOFs must have *
* previously been computed. The following types are *
* legal: *
* *
* construct_forces ( ) -> matrix (force vector) *
* construct_forces (double) -> matrix (force vector) *
* *
* An attempt is first made to coerce the argument to a *
* double value. *
************************************************************************/
int construct_forces_func (n)
int n;
{
double t;
Matrix v;
descriptor *arg;
descriptor *result;
descriptor temp;
unsigned nrows;
int type_error;
int status;
result = top ( );
temp = *result;
arg = &temp;
status = 0;
type_error = F_False;
arg = CoerceData (deref (arg), T_Double);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in construct_forces()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in construct_forces()");
status = 1;
} else {
if (is_static ( )) {
v = ConstructForceVector ( );
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = v;
} else if (D_Type (arg) == T_Double || D_Type (arg) == T_Null) {
nrows = problem.num_nodes * problem.num_dofs;
CreateData (result, NULL, NULL, T_Matrix, nrows, 1);
v = D_Matrix (result);
t = D_Type (arg) == T_Double ? *D_Double (arg) : 0;
AssembleTransientForce (t, v);
} else
type_error = F_True;
}
if (type_error == F_True)
TypeError ("construct_forces", arg, NULL, NULL, F_True);
RecycleData (arg);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: find_dofs_func *
* *
* Description: Interface to FindDOFS(). The return value is the *
* number of degrees of freedom. *
************************************************************************/
int find_dofs_func (n)
int n;
{
ste *s;
descriptor *dofs;
descriptor *result;
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in find_dofs()");
return 1;
}
FindDOFS ( );
s = st_lookup (&var_st, "dofs_num");
dofs = global (s -> idx);
D_Array (dofs) -> length = problem.num_dofs;
result = push ( );
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = problem.num_dofs;
return 0;
}
/************************************************************************
* Function: global_dof_func *
* *
* Description: Interface to GlobalDOF(). The active DOFs must have *
* previously been computed and the local DOF must be in *
* range. The following types are legal: *
* *
* global_dof (int, int) -> int (global DOF) *
* global_dof (node, int) -> int (global DOF) *
* *
* An attempt is first made to coerce both arguments to *
* integer values. *
************************************************************************/
int global_dof_func (n)
int n;
{
descriptor *arg1;
descriptor *arg2;
descriptor *result;
descriptor temp;
int type_error;
int status;
int dof;
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_Int);
arg2 = CoerceData (deref (arg2), T_Int);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in global_dof()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in global_dof()");
status = 1;
} else if (D_Type (arg2) == T_Int) {
dof = *D_Int (arg2);
if (dof < 1 || dof > 6) {
rterror ("illegal DOF in global_dof()");
status = 1;
} else if (D_Type (arg1) == T_Int) {
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = GlobalDOF (*D_Int (arg1), dof);
} else if (D_Type (arg1) == T_Node) {
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = GlobalDOF ((*D_Node (arg1)) -> number, dof);
} else
type_error = F_True;
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("global_dof", arg1, arg2, NULL, F_True);
RecycleData (arg1);
RecycleData (arg2);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: integrate_hyperbolic_func *
* *
* Description: Interface to IntegrateHyperbolicDE(). The active DOFs *
* must have previously been computed and the matrices *
* must be of the proper size and either be compact *
* matrices or be able to be compacted. Additionally, the *
* renumbering vector must be a proper permutation vector. *
* The following types are legal: *
* *
* integrate_hyperbolic (mtx, mtx, mtx) -> matrix *
* integrate_hyperbolic (mtx, mtx, mtx, array) -> matrix *
* *
* An attempt is first made to coerce the first three *
* arguments to matrices and the last argument to an array *
* of integers. *
************************************************************************/
int integrate_hyperbolic_func (n)
int n;
{
Matrix K;
Matrix M;
Matrix C;
Matrix D;
descriptor *arg1;
descriptor *arg2;
descriptor *arg3;
descriptor *result;
descriptor temp;
int type_error;
int status;
arg3 = pop ( );
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_Matrix);
arg2 = CoerceData (deref (arg2), T_Matrix);
arg3 = CoerceData (deref (arg3), T_Matrix);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in integrate_hyperbolic()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in integrate_hyperbolic()");
status = 1;
} else if (!check_analysis ("integrate_hyperbolic")) {
status = 1;
} else {
if (D_Type (arg1) != T_Matrix)
type_error = F_True;
else if (D_Type (arg2) != T_Matrix)
type_error = F_True;
else if (D_Type (arg3) != T_Matrix)
type_error = F_True;
}
if (type_error == F_False && status == 0) {
K = check_matrix ("integrate_hyperbolic", D_Matrix (arg1));
M = K ? check_matrix ("integrate_hyperbolic", D_Matrix (arg2)) : NULL;
C = K && M ? check_matrix ("integrate_hyperbolic", D_Matrix (arg3)) : NULL;
if (K && M && C) {
if ((D = IntegrateHyperbolicDE (K, M, C))) {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = D;
} else
status = 1;
} else
status = 1;
if (K && D_Matrix (arg1) != K)
DestroyMatrix (K);
if (M && D_Matrix (arg2) != M)
DestroyMatrix (M);
if (C && D_Matrix (arg3) != C)
DestroyMatrix (C);
}
if (type_error == F_True)
TypeError ("integrate_hyperbolic", arg1, arg2, arg3, F_True);
RecycleData (arg1);
RecycleData (arg2);
RecycleData (arg3);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: integrate_parabolic_func *
* *
* Description: Interface to IntegrateParabolicDE(). The active DOFs *
* must have previously been computed and the matrices *
* must be of the proper size and either be compact *
* matrices or be able to be compacted. Additionally, the *
* renumbering vector must be a proper permutation vector. *
* The following types are legal: *
* *
* integrate_parabolic (matrix, matrix) -> matrix *
* integrate_parabolic (matrix, matrix, array) -> matrix *
* *
* An attempt is first made to coerce the first two *
* arguments to matrices and the last argument to an array *
* of integers. *
************************************************************************/
int integrate_parabolic_func (n)
int n;
{
Matrix K;
Matrix M;
Matrix D;
descriptor *arg1;
descriptor *arg2;
descriptor *result;
descriptor temp;
int type_error;
int status;
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_Matrix);
arg2 = CoerceData (deref (arg2), T_Matrix);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in integrate_parabolic()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in integrate_parabolic()");
status = 1;
} else if (!check_analysis ("integrate_parabolic")) {
status = 1;
} else {
if (D_Type (arg1) != T_Matrix)
type_error = F_True;
else if (D_Type (arg2) != T_Matrix)
type_error = F_True;
}
if (type_error == F_False && status == 0) {
K = check_matrix ("integrate_parabolic", D_Matrix (arg1));
M = K ? check_matrix ("integrate_parabolic", D_Matrix (arg2)) : NULL;
if (K && M) {
if ((D = IntegrateParabolicDE (K, M))) {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = D;
} else
status = 1;
} else {
rterror ("compact matrices required in integrate_parabolic()");
status = 1;
}
if (K && D_Matrix (arg1) != K)
DestroyMatrix (K);
if (M && D_Matrix (arg2) != M)
DestroyMatrix (M);
}
if (type_error == F_True)
TypeError ("integrate_parabolic", arg1, arg2, NULL, F_True);
RecycleData (arg1);
RecycleData (arg2);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: local_dof_func *
* *
* Description: Interface to LocalDOF(). The active DOFs must have *
* previously been computed and the global DOF must be in *
* range. The number of the node is returned instead of *
* the node object since the nodes may have been *
* renumbered. The local DOF can be retrieved by *
* specifying it as the remaining parameter. The *
* following types are legal: *
* *
* local_dof (int) -> int (node number) *
* local_dof (int, int) -> int (node number) *
* *
* An attempt is first made to coerce the first argument *
* to an integer value. *
************************************************************************/
int local_dof_func (n)
int n;
{
descriptor *arg1;
descriptor *arg2;
descriptor *result;
descriptor temp;
unsigned node;
unsigned local_dof;
int global_dof;
int type_error;
int status;
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_Int);
if (!assignable (arg2)) {
RecycleData (arg2);
arg2 = NULL;
} else
arg2 = deref (arg2);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in local_dof()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in local_dof()");
status = 1;
} else if (D_Type (arg1) == T_Int) {
global_dof = *D_Int (arg1);
if (global_dof < 1 || global_dof > problem.num_nodes*problem.num_dofs) {
rterror ("illegal DOF in local_dof()");
status = 1;
} else {
LocalDOF (global_dof, &node, &local_dof);
CreateData (result, NULL, NULL, T_Double);
*D_Double (result) = node;
if (arg2)
status = AssignObject (arg2, T_Int, F_False, &local_dof);
}
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("local_dof", arg1, NULL, NULL, F_True);
RecycleData (arg1);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: remove_constrained_func *
* *
* Description: Interface to RemoveConstrainedMatrixDOF(). The active *
* DOFs must have previously been computed and the matrix *
* must be of the proper size and either be compact or be *
* able to be compacted. The following types are legal: *
* *
* remove_constrained (matrix) -> matrix *
* *
* An attempt is first made to coerce the argument to a *
* matrix. *
************************************************************************/
int remove_constrained_func (n)
int n;
{
Matrix a;
Matrix b;
descriptor *arg;
descriptor *result;
descriptor temp;
int type_error;
int status;
result = top ( );
temp = *result;
arg = &temp;
status = 0;
type_error = F_False;
arg = CoerceData (deref (arg), T_Matrix);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in remove_constrained()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in remove_constrained()");
status = 1;
} else if (D_Type (arg) == T_Matrix) {
if ((a = check_matrix ("remove_constrained", D_Matrix (arg)))) {
if ((b = RemoveConstrainedMatrixDOF (a))) {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = b;
} else
status = 1;
if (a != D_Matrix (arg))
DestroyMatrix (a);
} else
status = 1;
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("remove_constrained", arg, NULL, NULL, F_True);
RecycleData (arg);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: renumber_nodes_func *
* *
* Description: Interface to RenumberNodes(). The return value is the *
* permutation vector of node numbers. *
************************************************************************/
int renumber_nodes_func (n)
int n;
{
unsigned *v;
unsigned num_elts;
unsigned num_nodes;
descriptor *result;
int handler;
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in renumber_nodes()");
return 1;
}
num_elts = problem.num_elements;
num_nodes = problem.num_nodes;
v = RenumberNodes (problem.nodes, problem.elements, num_nodes, num_elts);
result = push ( );
handler = AddTrap (strict_assignment);
CreateData (result, NULL, NULL, T_Array, v, T_Int, num_nodes, handler);
D_Trapped (result) = AddTrap (array_assignment);
return 0;
}
/************************************************************************
* Function: restore_numbers_func *
* *
* Description: Interface to RestoreNodeNumbers(). The array must be *
* proper permutation vector. The following types *
* are legal: *
* *
* restore_numbers (array) -> null (renumber nodes) *
* *
* An attempt is first made to coerce the argument to an *
* array of integer values. *
************************************************************************/
int restore_numbers_func (n)
int n;
{
descriptor *arg;
descriptor *result;
descriptor temp;
unsigned *ptr;
int type_error;
int status;
result = top ( );
temp = *result;
arg = &temp;
status = 0;
type_error = F_False;
arg = CoerceToArray (deref (arg), T_Int);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in restore_numbers()");
status = 1;
} else if (D_Type (arg) == T_Array) {
if ((status = is_permutation ("restore_numbers", D_Array (arg)))) {
ptr = (unsigned *) D_Array (arg) -> ptr;
CreateData (result, NULL, NULL, T_Null);
RestoreNodeNumbers (problem.nodes, ptr, problem.num_nodes);
}
} else
type_error = F_True;
RecycleData (arg);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: set_up_func *
* *
* Description: Not available yet. *
************************************************************************/
int set_up_func (n)
int n;
{
rterror ("set_up() function is not available");
return 1;
}
/************************************************************************
* Function: solve_displacements_func *
* *
* Description: Interface to SolveForDisplacements(). The active DOFs *
* must have previously been computed and the matrices *
* must be of the proper size and either be compact *
* matrices of be able to be compacted. The following *
* types are legal: *
* *
* solve_displacements (matrix, matrix) -> matrix *
* *
* An attempt is first made to coerce both arguments to *
* matrices. *
************************************************************************/
int solve_displacements_func (n)
int n;
{
Matrix K;
Matrix f;
Matrix d;
descriptor *arg1;
descriptor *arg2;
descriptor *result;
descriptor temp;
unsigned size;
int type_error;
int status;
arg2 = pop ( );
result = top ( );
temp = *result;
arg1 = &temp;
status = 0;
type_error = F_False;
arg1 = CoerceData (deref (arg1), T_Matrix);
arg2 = CoerceData (deref (arg2), T_Matrix);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in solve_displacements()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in solve_displacements()");
status = 1;
} else if (D_Type (arg1) == T_Matrix && D_Type (arg2) == T_Matrix) {
if ((K = check_matrix ("solve_displacements", D_Matrix (arg1)))) {
f = D_Matrix (arg2);
size = problem.num_nodes * problem.num_dofs;
if (Mcols (f) != 1 || Mrows (f) != size) {
MatrixError ("solve_displacements", f, NULL, M_SIZEMISMATCH, F_True);
status = 1;
} else {
f = CreateCopyMatrix (f); /* f is overwritten with d */
if ((d = SolveForDisplacements (K, f))) {
D_Type (result) = T_Matrix;
D_Temp (result) = F_True;
D_Trapped (result) = F_False;
D_Matrix (result) = d;
} else {
DestroyMatrix (f);
status = 1;
}
}
if (K != D_Matrix (arg1))
DestroyMatrix (K);
} else
status = 1;
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("solve_displacements", arg1, arg2, NULL, F_False);
RecycleData (arg1);
RecycleData (arg2);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: zero_constrained_func *
* *
* Description: Interface to ZeroConstrainedMatrixDOF(). The active *
* DOFs must have previously been computed and the matrix *
* must be of the proper size. The following types are *
* legal: *
* *
* zero_constrained (col-vector) -> col_vector *
* zero_constrained (matrix) -> matrix *
* *
* An attempt is first made to coerce the argument to a *
* matrix. *
************************************************************************/
int zero_constrained_func (n)
int n;
{
Matrix a;
descriptor *arg;
descriptor *result;
descriptor temp;
unsigned size;
int type_error;
int status;
result = top ( );
temp = *result;
arg = &temp;
status = 0;
type_error = F_False;
arg = CoerceData (deref (arg), T_Matrix);
if (!problem.num_nodes) {
rterror ("no FElt problem loaded in zero_constrained()");
status = 1;
} else if (!problem.num_dofs) {
rterror ("no active DOFs in zero_constrained()");
status = 1;
} else if (D_Type (arg) == T_Matrix) {
a = D_Matrix (arg);
size = problem.num_nodes * problem.num_dofs;
if (IsColumnVector (a) && Mrows (a) == size) {
CreateData (result, arg, NULL, T_Matrix, size, 1);
ZeroConstrainedMatrixDOF (D_Matrix (result), a);
} else if (IsSymmetricMatrix (a) && Mrows (a) == size) {
CreateData (result, arg, NULL, T_Matrix, size, size);
ZeroConstrainedMatrixDOF (D_Matrix (result), a);
} else {
MatrixError ("zero_constrained", a, NULL, M_SIZEMISMATCH, F_True);
status = 1;
}
} else
type_error = F_True;
if (type_error == F_True)
TypeError ("zero_constrained", arg, NULL, NULL, F_True);
RecycleData (arg);
return type_error == F_True || status != 0;
}
syntax highlighted by Code2HTML, v. 0.9.1