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

loom.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:         loom.c
 *
 * Description:  Contains code for the driver application for the 
 *           FElt application loom, which is the command line
 *           based solution and post-processing engine intended
 *           mainly as the mathematical engine to be called from
 *           encapsulators such as WinFElt.
 *          
 *           It is not, nor is it meant to be, particularly 
 *           user-friendly.  It's meant to be called by other
 *           application using shell() or system() commands
 *           and has such has very mechanical, very overdone
 *           command-line specification which is probably built
 *           by the calling application based on check boxes
 *           and text entries ...
 *
 *****************************************************************************/

# include <stdio.h>
# include <string.h>
# include <unistd.h>
# include "problem.h"
# include "fe.h"
# include "error.h"
# include "allocate.h"
# include "definition.h"
# include "options.h"

typedef int Boolean;

static    FILE    *fp_err = NULL;
static    FILE    *fp_out = NULL;

static char *in = NULL;
static char *out = NULL;
static char *err_out = NULL;
static Boolean  unlink_input = 0;

      /*
       * Basic felt options
       */

static Boolean     solve = 0;
static Boolean   details = 0;
static Boolean   debug = 0;
static Boolean   no_err = 0;
static Boolean   summary = 0;
static Boolean   matrices = 0;
static Boolean   transfer = 0;
static Boolean   eigen = 0;
static Boolean   orthonormal = 0;
static Boolean   table = 0;
static Boolean   renumber = 0;

static char       *structure_out = NULL;
static char       *displaced_out = NULL;
static double   displaced_mag;
static double   xrot;
static double   yrot;
static double   zrot;
static double   zsc;
static int  contour_width = 0;
static int  contour_height = 0;

static char *graph_out = NULL;

static char    **contour_out = NULL;
static int       num_contour;
static char     *contour_type;
static char *contour_result;
static char *contour_format;
static int      *contour_component;
static Boolean     contour_hequal = 0;
static Boolean     contour_overlay = 0;

# define Check(n)         if (n < 0) return 1             /* simple error */
# define ReqCheck(n)      if (n < 1) return 1             /* required     */
# define OptCheck(n, opt) if (n < 1 && (opt)) return 1      /* req if opt   */
# define NeedCheck(n, opt) if (n == 1 && !(opt)) return 1   /* requires opt */
 
/************************************************************************
 * Function:      ParseOptions                                    *
 *                                                    *
 * Description:   Parses the specific command line options.       *
 ************************************************************************/

