OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removal of obsolete prototype code

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@815 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-02-22 17:19:30 +00:00 committed by Knut Morten Okstad
parent a91f3b90c1
commit 7be12bd1c1
12 changed files with 0 additions and 6512 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1510
src/Doxyfile.proto
View File

File diff suppressed because it is too large Load Diff

26
src/LinEqSystem.C
View File

@ -1,26 +0,0 @@
// $Id: LinEqSystem.C,v 1.1 2009-06-29 09:55:17 kmo Exp $
//==============================================================================
//!
//! \file LinEqSystem.C
//!
//! \date Apr 17 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Storage of a standard linear system of equations for structural FEM.
//!
//==============================================================================
#include "LinEqSystem.h"
void LinEqSystem::clear()
{
RHS.clear();
if (K) delete K;
if (M) delete M;
K = 0;
M = 0;
}

41
src/LinEqSystem.h
View File

@ -1,41 +0,0 @@
// $Id: LinEqSystem.h,v 1.4 2010-05-06 17:31:18 kmo Exp $
//==============================================================================
//!
//! \file LinEqSystem.h
//!
//! \date Apr 17 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Storage of a standard linear system of equations for structural FEM.
//!
//==============================================================================
#ifndef _LIN_EQ_SYSTEM_H
#define _LIN_EQ_SYSTEM_H
#include "SystemMatrix.h"
/*!
\brief Class for storage of a linear system of equations.
*/
class LinEqSystem
{
public:
SystemMatrix* K; //!< System stiffness matrix
SystemMatrix* M; //!< System mass matrix
StdVector RHS; //!< System right-hand-side vector
//! \brief Default constructor.
LinEqSystem() { K = 0; M = 0; }
//! \brief The destructor frees the dynamically allocated objects.
~LinEqSystem() { if (K) delete K; if (M) delete M; }
//! \brief Erases the system matrices and frees dynamically allocated storage.
void clear();
};
#endif

1376
src/LinearEl.C
View File

File diff suppressed because it is too large Load Diff

244
src/LinearEl.h
View File

