/*
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;
}
syntax highlighted by Code2HTML, v. 0.9.1