static int ParseOptions (argc, argv)
    int   argc;
    char *argv [ ];
{
    int           n;

    n = GetSoloStringOption(argc, argv, "in", &in);               ReqCheck(n);
    n = GetSoloStringOption(argc, argv, "out", &out);             ReqCheck(n);
    n = GetSoloStringOption(argc, argv, "err_out", &err_out);     ReqCheck(n);
    n = GetBooleanOption(argc, argv, "unlink", &unlink_input);    Check(n);
    
    n = GetBooleanOption(argc, argv, "solve",       &solve);      Check(n);
    n = GetBooleanOption(argc, argv, "no_err", &no_err); Check(n);
    n = GetBooleanOption(argc, argv, "debug",       &debug);      Check(n);
    n = GetBooleanOption(argc, argv, "details",     &details);    Check(n);
    n = GetBooleanOption(argc, argv, "summary",     &summary);    Check(n);
    n = GetBooleanOption(argc, argv, "matrices",    &matrices);   NeedCheck(n, solve);
    n = GetBooleanOption(argc, argv, "transfer",    &transfer);   NeedCheck(n, solve);
    n = GetBooleanOption(argc, argv, "eigen",       &eigen);      NeedCheck(n, solve);
    n = GetBooleanOption(argc, argv, "table",       &table);      NeedCheck(n, solve);
    n = GetBooleanOption(argc, argv, "orthonormal", &orthonormal);  NeedCheck(n, solve);
    n = GetBooleanOption(argc, argv, "renumber",    &renumber);   Check(n);
    
    n = GetSoloStringOption(argc, argv, "graph_out", &graph_out); 
    NeedCheck(n, solve);
   
    n = GetSoloStringOption(argc, argv, "structure_out", &structure_out);  
    Check(n);

    n = GetSoloStringOption(argc, argv, "displaced_out", &displaced_out);  
    NeedCheck(n, solve);
    n = GetSoloDoubleOption(argc, argv, "displaced_mag", &displaced_mag);  
    OptCheck(n, displaced_out);


    n = GetSoloDoubleOption(argc, argv, "xrot", &xrot);
    OptCheck(n, displaced_out || structure_out);
    n = GetSoloDoubleOption(argc, argv, "yrot", &yrot);
    OptCheck(n, displaced_out || structure_out);
    n = GetSoloDoubleOption(argc, argv, "zrot", &zrot);
    OptCheck(n, displaced_out || structure_out);
    n = GetSoloDoubleOption(argc, argv, "zscale", &zsc);
    OptCheck(n, displaced_out || structure_out);

    n = GetStringOption(argc, argv, "contour_out", &contour_out);
    NeedCheck(n, solve);
    num_contour = n;

    n = GetIntegerOption(argc, argv, "contour_component", &contour_component);
    OptCheck(n, contour_out);

    if (n != num_contour)
       return 1;
       
    n = GetBooleanOption(argc, argv, "contour_overlay", &contour_overlay);
    NeedCheck(n, contour_out);
    n = GetBooleanOption(argc, argv, "contour_hequal", &contour_hequal);
    NeedCheck(n, contour_out);
    n = GetSoloStringOption(argc, argv, "contour_result", &contour_result);
    OptCheck(n, contour_out);
    n = GetSoloStringOption(argc, argv, "contour_format", &contour_format);
    NeedCheck(n, contour_out);
    n = GetSoloStringOption(argc, argv, "contour_type", &contour_type);
    NeedCheck(n, contour_out);
    n = GetSoloIntegerOption(argc, argv, "contour_width", &contour_width);
    NeedCheck(n, contour_out);
    n = GetSoloIntegerOption(argc, argv, "contour_height", &contour_height);
    NeedCheck(n, contour_out);
    
    if (ArgsUsed ( ) + 1 < argc)
       return 1;

    return 0;
}

static void DumpOptions ( )
{
   printf ("in = %s\n", in);
   printf ("out = %s\n", out);
   printf ("err_out = %s\n\n", err_out);

   printf ("graph_out = %s\n", graph_out);
   printf ("contour_out = %s\n", contour_out);
   printf ("structure_out = %s\n", structure_out);
   printf ("displaced_out = %s\n\n", displaced_out);

   printf ("xrot = %g\n", xrot);
   printf ("yrot = %g\n", yrot);
   printf ("zrot = %g\n", zrot);
   printf ("zsc = %g\n", zsc);

   printf ("solve = %d\n\n", solve);

   printf ("debug = %d\n", debug);
   printf ("details = %d\n", details);
   printf ("summary = %d\n", summary);
   printf ("matrices = %d\n", matrices);
   printf ("transfer = %d\n", transfer);
   printf ("eigen = %d\n", eigen);
   printf ("orthonormal = %d\n", orthonormal);
   printf ("renumber = %d\n", renumber);
   printf ("table = %d\n", table);
}

/************************************************************************
 * Function:       ExitLoom                                 *
 *                                                    *
 * Description:                                             *
 ************************************************************************/

void ExitLoom(status)
    int      status;
{

    if (fp_out)
        fclose(fp_out);

    if (fp_err)
        fclose(fp_err);

    if (unlink_input)
       unlink(in); 

    exit(status);
}

/************************************************************************
 * Function:       main                                     *
 *                                                    *
 * Description:    Main is the driver function for loom.                *
 ************************************************************************/