@ -1,244 +0,0 @@
// $Id: LinearEl.h,v 1.19 2009-10-30 12:12:40 kmo Exp $
//==============================================================================
//!
//! \file LinearEl.h
//!
//! \date Jan 27 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Solver for linear elasticity problems using NURBS-based FEM.
//!
//==============================================================================
#ifndef _LINEAR_EL_H
#define _LINEAR_EL_H
#include "Vec3.h"
#include "Function.h"
#include "PressureLoad.h"
#include "LinEqSystem.h"
#include <map>
class VolumePatch;
class LocalSystem;
class SAMSpline;
class VTF;
/*!
\brief Struct for storage of data associated with one mode shape.
*/
struct Mode
{
double eigVal; //!< Eigenvalue associated with this mode
Vector eigVec; //!< Eigenvector associated with this mode
// \brief Default constructor setting \a eigVal to zero.
Mode() { eigVal = 0.0; }
};
/*!
\brief Result container for passing results to GLview export module.
*/
struct Result
{
Vector load; //!< The external load vector
Vector displ; //!< The displacement vector due to the external load
Matrix norms; //!< Element norms
std::vector<Mode> modes; //!< Mode shapes
bool freq; //!< Indicates whether eigenvalues are frequencies or not
// \brief Default constructor setting \a freq to \e false.
Result() { freq = false; }
};
/*!
\brief Driver class for the NURBS-based linear elasticity solver.
\details The class incapsulates data and methods for solving linear elasticity
problems using NURBS-based finite elements.
*/
class LinearEl
{
public:
//! \brief The constructor generates the model from the given input file.
//! \param[in] fileName Name of input file with model description
//! \param[in] checkRHS If \e true, ensure the model is in a right-hand system
//! \param[in] free If \e true, any specified boundary conditions are ignored
LinearEl(const char* fileName = 0, bool checkRHS = false, bool free = false);
//! \brief The destructor frees the dynamically allocated objects.
~LinearEl();
//! \brief Reads model data from the specified input file \a *fileName.
bool read(const char* fileName, bool checkRHS = false, bool free = false);
//! \brief Performs some preprocessing tasks on the finite element model.
//! \details This method should be invoked after \a readModel and before
//! any of the assembleK.... methods.
//! The main purpose of this method is to fill the SAMSpline object with data.
bool preprocess(const std::vector<int>& ignoredPatches = std::vector<int>(),
bool fixDup = false);
//! \brief Administers assembly of system stiffness matrix and load vector.
//! \param[in] solver Which equation solver (matrix format) to use
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
//! \param[out] R The external load vector in DOF-order (for visualization)
bool assembleKandR(SystemMatrix::Type solver, int nGauss, Vector* R);
//! \brief Administers assembly of the system stiffness- and mass matrices.
//! \param[in] solver Which equation solver (matrix format) to use
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
bool assembleKandM(SystemMatrix::Type solver, int nGauss);
//! \brief Administers assembly of the system stiffness matrices.
//! \param[in] solver Which equation solver (matrix format) to use
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
//! \param[in] dis Solution vector in DOF-order (for geometric stiffness)
bool assembleKandKg(SystemMatrix::Type solver, int nGauss, const Vector& dis);
//! \brief Administers assembly of the system stiffness matrix.
//! \param[in] solver Which equation solver (matrix format) to use
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
bool assembleKonly(SystemMatrix::Type solver, int nGauss);
//! \brief Solves the assembled linear system of equations for a given load.
//! \param[out] solution Solution vector, displacements at nodal points
bool solve(Vector& solution);
//! \brief Integrates and prints out some solution norm quantities.
//! \details If an analytical solution is provided, norms of the exact
//! error in the solution is computed as well.
//! \param[out] eNorm Element-wise norm quantities
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
//! \param[in] dis Solution vector in DOF-order, displacements at nodal points
bool solutionNorms(Matrix& eNorm, int nGauss, const Vector& dis);
//! \brief Performs a generalized eigenvalue analysis of the assembled system.
//! \param[in] iop Which eigensolver method to use
//! \param[in] nev Number of eigenvalues/vector (see ARPack documentation)
//! \param[in] ncv Number of Arnoldi vectors (see ARPack documentation)
//! \param[in] shift Eigenvalue shift
//! \param[out] solution Computed eigenvalues and associated eigenvectors
bool modes(int iop, int nev, int ncv, double shift,
std::vector<Mode>& solution);
//! \brief Writes the global grid and boundary conditions to a file.
//! \param[in] inputFile File name used to construct the grid file name from
//! \param[in] n Number of FE nodes over each knot-span
//! \param[in] nenod Number of nodes in each element to generate
bool writeGlobalGrid(const char* inputFile, const int* n,
int nenod = 8) const;
//! \brief Writes a VTF-file with the geometry and solution fields.
//! \param[in] inputFile File name used to construct the VTF-file name from
//! \param[in] solution The solution vector(s) to output
//! \param[in] nViz Number of visualization points over a knot-span
//! \param[in] format Format of VTF-file (0=ASCII, 1=BINARY)
//! \param[in] debug If \e true, output some additional debug quantities
//!
//! \details The NURBS patches are tesselated into linear hexahedrons with
//! a fixed number of HEX8-elements within each knot-span of non-zero length.
//! The solution fields (displacements and stresses) are then evaluated at the
//! nodal points of the generated HEX8-mesh and written to the VTF-file as
//! vector and scalar result fields.
bool writeGlv(const char* inputFile, const Result& solution,
const int* nViz, int format = 0, bool debug = false) const;
//! \brief Writes the geometry to g2-file for external viewing.
//! \param[in] g2file Filename for GoTools output file
void dumpGeometry(const char* g2file) const;
//! \brief Writes the solution vector to files for external viewing.
//! \param[in] solfile Filename prefix
//! \param[in] solution The global solution vector in DOF-order
//!
//! \details The solution vector is written to the ASCII file \a solfile.vec
//! in the same order as the control points are ordered in the NURBS patches.
//! In addition, each component is written to files \a solfile.u, \a solfile.v
//! and \a solfile.w respectively.
void dumpSolution(const char* solfile, const Vector& solution) const;
//! \brief Writes the system matrices to files for debugging.
//! \param[in] Kfile Filename for stiffness matrix
//! \param[in] Mfile Filename for mass matrix
//! \param[in] Rfile Filename for load vector
void dumpMat(const char* Kfile,
const char* Mfile = 0,
const char* Rfile = 0) const;
protected:
//! \brief Creates the computational FEM model.
bool createFEMmodel();
//! \brief Writes the grid geometry to the VTF-file.
//! \param vtf The VTF-file object to receive the geometry data
//! \param[in] nViz Number of visualization points over each knot-span
bool writeGlvG(VTF& vtf, const int* nViz) const;
//! \brief Writes the surface tractions for a given time step to VTF-file.
//! \param vtf The VTF-file object to receive the tractions
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter
bool writeGlvT(VTF& vtf, int iStep, int& nBlock) const;
//! \brief Writes boundary condition codes as scalar fields to the VTF-file.
//! \param vtf The VTF-file object to receive the boundary condition data
//! \param[in] nViz Number of visualization points over each knot-span
//! \param nBlock Running result block counter
//! \param[in] debug If \e true, output some additional debug quantities
bool writeGlvBC(VTF& vtf, const int* nViz, int& nBlock, bool debug) const;
//! \brief Writes the load vector for a given time step to VTF-file.
//! \param vtf The VTF-file object to receive the load vector
//! \param[in] load The load vector to output
//! \param[in] nViz Number of visualization points over each knot-span
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter
bool writeGlvR(VTF& vtf, const Vector& load,
const int* nViz, int iStep, int& nBlock) const;
//! \brief Writes solution fields for a given load/time step to the VTF-file.
//! \details If an analytical solution is provided, the exact stress fields
//! are written to the VTF-file as well.
//! \param vtf The VTF-file object to receive the solution fields
//! \param[in] solution The solution vector to derive the result fields from
//! \param[in] nViz Number of visualization points over each knot-span
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter
bool writeGlvS(VTF& vtf, const Vector& solution,
const int* nViz, int iStep, int& nBlock) const;
//! \brief Writes an eigenvector and associated eigenvalue to the VTF-file.
//! \details The eigenvalue is used as a label on the step state info
//! that is associated with the eigenvector.
//! \param vtf The VTF-file object to receive the eigenvector
//! \param[in] mode The eigenvector to output
//! \param[in] freq \e true if the eigenvalue is a frequency
//! \param[in] nViz Number of visualization points over each knot-span
//! \param[in] iMode Mode shape identifier
//! \param nBlock Running result block counter
bool writeGlvM(VTF& vtf, const Mode& mode, bool freq,
const int* nViz, int iMode, int& nBlock) const;
//! \brief Writes element norms for a given load/time step to the VTF-file.
//! \details This method can be used only when the number of visualization
//! points over each knot-span equals 2 (that is, no additonal points).
//! \param vtf The VTF-file object to receive the element norms
//! \param[in] norms The element norms to output
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter
bool writeGlvN(VTF& vtf, const Matrix& norms, int iStep, int& nBlock) const;
private:
std::vector<VolumePatch*> model; //!< The actual NURBS/spline model
std::vector<PLoad> load; //!< Surface pressure loads
Vec3 gravity; //!< Gravitation vector
LocalSystem* rCS; //!< Local coordinate system for results
std::map<Vec3,Vec3> trac; //!< Evaluated surface tractions
TensorFunc* asol; //!< Analytical stress field
SAMSpline* sam; //!< Data for FE assembly management
LinEqSystem sys; //!< The linear equation system
};
#endif

37
src/PressureLoad.h
View File

@ -1,37 +0,0 @@
// $Id: PressureLoad.h,v 1.2 2009-05-05 09:25:09 kmo Exp $
//==============================================================================
//!
//! \file PressureLoad.h
//!
//! \date Jan 27 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Representation of a constant pressure load.
//!
//==============================================================================
#ifndef _PRESSURELOAD_H
#define _PRESSURELOAD_H
/*!
\brief Struct for representing a constant pressure load.
*/
struct PLoad
{
int volp; //!< Volume spline patch index [0,nvol]
int face; //!< Local face index on spline volume [-3,3]
int pdir; //!< Pressure direction [0,3] (0=normal pressure)
double p; //!< Actual pressure value
//! \brief Constructor creating an initialized pressure instance.
//! \param[in] s Patch index
//! \param[in] f Local face index
//! \param[in] d Pressure direction
//! \param[in] v Pressure value
PLoad(int s, int f, int d, double v) : volp(s), face(f), pdir(d), p(v) {}
};
#endif

37
src/README
View File

