/*
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: felt.c
*
* Description: Contains code for the driver application for the Finite
* ELemenT (FELT) package currently under development.
*
*****************************************************************************/
# include <stdio.h>
# include <string.h>
# include "problem.h"
# include "fe.h"
# include "error.h"
# include "allocate.h"
# include "version.h"
# include "definition.h"
# define streq(a,b) !strcmp(a,b)
# ifndef DOS
static char *usage = "\
usage: felt [options] [filename]\n\
-debug write debugging output\n\
-preview write a simple ASCII rendering of the problem\n\
+table do not print tabular dynamic results\n\
-plot plot dynamic, load case, or load range results\n\
-transfer only show transfer functions for spectral results\n\
-eigen only compute eigen results for modal analysis\n\
-orthonormal use orthonormal mode shapes for modal matrices\n\
-renumber force automatic node renumbering\n\
-summary include material summary statistics\n\
-matrices print the global matrices\n\
-details print ancillary analysis details\n\
-matlab filename write the global matrices to a file\n\
-graphics filename create file for structure visualization\n\
-version print version information and exit\n\
-nocpp do not use a preprocessor\n\
-cpp filename preprocessor to use\n\
-Dname[=value] define a macro\n\
-Uname undefine a macro\n\
-Idirectory specify include directory\n\
";
# else
static char *usage = "\
usage: felt [options] [filename]\n\
-debug write debugging output\n\
-preview write a simple ASCII rendering of the problem\n\
+table do not print tabular dynamic results\n\
-plot do graphical plot for dynamic results\n\
-transfer only show transfer functions for spectral results\n\
-eigen only compute eigen results for modal analysis\n\
-orthonormal use orthonormal mode shapes for modal matrices\n\
-renumber force automatic node renumbering\n\
-summary include material summary statistics\n\
-matrices print the global matrices\n\
-details print ancillary analysis details\n\
-matlab filename write the global matrices to a file\n\
-graphics filename create file for structure visualization\n\
-version print version information and exit\n\
";
# endif
static int debug = 0;
static int preview = 0;
static int summary = 0;
static int matrices = 0;
static int dospectra = 1;
static int domodal = 1;
static int orthonormal = 0;
static int doplot = 0;
static int dotable = 1;
static int renumber = 0;
static int details = 0;
static char *graphics = NULL;
static char *matlab = NULL;
/************************************************************************
* Function: ParseFeltOptions *
* *
* Description: Parses the felt specific command line options. *
************************************************************************/
static int ParseFeltOptions (argc, argv)
int *argc;
char *argv [ ];
{
int i;
int j;
char *arg;
j = 1;
for (i = 1; i < *argc; i ++)
if (streq ((arg = argv [i]), "-help")) {
fputs (usage, stderr);
exit (0);
} else if (streq (arg, "-version")) {
fprintf (stderr, "felt %s\n", VERSION);
exit (0);
} else if (streq (arg, "-debug")) {
debug = 1;
} else if (streq (arg, "-preview")) {
preview = 1;
} else if (streq (arg, "-matrices")) {
matrices = 1;
} else if (streq (arg, "-summary")) {
summary = 1;
} else if (streq (arg, "-plot")) {
doplot = 1;
} else if (streq (arg, "+plot")) {
doplot = 0;
} else if (streq (arg, "-table")) {
dotable = 1;
} else if (streq (arg, "+table")) {
dotable = 0;
} else if (streq (arg, "-transfer")) {
dospectra = 0;
} else if (streq (arg, "-eigen")) {
domodal = 0;
} else if (streq (arg, "-orthonormal")) {
orthonormal = 1;
} else if (streq (arg, "-renumber")) {
renumber = 1;
} else if (streq (arg, "-details")) {
details = 1;
} else if (streq (arg, "-matlab")) {
if (++ i == *argc) {
fputs (usage, stderr);
return 1;
}
matlab = argv [i];
} else if (streq (arg, "-graphics")) {
if (++ i == *argc) {
fputs (usage, stderr);
return 1;
}
graphics = argv [i];
} else
argv [j ++] = arg;
argv [*argc = j] = NULL;
return 0;
}
/************************************************************************
* Function: main *
* *
* Description: Main is the driver function for the felt package. *
************************************************************************/
int main (argc, argv)
int argc;
char *argv [ ];
{
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 vectors */
Matrix x; /* eigenvectors */
Vector lambda; /* eigenvalues */
NodeDOF *forced; /* array of forced nodal DOF */
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 */
Matrix ttable; /* time step table */
unsigned *old_numbers; /* original node numbering */
AnalysisType mode; /* current analysis type */
/*
* Do everything to setup the problem
*/
# ifndef DOS
if (ParseCppOptions (&argc, argv)) {
fputs (usage, stderr);
exit (1);
}
# endif
if (ParseFeltOptions (&argc, argv)) {
fputs (usage, stderr);
exit (1);
}
if (argc > 2) {
fputs (usage, stderr);
exit (1);
}
add_all_definitions ( );
if (ReadFeltFile (argc == 2 ? argv [1] : "-"))
exit (1);
title = problem.title;
/*
* If debugging write the problem as we understand it to stdout
*/
if (debug)
WriteFeltFile ("-");
if (problem.num_nodes == 0 || problem.num_elements == 0)
Fatal ("nothing to do");
if (preview)
DrawStructureASCII (stdout, 78, 22);
/*
* Write a graphics file if necessary
*/
if (graphics != NULL) {
status = WriteGraphicsFile (graphics, 0.0);
if (status)
Fatal ("could not open graphics file for output");
}
/*
* if the user wanted analysis details we need
* to set the detail stream
*/
if (details)
SetDetailStream (stdout);
/*
* 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 (transient or static)
*/
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 (stdout, M, C, K);
if (matlab)
MatlabGlobalMatrices (matlab, M, C, K);
if (status)
Fatal ("%d fatal errors in stiffness and mass definitions",status);
if (analysis.step > 0.0) {
dtable = IntegrateHyperbolicDE (K, M, C);
ttable = NullMatrix;
}
else
dtable = RosenbrockHyperbolicDE (K, M, C, &ttable);
if (old_numbers != NULL)
RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);
if (dtable == NullMatrix)
Fatal ("fatal error in integration (probably a singularity).");
if (dotable)
WriteTransientTable (dtable, ttable, stdout);
if (doplot)
PlotTransientTable (dtable, ttable, analysis.step, stdout);
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 (stdout, M, NullMatrix, K);
if (matlab)
MatlabGlobalMatrices (matlab, 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 (old_numbers != NULL)
RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);
if (dotable)
WriteTransientTable (dtable, NullMatrix, stdout);
if (doplot)
PlotTransientTable (dtable, NullMatrix, analysis.step, stdout);
break;
case Static:
K = ConstructStiffness(&status);
if (status)
Fatal ("%d Fatal errors in element stiffness definitions", status);
if (matrices)
PrintGlobalMatrices (stdout, NullMatrix, NullMatrix, K);
if (matlab)
MatlabGlobalMatrices (matlab, NullMatrix, NullMatrix, K);
F = ConstructForceVector ( );
ZeroConstrainedDOF (K, F, &Kcond, &Fcond);
d = SolveForDisplacements (Kcond, Fcond);
if (d == NullMatrix)
Fatal ("could not solve for global displacements");
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);
WriteStructuralResults (stdout, title, R, numreactions);
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 (stdout, NullMatrix, NullMatrix, K);
if (matlab)
MatlabGlobalMatrices (matlab, 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 (dotable)
if (mode == StaticLoadCases)
WriteLoadCaseTable (dtable, stdout);
else
WriteLoadRangeTable (dtable, stdout);
if (doplot)
if (mode == StaticLoadCases)
PlotLoadCaseTable (dtable, stdout);
else
PlotLoadRangeTable (dtable, stdout);
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 (dotable)
WriteLoadRangeTable (dtable, stdout);
if (doplot)
PlotLoadRangeTable (dtable, stdout);
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);
WriteStructuralResults (stdout, title, NULL, 0);
break;
case StaticThermal:
K = ConstructStiffness(&status);
if (status)
Fatal ("%d Fatal errors in element stiffness definitions", status);
if (matrices)
PrintGlobalMatrices (stdout, NullMatrix, NullMatrix, K);
if (matlab)
MatlabGlobalMatrices (matlab, NullMatrix, NullMatrix, K);
F = ConstructForceVector ( );
ZeroConstrainedDOF (K, F, &Kcond, &Fcond);
d = SolveForDisplacements (Kcond, Fcond);
if (d == NullVector)
exit (1);
if (old_numbers != NULL)
RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);
WriteTemperatureResults (stdout, title);
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 (stdout, Mcond, Ccond, Kcond);
if (matlab)
MatlabGlobalMatrices (matlab, 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, stdout);
if (doplot)
PlotModeShapes (x, stdout);
if (domodal) {
FormModalMatrices (x, Mcond, Ccond, Kcond, &Mm, &Cm, &Km, orthonormal);
WriteModalResults (stdout, 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 (stdout, Mcond, Ccond, Kcond);
if (matlab)
MatlabGlobalMatrices (matlab, Mcond, Ccond, Kcond);
FindForcedDOF (&forced, &numforced);
H = ComputeTransferFunctions (Mcond, Ccond, Kcond, forced, numforced);
if (dospectra) {
S = ComputeOutputSpectra (H, forced, numforced);
if (S == NullMatrix)
break;
}
if (old_numbers != NULL)
RestoreNodeNumbers (problem.nodes, old_numbers, problem.num_nodes);
if (!dospectra) {
if (dotable)
WriteTransferFunctions (H, forced, numforced, stdout);
if (doplot)
PlotTransferFunctions (H, forced, numforced, stdout);
break;
}
if (dotable)
WriteOutputSpectra (S, stdout);
if (doplot)
PlotOutputSpectra (S, stdout);
break;
}
if (summary)
WriteMaterialStatistics (stdout);
exit (0);
}
syntax highlighted by Code2HTML, v. 0.9.1