int main (argc, argv)
    int   argc;
    char *argv [ ];
{
    int            i;
    char    *title;                 /* title of problem           */
    Matrix   M, K, C;         /* global matrices            */
    Matrix   Mcond, Ccond, Kcond;   /* condensed matrices         */
    Matrix   Mm, Km, Cm;            /* modal matrices       */
    Matrix   *H;              /* transfer function matrices   */
    Matrix   S;               /* output spectra       */
    Vector   F,               /* force vector               */
             Fcond;                 /* condensed force vector     */
    Vector   d;               /* displacement vector        */
    Matrix   x;               /* eigenvectors               */
    Vector   lambda;          /* eigenvalues                */
    NodeDOF  *forced;         /* array of forced DOF numbers  */
    unsigned       numforced;       /* number of forced DOF       */
    int            status;          /* return status        */
    Reaction       *R;              /* reaction force vector      */
    unsigned       numreactions;          /* the number of reactions    */
    Matrix   dtable;          /* time-displacement table    */
    unsigned       *old_numbers;          /* original node numbering    */
    AnalysisType mode;


      /*
       * if there are any option errors then do nothing 
       * at all because there must be a bug ...
       */

    if (ParseOptions (argc, argv)) {
        error ("options error");
        ExitLoom (1);
    }
/*
    DumpOptions ( );
*/
    if (in == NULL || out == NULL || err_out == NULL) {
        error ("file error");
        ExitLoom (1);
    }

    if (!no_err)
       fp_err = freopen(err_out, "w", stderr);
    else 
       fp_err = stderr;

      /*
       * we just have to bail if this didn't work because we don't
       * know where any error output should go
       */

    if (fp_err == NULL) {
        error ("file redirect error");
        ExitLoom (1);
    }

    fp_out = fopen(out, "w");
    if (fp_out == NULL)
       Fatal ("temporary file error -> output");

    add_all_definitions ( );

      /*
       * this will generate its own syntax errors so we don't need to
       * say anything
       */

    if (ReadFeltFile (in)) 
        ExitLoom (1);

    title    = problem.title;

      /*
       * If debugging write the problem as we understand it
       */

    if (debug) 
      fDumpFeltFile (fp_out);

    if (problem.num_nodes == 0 || problem.num_elements == 0) 
        Fatal ("nothing to do");

      /*
       * Write a graphics file if necessary
       */

    if (structure_out) 
       WriteWireframeFile (structure_out, 0.0, xrot, yrot, zrot, zsc);

      /*
       * compute material statistics if desired
       */
    
    if (summary)
       WriteMaterialStatistics (fp_out);

      /*
       * if we're not actually going to solve for anything
       * we can bail out right here
       */

    if (!solve)
       ExitLoom (0);

        /*
         * if the user wanted analysis details we need
         * to set the detail stream
         */

    if (details)
        SetDetailStream (fp_out);

      /*
       * find all the active DOFs in this problem     
       */

    FindDOFS ( );
       
      /*
       * renumber the nodes if desired
       */

    old_numbers = NULL;

    if (renumber) 
      old_numbers = RenumberNodes (problem.nodes, problem.elements, 
                                     problem.num_nodes, problem.num_elements);

      /*
       * switch on the problem type
       */

    mode = SetAnalysisMode ( );

    switch (mode) {

       case Transient:
          status = CheckAnalysisParameters (Transient);

          if (status) 
             Fatal ("%d Errors found in analysis parameters.", status);

          status = ConstructDynamic (&K, &M, &C);
          if (matrices)
             PrintGlobalMatrices (fp_out, M, C, K);
 
          if (status)
             Fatal ("%d fatal errors in stiffness and mass definitions",status);

          dtable = IntegrateHyperbolicDE (K, M, C);

          if (dtable == NullMatrix)
             Fatal ("fatal error in integration (probably a singularity).");

          if (table)
             WriteTransientTable (dtable, NullMatrix, fp_out);
    
          if (graph_out)
             WriteLineGraph (dtable, "Nodal Time-Displacement", "time", "dx", graph_out);

          break;
       
       case TransientThermal:
          analysis.dofs [1] = Tx;
          analysis.numdofs = 1;
          status = CheckAnalysisParameters (TransientThermal);

          if (status) 
             Fatal ("%d Errors found in analysis parameters.", status);

          status = ConstructDynamic (&K, &M, &C);
          
          if (matrices) 
             PrintGlobalMatrices (fp_out, M, NullMatrix, K);
 
          if (status)
             Fatal ("%d fatal errors in stiffness and mass definitions",status);

          dtable = IntegrateParabolicDE (K, M);

          if (dtable == NullMatrix)
             Fatal ("fatal error in integration (probably a singularity).");

          if (table)
             WriteTransientTable (dtable, NullMatrix, fp_out);
    
          if (graph_out)
             WriteLineGraph (dtable, "Nodal Time-Temperature", "time", "T", graph_out);

          break;

       case Static:

          K = ConstructStiffness(&status);
          if (status)
             Fatal ("%d Fatal errors in element stiffness definitions", status);

          if (matrices)
             PrintGlobalMatrices (fp_out, NullMatrix, NullMatrix, K);

          F = ConstructForceVector ( );
          
          ZeroConstrainedDOF (K, F, &Kcond, &Fcond);
     
          d = SolveForDisplacements (Kcond, Fcond);
          if (d == NullVector)
             Fatal("could not solve for displacements, probably a singularity");
    
          status = ElementStresses ( );

          if (status) 
             Fatal ("%d Fatal errors found computing element stresses", status);

          numreactions = SolveForReactions (K, d, old_numbers, &R);

          if (old_numbers != NULL)
             RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);

          if (table)
             WriteStructuralResults (fp_out, title, R, numreactions);

          if (displaced_out) 
             WriteWireframeFile (displaced_out, displaced_mag, 
                                 xrot, yrot, zrot, zsc);

          if (contour_out) {
             for (i = 0 ; i < num_contour ; i++) {
                if (strcmp(contour_result, "stress") == 0)
                   PlotStressField (contour_out [i], problem.elements, 
                            problem.num_elements, contour_component [i],
                            contour_hequal, contour_overlay,
                                    contour_width, contour_height);
                else if (strcmp(contour_result, "displacement") == 0)
                   PlotDisplacementField (contour_out [i], problem.nodes, 
                                problem.num_nodes, problem.elements, 
                                  problem.num_elements,
                                          contour_component [i], 
                                contour_hequal, contour_overlay,
                                          contour_width, contour_height);
             }
          }

          break;

       case StaticThermal:

          K = ConstructStiffness(&status);

          if (status)
             Fatal ("%d Fatal errors in element stiffness definitions", status);

          if (matrices)
             PrintGlobalMatrices (fp_out, NullMatrix, NullMatrix, K);

          F = ConstructForceVector ( );
          
          ZeroConstrainedDOF (K, F, &Kcond, &Fcond);
     
          d = SolveForDisplacements (Kcond, Fcond);
          if (d == NullVector)
             Fatal("could not solve for displacements, probably a singularity");

          if (old_numbers != NULL)
             RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);

          if (table)
             WriteTemperatureResults (fp_out, title);

          if (contour_out) {
             if (strcmp(contour_result, "displacement") == 0)
                PlotDisplacementField (contour_out [0], problem.nodes, 
                               problem.num_nodes, problem.elements, 
                               problem.num_elements, 1,
                                       contour_hequal, contour_overlay,
                                       contour_width, contour_height);
          }

          break;

       case StaticLoadCases:
       case StaticLoadRange:
          status = CheckAnalysisParameters (mode);
          if (status)
             Fatal ("%d errors found in analysis parameters.", status);

          K = ConstructStiffness(&status);
          if (status)
             Fatal ("%d Fatal errors in element stiffness definitions", status);

          if (matrices)
             PrintGlobalMatrices (fp_out, NullMatrix, NullMatrix, K);

          F = ConstructForceVector ( );
          
          ZeroConstrainedDOF (K, F, &Kcond, &Fcond);
           
          if (mode == StaticLoadCases)
             dtable = SolveStaticLoadCases (Kcond, Fcond);
          else
             dtable = SolveStaticLoadRange (Kcond, Fcond);

          if (dtable == NullMatrix)
             Fatal ("could not solve for global displacements");
            
          if (old_numbers != NULL)
             RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);

          if (table)
             if (mode == StaticLoadCases)
                WriteLoadCaseTable (dtable, fp_out);
             else 
                WriteLoadRangeTable (dtable, fp_out);

          if (graph_out && mode == StaticLoadRange)
             WriteLineGraph (dtable, "Displacement vs. Force Level", "force", "dx", graph_out);

          break; 

       case StaticSubstitutionLoadRange:
       case StaticIncrementalLoadRange:
          status = CheckAnalysisParameters (StaticSubstitution);
          status += CheckAnalysisParameters (StaticLoadRange);
          if (status) 
             Fatal ("%d errors found in analysis parameters.", status);

          K = CreateNonlinearStiffness (&status);
          if (status)
             Fatal ("could not create global stiffness matrix");
         
          F = ConstructForceVector ( );
 
          if (mode == StaticSubstitutionLoadRange)
             dtable = SolveNonlinearLoadRange (K, F, 0);
          else
             dtable = SolveNonlinearLoadRange (K, F, 0); /* should be 1*/
             
          if (dtable == NullMatrix)
             Fatal ("did not converge on a solution");
         
          if (old_numbers != NULL)
             RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);
                
          if (table)
             WriteLoadRangeTable (dtable, fp_out);

          if (graph_out)
             WriteLineGraph (dtable, "Displacement vs. Force Level", "force", "dx", graph_out);

          break;

       case StaticSubstitution:
       case StaticIncremental:
          status = CheckAnalysisParameters (mode);
          if (status) 
             Fatal ("%d errors found in analysis parameters.", status);

          K = CreateNonlinearStiffness (&status);
          if (status)
             Fatal ("could not create global stiffness matrix");
         
          F = ConstructForceVector ( );
 
          if (mode == StaticSubstitution)
             d = StaticNonlinearDisplacements (K, F, 0);
          else
             d = StaticNonlinearDisplacements (K, F, 0); /* should be 1 */
             
          if (d == NullMatrix)
             Fatal ("did not converge on a solution");
         
          if (old_numbers != NULL)
             RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);
                
          if (table)
             WriteStructuralResults (fp_out, title, NULL, 0);

          if (contour_out) {
             if (strcmp(contour_result, "displacement") == 0)
                for (i = 0 ; i < num_contour ; i++)
                   PlotDisplacementField (contour_out [i], problem.nodes, 
                                  problem.num_nodes, problem.elements, 
                                  problem.num_elements, 
                                contour_component [i],
                                          contour_hequal, contour_overlay,
                                          contour_width, contour_height);
          }

          break;

       case Modal:
          status = CheckAnalysisParameters (Modal);

          if (status) 
             Fatal ("%d Errors found in analysis parameters.", status);

          status = ConstructDynamic (&K, &M, &C);
          
          if (status)
             Fatal ("%d fatal errors in stiffness and mass definitions",status);

          RemoveConstrainedDOF (K, M, C, &Kcond, &Mcond, &Ccond);

          if (matrices)
             PrintGlobalMatrices (fp_out, Mcond, Ccond, Kcond);
 
          status = ComputeEigenModes (Kcond, Mcond, &lambda, &x);

          if (status == M_NOTPOSITIVEDEFINITE)
             Fatal ("coefficient matrix is not positive definite.");
          else if (status)
             Fatal ("could not compute eigenmodes (report status code %d).", status);

          NormalizeByFirst (x, x);
          WriteEigenResults (lambda, x, title, fp_out);