@ -1,37 +0,0 @@
$Id: README,v 1.3 2009-12-03 12:20:44 kmo Exp $
This is the source of the prototype spline/NURBS-based FEM solver
developed for the ICADA project at SINTEF ICT in Trondheim.
The Makefile is for compilation on Linux.
It assumes that the libraries liblapack and libblas are installed.
In addition, the following dependencies on 3rd party libraries are present:
Required:
GoTools - Download and install from the source code repository at
http://svn.math.sintef.no. See also the wiki documentation at
http://svn.math.sintef.no/index.php/Some_info_about_the_GoTools_and_SISL_SVN_repositories
The libraries are assumed found in the sub-folder GoTools/lib
The header files are assumed found in the sub-folder GoTools/include
The GoTools packages require that the libboost package is installed.
ARPACK - Download and compile from http://www.caam.rice.edu/software/ARPACK
Optional:
SAM - A few fortran subroutine from the SAM library by Kolbein Bell
can be used for the assembly of element matrices into system matrices.
In addition, the sparse solver SPR from this library can be used.
It is possible to compile without the SAM library, by commenting out
the symbols SAMOPT and SAMLIB in the Makefile.
SuperLU - This is a public domain package for direct solution of non-symmetric
sparse equation systems. The source code is available from
http://crd.lbl.gov/~xiaoye/SuperLU. To compile without it, comment out
the symbol definitions SLUOPT and SLULIB in the Makefile
SAMG - This is a commercial package (algebraic multigrid solver) and is
included for testing only. To compile without it, comment out the
symbol definitions SAMGOPT and SAMGLIB in the Makefile.
GLviewExpressWriter - Ceetron's VTF API. If you don't have it, comment out the
symbol defintions VTFOPT and VTFLIB in the Makefile.

211
src/SAMSpline.C
View File

@ -1,211 +0,0 @@
// $Id: SAMSpline.C,v 1.3 2009-08-18 08:49:56 kmo Exp $
//==============================================================================
//!
//! \file SAMSpline.C
//!
//! \date Dec 10 2008
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Assembly of FE matrices into system matrices for SplineVolume models.
//!
//==============================================================================
#include "SAMSpline.h"
#include "VolumePatch.h"
#include "LinEqSystem.h"
#include "MPC.h"
bool SAMSpline::init (const VolumePatch& patch)
{
// Initialize some model size parameters
nnod = patch.getNoNodes();
nel = patch.getNoElms();
nceq = patch.getNoMPCs();
ndof = 3*nnod;
std::vector<VolumePatch*> smod(1,(VolumePatch*)&patch);
// Initialize the node/dof arrays (madof,msc)
initNodeDofs(smod);
// Initialize the element connectivity arrays (mpmnpc,mmnpc)
initElementConn(smod);
// Initialize the constraint equation arrays (mpmceq,mmceq,ttcc)
initConstraintEqs(smod);
// Initialize the dof-to-equation connectivity array (meqn)
return initSystemEquations();
}
bool SAMSpline::init (const std::vector<VolumePatch*>& smod, int numNod)
{
// Initialize some model size parameters
nnod = numNod;
for (size_t i = 0; i < smod.size(); i++)
{
if (numNod == 0) nnod += smod[i]->getNoNodes();
nel += smod[i]->getNoElms();
nceq += smod[i]->getNoMPCs();
}
ndof = 3*nnod;
std::cout <<"\n\n >>> SAM model summary <<<"
<<"\nNumber of elements "<< nel
<<"\nNumber of nodes "<< nnod
<<"\nNumber of dofs "<< ndof << std::endl;
// Initialize the node/dof arrays (madof,msc)
initNodeDofs(smod);
// Initialize the element connectivity arrays (mpmnpc,mmnpc)
initElementConn(smod);
// Initialize the constraint equation arrays (mpmceq,mmceq,ttcc)
initConstraintEqs(smod);
// Initialize the dof-to-equation connectivity array (meqn)
bool status = initSystemEquations();
std::cout <<"Number of unknowns "<< neq << std::endl;
return status;
}
void SAMSpline::initNodeDofs (const std::vector<VolumePatch*>& smod)
{
if (nnod < 1) return;
// Initialize the node and dof arrays
madof = new int[nnod+1];
msc = new int[ndof];
int n; madof[0] = 1;
for (n = 0; n < nnod; n++)
madof[n+1] = madof[n] + 3;
for (n = 0; n < ndof; n++)
msc[n] = 1;
std::vector<VolumePatch::BC>::const_iterator bit;
for (size_t j = 0; j < smod.size(); j++)
for (bit = smod[j]->begin_BC(); bit != smod[j]->end_BC(); bit++)
{
n = bit->node;
msc[3*n-3] *= bit->CX;
msc[3*n-2] *= bit->CY;
msc[3*n-1] *= bit->CZ;
}
}
void SAMSpline::initElementConn (const std::vector<VolumePatch*>& smod)
{
if (nel < 1) return;
// Find the size of the element connectivity array
size_t j;
IntMat::const_iterator eit;
for (j = 0; j < smod.size(); j++)
for (eit = smod[j]->begin_elm(); eit != smod[j]->end_elm(); eit++)
nmmnpc += eit->size();
// Initialize the element connectivity arrays
mpmnpc = new int[nel+1];
mmnpc = new int[nmmnpc];
int ip = mpmnpc[0] = 1;
for (j = 0; j < smod.size(); j++)
for (eit = smod[j]->begin_elm(); eit != smod[j]->end_elm(); eit++, ip++)
{
mpmnpc[ip] = mpmnpc[ip-1];
for (size_t i = 0; i < eit->size(); i++)
mmnpc[(mpmnpc[ip]++)-1] = smod[j]->getNodeID(1+(*eit)[i]);
}
}
void SAMSpline::initConstraintEqs (const std::vector<VolumePatch*>& smod)
{
// Estimate the size of the constraint equation array
size_t j;
MPCIter cit;
for (j = 0; j < smod.size(); j++)
for (cit = smod[j]->begin_MPC(); cit != smod[j]->end_MPC(); cit++)
nmmceq += 1 + (*cit)->getNoMaster();
// Initialize the constraint equation arrays
mpmceq = new int[nceq+1];
int ip = mpmceq[0] = 1;
if (nceq < 1) return;
mmceq = new int[nmmceq];
ttcc = new real[nmmceq];
for (j = 0; j < smod.size(); j++)
for (cit = smod[j]->begin_MPC(); cit != smod[j]->end_MPC(); cit++, ip++)
{
mpmceq[ip] = mpmceq[ip-1];
// Slave dof ...
int idof = madof[(*cit)->getSlave().node-1] + (*cit)->getSlave().dof - 1;
if (msc[idof-1] == 0)
{
std::cerr <<"SAM: Ignoring constraint equation for dof "
<< idof <<" ("<< (*cit)->getSlave()
<<").\n This dof is already marked as FIXED."<< std::endl;
ip--;
nceq--;
continue;
}
else if (msc[idof-1] < 0)
{
std::cerr <<"SAM: Ignoring constraint equation for dof "
<< idof <<" ("<< (*cit)->getSlave()
<<").\n This dof is already marked as SLAVE."<< std::endl;
ip--;
nceq--;
continue;
}
int ipslv = (mpmceq[ip]++) - 1;
mmceq[ipslv] = idof;
ttcc[ipslv] = (*cit)->getSlave().coeff;
msc[idof-1] = -ip;
// Master dofs ...
for (size_t i = 0; i < (*cit)->getNoMaster(); i++)
{
idof = madof[(*cit)->getMaster(i).node-1] + (*cit)->getMaster(i).dof-1;
if (msc[idof-1] > 0)
{
int ipmst = (mpmceq[ip]++) - 1;
mmceq[ipmst] = idof;
ttcc[ipmst] = (*cit)->getMaster(i).coeff;
}
else if (msc[idof-1] < 0)
{
// Master dof is constrained (unresolved chaining)
std::cerr <<"SAM: Chained MPCs detected"
<<", slave "<< (*cit)->getSlave()
<<", master "<< (*cit)->getMaster(i)
<<" (ignored)."<< std::endl;
mpmceq[ip] = mpmceq[ip-1];
ip--;
nceq--;
break;
}
}
}
// Reset the negative values in msc before calling SYSEQ
for (ip = 0; ip < ndof; ip++)
if (msc[ip] < 0) msc[ip] = 0;
}
bool SAMSpline::initForAssembly (LinEqSystem& sys, bool withRHS) const
{
if (sys.K) sys.K->initAssembly(*this);
if (sys.M) sys.M->initAssembly(*this);
if (withRHS) sys.RHS.resize(neq,true);
return neq > 0 ? true : false;
}

