Logo Search packages:      
Sourcecode: felt version File versions  Download package

fefunc.c

/*
    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;
          D_Element (arg1) = &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;
          D_Element (arg1) = &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;
}

Generated by  Doxygen 1.6.0   Back to index