/*
          if (doplot)
             PlotModeShapes (x, fp_out);
*/           
          if (!eigen) {
             FormModalMatrices (x, Mcond, Ccond, Kcond, &Mm, &Cm, &Km, orthonormal);
             WriteModalResults (fp_out, Mm, Cm, Km, lambda);
          }
           
          break;

       case Spectral:
          status = CheckAnalysisParameters (Spectral);

          if (status) 
             Fatal ("%d Errors found in analysis parameters.", status);

          status = ConstructDynamic (&K, &M, &C);

          if (status)
             Fatal ("%d fatal errors in stiffness and mass definitions",status);
 
          ZeroConstrainedDOF (K, NullMatrix, &Kcond, NULL);
          ZeroConstrainedDOF (M, NullMatrix, &Mcond, NULL);
          ZeroConstrainedDOF (C, NullMatrix, &Ccond, NULL);

          if (matrices)
             PrintGlobalMatrices (fp_out, Mcond, Ccond, Kcond);

          FindForcedDOF (&forced, &numforced);

          H = ComputeTransferFunctions (Mcond, Ccond, Kcond, forced, numforced);

          if (transfer) {
             if (table)
                WriteTransferFunctions (H, forced, numforced, fp_out);
             if (graph_out)
                WriteLineGraphTransferFunctions (H, forced, numforced, graph_out);
          }
          else {
             S = ComputeOutputSpectra (H, forced, numforced);
             if (S == NullMatrix)
                break; 

             if (table)  
                WriteOutputSpectra (S, fp_out);
       
             if (graph_out)
                WriteLineGraph (S, "Output Power Spectra", "frequency", "S", graph_out);
          }

          break;
    }

    ExitLoom (0);
  
    return 0;
}

Generated by  Doxygen 1.6.0   Back to index