63
src/SAMSpline.h
View File

@ -1,63 +0,0 @@
// $Id: SAMSpline.h,v 1.2 2009-05-27 12:52:34 kmo Exp $
//==============================================================================
//!
//! \file SAMSpline.h
//!
//! \date Dec 10 2008
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Assembly of FE matrices into system matrices for SplineVolume models.
//!
//==============================================================================
#ifndef _SAM_SPLINE_H
#define _SAM_SPLINE_H
#include "SAM.h"
class VolumePatch;
class LinEqSystem;
/*!
\brief This is a sub-class of SAM with added functionality for spline models.
\details It contains some additional functions for initializing the SAM-data
for an FE model based on SplineVolume patches.
*/
class SAMSpline : public SAM
{
public:
//! \brief Default constructor which initializes an empty SAMSpline object.
SAMSpline() : SAM() {}
//! \brief Empty destructor.
virtual ~SAMSpline() {}
//! \brief Allocates the dynamic arrays and populates them with data
//! for a single-patch spline model.
bool init(const VolumePatch& model);
//! \brief Allocates the dynamic arrays and populates them with data
//! for a multi-patch spline model.
bool init(const std::vector<VolumePatch*>& model, int numNod = 0);
//! \brief Performs the necessary initialization of the system matrices
//! prior to the element assembly.
//! \details This method must be called once before the first call to
//! \a assembleSystem for a given load case or time step.
//! \param sys The system left-hand-side matrices and right-hand-side vector
//! \param withRHS If \a true, initialize the right-hand-side vector too
//! \return \e false if no free DOFs in the system, otherwise \e true
bool initForAssembly(LinEqSystem& sys, bool withRHS = false) const;
protected:
//! \brief Initializes the nodal arrays \a MINEX, \a MADOF and \a MSC.
void initNodeDofs(const std::vector<VolumePatch*>& p);
//! \brief Initializes the element topology arrays \a MPMNPC and \a MMNPC.
void initElementConn(const std::vector<VolumePatch*>& p);
//! \brief Initializes the multi-point constraint arrays
//! \a MPMCEQ, \a MMCEQ and \a TTCC.
void initConstraintEqs(const std::vector<VolumePatch*>& p);
};
#endif

2241
src/VolumePatch.C
View File

File diff suppressed because it is too large Load Diff

448
src/VolumePatch.h
View File

@ -1,448 +0,0 @@
// $Id: VolumePatch.h,v 1.22 2010-06-22 11:38:05 kmo Exp $
//==============================================================================
//!
//! \file VolumePatch.h
//!
//! \date Dec 10 2008
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Representation of a topological cube as a SplineVolume.
//! \details The class contains the necessary data and methods to establish
//! a finite element algebraic system of equations for the stiffness relation
//! of the volume, using NURBS as basis functions via the GoTools library.
//!
//==============================================================================
#ifndef _VOLUME_PATCH_H
#define _VOLUME_PATCH_H
#include "Function.h"
#include "MPCLess.h"
#include "MatVec.h"
#include <set>
#include <map>
namespace Go {
class SplineVolume;
}
struct ElementBlock;
class SystemVector;
class LinEqSystem;
class LocalSystem;
class SAM;
class Vec3;
class Tensor;
typedef std::vector<int> IntVec; //!< General integer vector
typedef std::vector<IntVec> IntMat; //!< General 2D integer matrix
typedef std::set<MPC*,MPCLess> MPCSet; //!< Sorted set of MPC equations
typedef MPCSet::const_iterator MPCIter; //!< Iterator over a MPC equation set
/*!
\brief Assembly of finite element contributions from a SplineVolume object.
\details This class stores data and provides methods for the calculation of
element stiffness and mass matrices, as well as load vectors for a structured
isogeometric grid, using NURBS or splines as basis functions.
*/
class VolumePatch
{
//! \brief Struct for nodal point data.
struct IJK
{
int I; //!< Index in first parameter direction
int J; //!< Index in second parameter direction
int K; //!< Index in third parameter direction
int global; //!< Global node number [1,nnod]
};
public:
//! \brief Struct for boundary condition codes.
struct BC
{
int node; //!< Global node number [1,nnod]
char CX; //!< Boundary condition code for X-translation
char CY; //!< Boundary condition code for Y-translation
char CZ; //!< Boundary condition code for Z-translation
//! \brief Constructor initializing a BC instance.
BC(int n, char x, char y, char z) : node(n), CX(x), CY(y), CZ(z) {}
};
//! \brief Constructor creating an instance by reading the given file.
VolumePatch(const char* fileName, bool checkRHS = false);
//! \brief Constructor creating an instance by reading the given input stream.
VolumePatch(std::istream& is, bool checkRHS = false);
//! \brief Default constructor creating an empty patch.
VolumePatch() { svol = 0; swapW = false; E = nu = rho = 0.0; }
//! \brief The Destructor frees the dynamically allocated SplineVolume object.
~VolumePatch();
//! \brief Creates an instance by reading the given input stream, \a is.
bool read(std::istream& is);
//! \brief Writes the geometry of the SplineVolume object to the given stream.
bool write(std::ostream& os) const;
//! \brief Check that the patch is modelled in a right-hand-side system.
//! \details If it isn't, the w-parameter direction is swapped.
bool checkRightHandSystem();
//! \brief Refine the parametrization by inserting extra knots uniformly.
//! \param[in] dir Parameter direction to refine
//! \param[in] nInsert Number of extra knots to insert in each knot-span
bool uniformRefine(int dir, int nInsert);
//! \brief Raise the order of the SplineVolume object for this patch.
//! \param[in] ru Number of times to raise the order i u-direction
//! \param[in] rv Number of times to raise the order i u-direction
//! \param[in] rw Number of times to raise the order i u-direction
bool raiseOrder(int ru, int rv, int rw);
//! \brief Generates the finite element topology data for the patch.
bool generateFEMTopology();
//! \brief Clears the contents of the patch, making it empty.
void clear();
//! \brief Checks if the patch is empty.
bool empty() const { return svol == 0; }
//! \brief Returns local 1-based index of the node with given global number.
//! \details If the given node number is not present, 0 is returned.
//! \param[in] globalNum Global node number
int getNodeIndex(int globalNum) const;
//! \brief Returns the global node number for the given node.
//! \param[in] inod 1-based node index local to current patch
int getNodeID(int inod) const;
//! \brief Returns the global coordinates for the given node.
//! \param[in] inod 1-based node index local to current patch
Vec3 getCoord(int inod) const;
//! \brief Returns the total number of nodes in the patch.
size_t getNoNodes() const { return nodeInd.size(); }
//! \brief Returns the total number of elements in the patch.
size_t getNoElms() const { return MNPC.size(); }
//! \brief Returns the total number of MPC equations in the patch.
size_t getNoMPCs() const { return mpcs.size(); }
//! \brief Returns the beginning of the BC array.
std::vector<BC>::const_iterator begin_BC() const { return BCode.begin(); }
//! \brief Returns the end of the BC array.
std::vector<BC>::const_iterator end_BC() const { return BCode.end(); }
//! \brief Returns the beginning of the MNPC array.
IntMat::const_iterator begin_elm() const { return MNPC.begin(); }
//! \brief Returns the end of the MNPC array.
IntMat::const_iterator end_elm() const { return MNPC.end(); }
//! \brief Returns the beginning of the MPC set.
MPCIter begin_MPC() const { return mpcs.begin(); }
//! \brief Returns the end of the MPC set.
MPCIter end_MPC() const { return mpcs.end(); }
//! \brief Merges a given node in this patch with a given global node.
//! \param[in] node 1-based node index local to current patch
//! \param[in] globalNum Global number of the node to merge \a node with
bool mergeNodes(int node, int globalNum);
//! \brief Renumbers the global node numbers referred by this patch.
//! \details The node numbers referred by boundary condition and
//! multi-point constraint objects in the patch are renumbered.
//! The noded themselves are assumed to already be up to date.
//! \param[in] old2new Old-to-new node number mapping
//! \param[in] silent If \e false, give error on missing node in \a old2new
bool renumberNodes(const std::map<int,int>& old2new, bool silent);
//! \brief Renumbers all global node numbers in the entire model.
//! \param[in] model All volume patches in the model
//! \return The number of unique nodes in the model
//!
//! \details After the renumbering, the global node numbers are in the range
//! [1,\a nnod ], where \a nnod is the number of unique nodes in the model.
static int renumberNodes(const std::vector<VolumePatch*>& model);
//! \brief Resolves (possibly multi-level) chaining in MPC equations.
//! \param[in] model All volume patches in the model
//!
//! \details If a master DOF in one MPC is specified as slave by another MPC,
//! it is replaced by the master(s) of that other equation.
//! Since an MPC-equation may couple nodes belonging to different patches,
//! this method must have all patches in the model available.
static void resolveMPCchains(const std::vector<VolumePatch*>& model);
//! \brief Defines material properties for current volume patch.
//! \param[in] Emod Young's modulus
//! \param[in] Poiss Poisson's ratio
//! \param[in] Density Mass density
void setMaterial(double Emod, double Poiss, double Density)
{ E = Emod; nu = Poiss; rho = Density; }
//! \brief Retrieves material properties for current volume patch.
//! \param[out] Emod Young's modulus
//! \param[out] Poiss Poisson's ratio
//! \param[out] Density Mass density
void getMaterial(double& Emod, double& Poiss, double& Density)
{ Emod = E; Poiss = nu; Density = rho; }
//! \brief Makes two opposite boundary faces periodic.
//! \param[in] dir Parameter direction defining the periodic faces
void closeFaces(int dir);
//! \brief Constrains all DOFs on a given boundary face.
//! \param[in] dir Parameter direction defining the face to constrain
//! \param[in] dof Which DOFs to constrain at each node on the face
//! \param[in] value Prescribed value of the constrained DOFs
void constrainFace(int dir, int dof = 123, double value = 0.0);
//! \brief Constrains all DOFs along a line on a given boundary face.
//! \param[in] fdir Parameter direction defining the face to constrain
//! \param[in] ldir Parameter direction defining the line to constrain
//! \param[in] xi Parameter value defining the line to constrain
//! \param[in] dof Which DOFs to constrain at each node along the line
//! \param[in] value Prescribed value of the constrained DOFs
//!
//! \details The parameter \a xi has to be in the domain [0.0,1.0], where
//! 0.0 means the beginning of the domain and 1.0 means the end. The line to
//! constrain goes along the parameter direction \a ldir in the face with
//! with normal in parameter direction \a fdir, and positioned along the third
//! parameter direction as indicated by \a xi. The actual value of \a xi
//! is converted to the integer value closest to \a xi*n, where \a n is the
//! number of nodes (control points) in that parameter direction.
void constrainLine(int fdir, int ldir, double xi,
int dof = 123, double value = 0.0);
//! \brief Constrains a corner node identified by the three parameter indices.
//! \param[in] I Parameter index in u-direction
//! \param[in] J Parameter index in v-direction
//! \param[in] K Parameter index in w-direction
//! \param[in] dof Which DOFs to constrain at the node
//! \param[in] value Prescribed value of the constrained DOFs
//!
//! \details The sign of the three indices is used to define whether we want
//! the node at the beginning or the end of that parameter direction.
//! The magnitude of the indices are not used.
void constrainCorner(int I, int J, int K, int dof = 123, double value = 0.0);
//! \brief Constrains a node identified by three relative parameter values.
//! \param[in] xi Parameter in u-direction
//! \param[in] eta Parameter in v-direction
//! \param[in] zeta Parameter in w-direction
//! \param[in] dof Which DOFs to constrain at the node
//! \param[in] value Prescribed value of the constrained DOFs
//!
//! \details The parameter values have to be in the domain [0.0,1.0], where
//! 0.0 means the beginning of the domain and 1.0 means the end. For values
//! in between, the actual index is taken as the integer value closest to
//! \a r*n, where \a r denotes the given relative parameter value,
//! and \a n is the number of nodes along that parameter direction.
void constrainNode(double xi, double eta, double zeta,
int dof = 123, double value = 0.0);
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch
//! \param[in] nface Local face index of neighbor patch, in range [1,6]
//! \param[in] norient Relative face orientation flag (see below)
//!
//! \details The face orientation flag \a norient must be in range [0,7].
//! When interpreted as a binary number, its 3 digits are decoded as follows:
//! - left digit = 1: The u and v parameters of the neighbor face are swapped
//! - middle digit = 1: Parameter \a u in neighbor patch face is reversed
//! - right digit = 1: Parameter \a v in neighbor patch face is reversed
bool connectPatch(int face, VolumePatch& neighbor,
int nface, int norient = 0);
//! \brief Assembles coefficient matrices and right-hand-side vector.
//! \param sys The linear system of equations
//! \param[in] sam Data for finite element assembly management
//! \param[in] gravity Gravitation vector
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
//! \param[in] displ Solution vector in DOF-order
bool assembleSystem(LinEqSystem& sys, const SAM& sam, const Vec3& gravity,
int nGauss = 4, const Vector& displ = Vector());
//! \brief Assembles right-hand-side vector due to surface traction on a face.
//! \param S The right-hand-side vector
//! \param[in] sam Data for finite element assembly management
//! \param[in] t The surface traction function
//! \param[in] dir Local index of the face subjected to the traction
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
//! \param[out] trac Evaluated tractions at integration points (optional)
bool assembleForces(SystemVector& S, const SAM& sam, const TractionFunc& t,
int dir, int nGauss = 4, std::map<Vec3,Vec3>* trac = 0);
//! \brief Evaluates some norms of the finite element solution.
//! \details The energy norm of the solution is computed by numerical
//! integration of \f$a({\bf u}^h,{\bf u}^h)\f$ over the patch.
//! If an analytical solution is available, the norm of the exact error
//! \f$a({\bf u}-{\bf u}^h,{\bf u}-{\bf u}^h)\f$ is computed as well.
//! \param[out] gNorm Global norms
//! \param[out] eNorm Element norms
//! \param[in] nGauss Numerical integration scheme (number of points in 1D)
//! \param[in] displ Solution vector in DOF-order
//! \param[in] sol Pointer to analytical stress field (optional)
bool solutionNorms(Vector& gNorm, Matrix& eNorm, int nGauss,
const Vector& displ, const TensorFunc* sol = 0);
//! \brief Creates a hexahedron element model of this patch for visualization.
//! \param[out] grid The generated hexahedron grid
//! \param[in] npe Number of visualization nodes over each knot span
//! \note The number of element nodes must be set in \a grid on input.
bool convertToElementBlock(ElementBlock& grid, const int* npe) const;
//! \brief Extracts element results for this patch from a global vector.
//! \param[in] globRes Global matrix of element results
//! \param[out] elmRes Element results for this patch
void extractElmRes(const Matrix& globRes, Matrix& elmRes) const;
//! \brief Extracts nodal displacements for this patch from the global vector.
//! \param[in] solution Global solution vector in DOF-order
//! \param[out] displ Nodal displacement vector for this patch
void extractSolution(const Vector& solution, Vector& displ) const;
//! \brief Evaluates the displacement field at all visualization points.
//! \param[out] dField Displacement field
//! \param[in] displ Solution vector in DOF-order
//! \param[in] npe Number of visualization nodes over each knot span
//! \param[in] cs Local coordinate system
bool evalDisplField(Matrix& dField, const Vector& displ,
const int* npe, const LocalSystem* cs = 0) const;
//! \brief Evaluates the stress field at all visualization points.
//! \param[out] sField Stress field
//! \param[in] displ Solution vector in DOF-order
//! \param[in] npe Number of visualization nodes over each knot span
//! \param[in] cs Local coordinate system
bool evalStressField(Matrix& sField, const Vector& displ,
const int* npe, const LocalSystem* cs = 0) const;
//! \brief Calculates von Mises stress field from a stress tensor field.
//! \param[in] sigma Stress tensor field
//! \param[out] vm von Mises stress field
static void vonMises(const Matrix& sigma, Vector& vm);
//! \brief If \e true, matching nodes in two adjacent patches are merged.
static bool mergeDuplNodes; //!< If \e false, MPC-equations are used.
static bool swapJac; //!< Should the sign of the Jacobian be swapped?
static int splineEvalMethod; //!< Which spline evaluation method to use
protected:
//! \brief Adds a general multi-point-constraint equation to the patch.
//! \param mpc Pointer to an MPC-object
bool addMPC(MPC* mpc);
//! \brief Creates and adds a single-point constraint to the patch.
//! \param[in] node Global node number of the node to constrain
//! \param[in] dir Which local DOF to constrain (1, 2, 3)
//! \param[in] value The prescribed value
bool addSPC(int node, int dir, double value = 0.0);
//! \brief Creates and adds a periodicity constraint to the patch.
//! \param[in] master Global node number of the master node
//! \param[in] slave Global node number of the slave node to constrain
//! \param[in] dir Which local DOF to constrain (1, 2, 3)
bool addPeriodicity(int master, int slave, int dir);
//! \brief Creates periodicity constraints between two nodes in the patch.
//! \param[in] master Global node number of the master node
//! \param[in] slave Global node number of the slave node to constrain
//! \param[in] code Which local DOFs to constrain (1, 2, 3, 12, 23, 123)
void makePeriodic(int master, int slave, int code = 123);
//! \brief Constrains DOFs in the given node to the given value.
//! \param[in] node 1-based node index local to current patch
//! \param[in] code Which local DOFs to constrain (1, 2, 3, 12, 23, 123)
//! \param[in] value The prescribed value
void prescribe(int node, int code = 123, double value = 0.0);
//! \brief Constrains DOFs in the given node to zero.
//! \param[in] node 1-based node index local to current patch
//! \param[in] code Which local DOFs to constrain (1, 2, 3, 12, 23, 123)
void fix(int node, int code = 123);
//! \brief Recursive method used by \a resolveMPCchains.
static bool resolveMPCchain(const MPCSet& allMPCs, MPC* mpc);
//! \brief Merges the given two nodes into one.
//! \details The global node number of the node with the highest number
//! is changed into the number of the other node.
static void mergeNodes(IJK& node1, IJK& node2);
//! \brief Calculates the parameter values for all visualization points in
//! one direction.
//! \param[out] par Parameter values for all visualization points.
//! \param[in] dir Parameter direction (0,1,2)
//! \param[in] nSegSpan Number of visualization segments over each non-zero
//! knot-spans
bool getGridParameters(RealArray& par, int dir, int nSegSpan) const;
//! \brief Returns the volume in the parameter space for an element.
//! \param[in] iel Element index
double getParametricVolume(int iel) const;
//! \brief Returns boundary face area in the parameter space for an element.
//! \param[in] iel Element index
//! \param[in] dir Local face index of the boundary face
double getParametricArea(int iel, int dir) const;
//! \brief Sets up the constitutive matrix for this patch.
//! \param[out] C 6\f$\times\f$6-matrix, representing the material tensor
//! \param[in] inverse If \e true, set up the inverse matrix instead
bool getMaterialMatrix(Matrix& C, bool inverse = false) const;
//! \brief Returns a matrix with nodal coordinates for an element.
//! \param[in] iel Element index
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
bool getElementCoordinates(Matrix& X, int iel) const;
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in the patch
void getNodalCoordinates(Matrix& X) const;
//! \brief Returns an index into the internal coefficient array for a node.
//! \param[in] inod 1-based node index local to current patch
int coeffInd(int inod) const;
//! \brief Evaluates the displacement field at all visualization points.
//! \details Memory intensive version, used when \a splineEvalMethod = 1.
//! \param[out] dField Displacement field
//! \param[in] displ Solution vector in DOF-order
//! \param[in] gpar Parameter values of the vizualization points
//! \param[in] cs Local coordinate system
bool evalDisplField1(Matrix& dField, const Vector& displ,
const RealArray* gpar, const LocalSystem* cs = 0) const;
//! \brief Evaluates the displacement field at all visualization points.
//! \details More efficient version, used when \a splineEvalMethod = 2.
//! \param[out] dField Displacement field
//! \param[in] displ Solution vector in DOF-order
//! \param[in] gpar Parameter values of the vizualization points
//! \param[in] cs Local coordinate system
bool evalDisplField2(Matrix& dField, const Vector& displ,
const RealArray* gpar, const LocalSystem* cs = 0) const;
//! \brief Evaluates the stress field at all visualization points.
//! \details Memory intensive version, used when \a splineEvalMethod = 1.
//! \param[out] sField Stress field
//! \param[in] displ Solution vector in DOF-order
//! \param[in] gpar Parameter values of the vizualization points
//! \param[in] cs Local coordinate system
bool evalStressField1(Matrix& sField, const Vector& displ,
const RealArray* gpar, const LocalSystem* cs = 0) const;
//! \brief Evaluates the stress field at all visualization points.
//! \details More efficient version, used when \a splineEvalMethod = 2.
//! \param[out] sField Stress field
//! \param[in] displ Solution vector in DOF-order
//! \param[in] gpar Parameter values of the vizualization points
//! \param[in] cs Local coordinate system
bool evalStressField2(Matrix& sField, const Vector& displ,
const RealArray* gpar, const LocalSystem* cs = 0) const;
//! \brief Calculates the Strain-displacement matrix \b B.
static void formBmatrix(Matrix& B, const Matrix& dNdX);
private:
// Geometry and physical properties
Go::SplineVolume* svol; //!< Pointer to the actual spline volume object
bool swapW; //!< Has the w-parameter direction been swapped?
double E; //!< Young's modulus
double nu; //!< Poisson's ratio
double rho; //!< Mass density
// Finite element data structures
IntMat MNPC; //!< Matrix of Nodal Point Correspondance
std::vector<int> MLGE; //!< Matrix of Local to Global Element numbers
std::vector<IJK> nodeInd; //!< IJK-triplets for the control points (nodes)
std::vector<BC> BCode; //!< Boundary condition codes
MPCSet mpcs; //!< All multi-point constraints
static int gEl; //!< Global element counter
static int gNod; //!< Global node counter
};
#endif

278
src/main_LinEl.C
View File

@ -1,278 +0,0 @@
// $Id: main_LinEl.C,v 1.20 2010-05-06 17:31:18 kmo Exp $
//==============================================================================
//!
//! \file main_LinEl.C
//!
//! \date Jan 28 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Main program for the NURBS-based linear elasticity solver.
//!
//==============================================================================
#include "LinearEl.h"
#include "VolumePatch.h"
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
/*!
\brief Main program for the NURBS-based linear elasticity solver.
The input to the program is specified through the following
command-line arguments. The arguments may be given in arbitrary order.
\arg \a input-file : Input file with model definition
\arg -dense : Use the dense LAPACK matrix equation solver
\arg -superlu : Use the sparse superLU equation solver
\arg -samg : Use the sparse algebraic multi-grid equation solver
\arg -swapJ : Swap the sign of the Jacobian to account for inverse ordering
\arg -splineMethod \a iop : Spline evaluation method (1 or 2)
\arg -nGauss \a n : Number of Gauss points over a knot-span in each direction
\arg -vtf \a format : VTF-file format (0=ASCII, 1=BINARY)
\arg -nviz \a nviz : Number of visualization points over each knot-span
\arg -nu \a nu : Number of visualization points per knot-span in u-direction
\arg -nv \a nv : Number of visualization points per knot-span in v-direction
\arg -nw \a nw : Number of visualization points per knot-span in w-direction
\arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
\arg -eig \a iop : Eigenproblem solver to use (1...6)
\arg -nev \a nev : Number of eigenvalues to compute
\arg -ncv \a ncv : Number of Arnoldi vectors to use in the eigenvalue analysis
\arg -shift \a shf : Shift value to use in the eigenproblem solver
\arg -free : Ignore all boundary conditions (use in free vibration analysis)
\arg -check : Data check only, read model and output to VTF (no solution)
\arg -checkRHS : Check that the patches are modelled in a right-hand system
\arg -fixDup : Resolve co-located nodes by merging them into a single node
\arg -vizRHS : Save the right-hand-side load vector on the VTF-file
\arg -dumpGeo : Dumps the (possibly refined) geometry to a g2-file
\arg -dumpMat : Dumps the system matrices to file
\arg -dumpDis : Dumps the solution vector to files, one for each component
\arg -dumpGrid : Dumps a grid file with 8-noded hexahedron elements
\arg -dumpGrid2 : Dumps a grid file with 27-noded hexahedron elements
*/
int main (const int argc, char** argv)
{
SystemMatrix::Type solver = SystemMatrix::SPR;
int nGauss = 4;
int format = 0;
int n[3] = { 2, 2, 2 };
std::vector<int> ignoredPatches;
int iop = 0;
int nev = 10;
int ncv = 20;
double shf = 0.0;
bool free = false;
bool checkRHS = false;
bool vizRHS = false;
bool fixDup = false;
bool dumpGeo = false;
bool dumpMat = false;
bool dumpDis = false;
int dumpGrid = 0;
char* infile = 0;
for (int i = 1; i < argc; i++)
if (!strcmp(argv[i],"-dense"))
solver = SystemMatrix::DENSE;
else if (!strcmp(argv[i],"-superlu"))
solver = SystemMatrix::SPARSE;
else if (!strcmp(argv[i],"-samg"))
solver = SystemMatrix::SAMG;
else if (!strcmp(argv[i],"-swapJ"))
VolumePatch::swapJac = true;
else if (!strcmp(argv[i],"-splineMethod") && i < argc-1)
VolumePatch::splineEvalMethod = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nGauss") && i < argc-1)
nGauss = atoi(argv[++i]);
else if (!strcmp(argv[i],"-vtf") && i < argc-1)
format = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nviz") && i < argc-1)
n[0] = n[1] = n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nu") && i < argc-1)
n[0] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nv") && i < argc-1)
n[1] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nw") && i < argc-1)
n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-ignore"))
while (i < argc-1 && isdigit(argv[i+1][0]))
{
char* endp = 0;
int endVal = 0;
ignoredPatches.push_back(strtol(argv[++i],&endp,10));
if (endp && *endp == ':')
endVal = strtol(endp+1,&endp,10);
while (ignoredPatches.back() < endVal)
ignoredPatches.push_back(ignoredPatches.back()+1);
}
else if (!strcmp(argv[i],"-eig") && i < argc-1)
iop = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nev") && i < argc-1)
nev = atoi(argv[++i]);
else if (!strcmp(argv[i],"-ncv") && i < argc-1)
ncv = atoi(argv[++i]);
else if (!strcmp(argv[i],"-shift"))
shf = atof(argv[++i]);
else if (!strcmp(argv[i],"-free"))
free = true;
else if (!strcmp(argv[i],"-check"))
iop = 100;
else if (!strcmp(argv[i],"-checkRHS"))
checkRHS = true;
else if (!strcmp(argv[i],"-vizRHS"))
vizRHS = true;
else if (!strcmp(argv[i],"-fixDup"))
fixDup = true;
else if (!strcmp(argv[i],"-dumpGeo"))
dumpGeo = true;
else if (!strcmp(argv[i],"-dumpMat"))
dumpMat = true;
else if (!strcmp(argv[i],"-dumpDis"))
dumpDis = true;
else if (!strcmp(argv[i],"-dumpGrid"))
dumpGrid = 1;
else if (!strcmp(argv[i],"-dumpGrid2"))
dumpGrid = 2;
else if (!infile)
infile = argv[i];
if (!infile)
{
std::cout <<"usage: "<< argv[0] <<" <inputfile> [-dense|-superlu|-samg]\n"
<<" [-splineMethod <iop>] [-nGauss <n>] [-vtf <format>]\n"
<<" [-nviz <nviz>] [-nu <nu>] [-nv <nv>] [-nw <nw>]\n"
<<" [-eig <iop>] [-nev <nev>] [-ncv <ncv] [-shift <shf>]\n"
<<" [-free] [-check] [-checkRHS] [-vizRHS] [-fixDup]\n"
<<" [-ignore <p1> <p2> ...] [-swapJ] [-dumpGeo]\n"
<<" [-dumpMat] [-dumpDis] [-dumpGrid|-dumpGrid2]\n";
return 0;
}
// Load vector visualization is not available when using additional viz-points
if (n[0] > 2 || n[1] > 2 || n[2] > 2) vizRHS = false;
// Boundary conditions can be ignored only in generalized eigenvalue analysis
if (abs(iop) != 4 && abs(iop) != 6) free = false;
std::cout <<"\n >>> ICADA Spline FEM prototype <<<"
<<"\n ==================================\n"
<<"\nInput file: "<< infile
<<"\nEquation solver: "<< solver
<<"\nspline evaluation method: "<< VolumePatch::splineEvalMethod
<<"\nNumber of Gauss points: "<< nGauss
<<"\nVTF file format: "<< (format ? "BINARY":"ASCII")
<<"\nNumber of visualization points: "
<< n[0] <<" "<< n[1] <<" "<< n[2];
if (iop != 0 && iop < 100)
std::cout <<"\nEigenproblem solver: "<< iop
<<"\nNumber of eigenvalues: "<< nev
<<"\nNumber of Arnoldi vectors: "<< ncv
<<"\nShift value: "<< shf;
if (free)
std::cout <<"\nSpecified boundary conditions are ignored";
if (VolumePatch::swapJac)
std::cout <<"\nSwapped Jacobian determinant";
if (fixDup)
std::cout <<"\nCo-located nodes will be merged";
if (checkRHS)
std::cout <<"\nChecking that each patch has a right-hand coordinate system";
if (!ignoredPatches.empty())
{
std::cout <<"\nIgnored patches:";
for (size_t i = 0; i < ignoredPatches.size(); i++)
std::cout <<" "<< ignoredPatches[i];
}
std::cout << std::endl;
// Read in model definitions and establish the FE structures
LinearEl model(infile,checkRHS,free);
if (!model.preprocess(ignoredPatches,fixDup))
return 1;
if (dumpGeo) // Write the refined geometry to g2-file
model.dumpGeometry(strcat(strtok(infile,"."),".g2"));
Result solution;
Matrix eNorm;
switch (iop) {
case 0:
case 5:
case 54:
case 56:
// Assemble and solve the linear system of equations
if (!model.assembleKandR(solver,nGauss,vizRHS ? &solution.load : 0))
return 2;
else if (!model.solve(solution.displ))
return 3;
else if (iop == 0)
if (!model.solutionNorms(solution.norms,nGauss,solution.displ))
return 4;
else
break;
// Linearized buckling: Assemble stiffness matrices [K] and [Kg]
if (!model.assembleKandKg(solver,nGauss,solution.displ))
return 2;
else if (dumpMat)
model.dumpMat("Km.mat","Kg.mat");
else if (!model.modes(iop > 50 ? iop%50 : iop, nev,ncv,shf,solution.modes))
return 3;
break;
case -1:
case -2:
case 1:
case 2:
// Assemble and solve the regular eigenvalue problem
if (!model.assembleKonly(solver,nGauss))
return 2;
else if (!model.modes(iop,nev,ncv,shf,solution.modes))
return 3;
break;
case 100:
case 101:
case 102:
break; // Model check
default:
// Assemble and solve the generalized eigenvalue problem
solution.freq = true;
if (!model.assembleKandM(solver,nGauss))
return 2;
else if (!model.modes(iop,nev,ncv,shf,solution.modes))
return 3;
}
// Write VTF-file with model and results
if (dumpGrid == 0 || iop < 100)
if (!model.writeGlv(infile,solution,n,format,iop==100)) return 4;
// Write a global FE grid file for external processing
if (dumpGrid > 0)
if (!model.writeGlobalGrid(infile, n, dumpGrid == 2 ? 27 : 8)) return 5;
// Write solution vectors to files for external viewing
if (dumpDis)
{
strcat(strtok(infile,"."),"_dis");
if (!solution.displ.empty())
model.dumpSolution(infile,solution.displ);
infile[strlen(infile)-3] = 'm';
infile[strlen(infile)-2] = '0';
infile[strlen(infile)-1] = '\0';
for (size_t j = 1; j < solution.modes.size() && j < 10; j++)
{
infile[strlen(infile)-1]++;
model.dumpSolution(infile,solution.modes[j-1].eigVec);
}
}
return 0;
}