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

Added usage of adaptive driver in linear elasticity solver. Some corrections on VTF export of adaptive simulation.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1205 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-09-25 13:27:20 +00:00 committed by Knut Morten Okstad
parent 5143525343
commit ba3d3a0fe5
12 changed files with 266 additions and 157 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
Apps/LinearElasticity/main_LinEl3D.C
View File

@ -13,6 +13,7 @@
#include "SIMLinEl3D.h"
#include "SIMLinElKL.h"
#include "AdaptiveSIM.h"
#include "LinAlgInit.h"
#include "HDF5Writer.h"
#include "XMLWriter.h"
@ -38,6 +39,7 @@
\arg -petsc : Use equation solver from PETSc library
\arg -lag : Use Lagrangian basis functions instead of splines/NURBS
\arg -spec : Use Spectral basis functions instead of splines/NURBS
\arg -LR : Use LR-spline basis functions instead of tensorial splines/NURBS
\arg -nGauss \a n : Number of Gauss points over a knot-span in each direction
\arg -vtf \a format : VTF-file format (-1=NONE, 0=ASCII, 1=BINARY)
\arg -nviz \a nviz : Number of visualization points over each knot-span
@ -60,6 +62,7 @@
\arg -2Dpstrain : Use two-parametric simulation driver (plane strain)
\arg -2Daxisymm : Use two-parametric simulation driver (axi-symmetric solid)
\arg -KL : Use two-parametric simulation driver for Kirchhoff-Love plate
\arg -adap : Use adaptive simulation driver with LR-splines discretization
*/
int main (int argc, char** argv)
@ -106,6 +109,8 @@ int main (int argc, char** argv)
SIMbase::discretization = SIMbase::Lagrange;
else if (!strncmp(argv[i],"-spec",5))
SIMbase::discretization = SIMbase::Spectral;
else if (!strncmp(argv[i],"-LR",3))
SIMbase::discretization = SIMbase::LRSpline;
else if (!strcmp(argv[i],"-nGauss") && i < argc-1)
nGauss = atoi(argv[++i]);
else if (!strcmp(argv[i],"-vtf") && i < argc-1)
@ -151,10 +156,11 @@ int main (int argc, char** argv)
twoD = SIMLinEl2D::axiSymmetry = true;
else if (!strncmp(argv[i],"-2D",3))
twoD = true;
else if (!strncmp(argv[i],"-lag",4))
SIMbase::discretization = SIMbase::Lagrange;
else if (!strncmp(argv[i],"-spec",5))
SIMbase::discretization = SIMbase::Spectral;
else if (!strncmp(argv[i],"-adap",5))
{
SIMbase::discretization = SIMbase::LRSpline;
iop = 10;
}
else if (!infile)
infile = argv[i];
else
@ -164,8 +170,8 @@ int main (int argc, char** argv)
{
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n "
<<" [-free] [-lag] [-spec] [-2D[pstrain|axisymm]|-KL]"
<<" [-nGauss <n>]\n [-vtf <format> [-nviz <nviz>]"
<<" [-free] [-lag] [-spec] [-LR] [-2D[pstrain|axisymm]|-KL]"
<<" [-adap] [-nGauss <n>]\n [-vtf <format> [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]] [-hdf5]\n"
<<" [-eig <iop> [-nev <nev>] [-ncv <ncv] [-shift <shf>]]\n"
<<" [-ignore <p1> <p2> ...] [-fixDup]"
@ -197,7 +203,7 @@ int main (int argc, char** argv)
if (!twoD) std::cout <<" "<< n[2];
}
if (iop > 0 && iop < 100)
if (iop > 0 && iop < 10)
std::cout <<"\nEigenproblem solver: "<< iop
<<"\nNumber of eigenvalues: "<< nev
<<"\nNumber of Arnoldi vectors: "<< ncv
@ -232,7 +238,12 @@ int main (int argc, char** argv)
else
model = new SIMLinEl3D(checkRHS);
if (!model->read(infile))
SIMinput* theSim = model;
AdaptiveSIM* aSim = 0;
if (iop == 10)
theSim = aSim = new AdaptiveSIM(model);
if (!theSim->read(infile))
return 1;
utl::profiler->stop("Model input");
@ -250,6 +261,7 @@ int main (int argc, char** argv)
Vectors displ(1), projs;
std::vector<Mode> modes;
std::vector<Mode>::const_iterator it;
int iStep = 1, nBlock = 0;
switch (iop) {
case 0:
@ -331,6 +343,16 @@ int main (int argc, char** argv)
return 6;
break;
case 10:
// Adaptive simulation
while (true)
if (!aSim->solveStep(infile,solver,iStep))
return 5;
else if (!aSim->writeGlv(infile,format,n,iStep,nBlock))
return 6;
else if (!aSim->adaptMesh(++iStep))
break;
case 100:
break; // Model check
@ -353,22 +375,22 @@ int main (int argc, char** argv)
if (linalg.myPid == 0)
std::cout <<"\nWriting HDF5 file "<< infile <<".hdf5"<< std::endl;
DataExporter exporter(true);
exporter.registerField("u","solution",DataExporter::SIM,DataExporter::SECONDARY);
exporter.registerField("u","solution",DataExporter::SIM,
DataExporter::SECONDARY);
exporter.setFieldValue("u",model,&displ.front());
exporter.registerWriter(new HDF5Writer(infile));
exporter.registerWriter(new XMLWriter(infile));
exporter.dumpTimeLevel();
}
if (format >= 0)
if (iop != 10 && format >= 0)
{
// Write VTF-file with model geometry
int iStep = 1, nBlock = 0;
if (!model->writeGlv(infile,n,format))
return 7;
// Write boundary tractions, if any
if (!model->writeGlvT(iStep,nBlock))
if (!model->writeGlvT(iStep,++nBlock))
return 8;
// Write Dirichlet boundary conditions
@ -386,20 +408,20 @@ int main (int argc, char** argv)
// Write projected solution fields to VTF-file
if (projs.size() > 0)
if (!model->writeGlvP(projs.front(),n,iStep,nBlock))
return 10;
return 11;
// Write eigenmodes
for (it = modes.begin(); it != modes.end(); it++)
if (!model->writeGlvM(*it, iop==3 || iop==4 || iop==6, n, nBlock))
return 11;
return 12;
// Write element norms
if (!model->writeGlvN(eNorm,iStep,nBlock))
return 12;
return 13;
model->writeGlvStep(1);
model->closeGlv();
}
model->closeGlv();
if (dumpASCII)
{
@ -431,6 +453,6 @@ int main (int argc, char** argv)
}
utl::profiler->stop("Postprocessing");
delete model;
delete theSim;
return 0;
}

28
src/Integrands/KirchhoffLovePlate.C
View File

@ -12,14 +12,15 @@
//==============================================================================
#include "KirchhoffLovePlate.h"
#include "FiniteElement.h"
#include "LinIsotropic.h"
#include "FiniteElement.h"
#include "Utilities.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3.h"
#include "Vec3Oper.h"
#include "AnaSol.h"
#include "VTF.h"
KirchhoffLovePlate::KirchhoffLovePlate ()
@ -77,6 +78,7 @@ void KirchhoffLovePlate::setMode (SIM::SolutionMode mode)
myMats->resize(1,1);
eK = &myMats->A[0];
eS = &myMats->b[0];
presVal.clear();
break;
case SIM::DYNAMIC:
@ -84,6 +86,7 @@ void KirchhoffLovePlate::setMode (SIM::SolutionMode mode)
eK = &myMats->A[0];
eM = &myMats->A[1];
eS = &myMats->b[0];
presVal.clear();
break;
case SIM::VIBRATION:
@ -107,6 +110,7 @@ void KirchhoffLovePlate::setMode (SIM::SolutionMode mode)
if (myMats->b.empty())
myMats->b.resize(1);
eS = &myMats->b[0];
presVal.clear();
break;
case SIM::RECOVERY:
@ -118,6 +122,7 @@ void KirchhoffLovePlate::setMode (SIM::SolutionMode mode)
default:
myMats->resize(0,0);
mySols.clear();
presVal.clear();
}
}
@ -170,6 +175,21 @@ bool KirchhoffLovePlate::haveLoads () const
}
bool KirchhoffLovePlate::writeGlvT (VTF* vtf, int iStep, int& nBlock) const
{
if (presVal.empty())
return true;
else if (!vtf)
return false;
// Write surface pressures as discrete point vectors to the VTF-file
if (!vtf->writeVectors(presVal,++nBlock))
return false;
return vtf->writeVblk(nBlock,"Pressure",1,iStep);
}
/*!
The strain-displacement matrix for a Kirchhoff-Love plate element is formally
defined as:
@ -236,9 +256,11 @@ void KirchhoffLovePlate::formBodyForce (Vector& ES, const Vector& N,
const Vec3& X, double detJW) const
{
double p = this->getPressure(X);
if (p != 0.0)
{
ES.add(N,p*detJW);
presVal[X].z = p; // Store pressure value for visualization
}
}

17
src/Integrands/KirchhoffLovePlate.h
View File

@ -15,6 +15,8 @@
#define _KIRCHHOFF_LOVE_PLATE_H
#include "IntegrandBase.h"
#include "Vec3.h"
#include <map>
class LocalSystem;
class Material;
@ -31,7 +33,6 @@ class KirchhoffLovePlate : public IntegrandBase
public:
//! \brief The default constructor initializes all pointers to zero.
KirchhoffLovePlate();
//! \brief The destructor frees the dynamically allocated data objects.
virtual ~KirchhoffLovePlate();
@ -84,12 +85,10 @@ public:
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] d2NdX2 Basis function 2nd derivatives at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
virtual bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
virtual bool evalSol(Vector& s, const Vector&, const Matrix&,
const Matrix3D& d2NdX2, const Vec3& X,
const std::vector<int>& MNPC) const;
@ -115,6 +114,14 @@ public:
//! \brief Returns whether an external load is defined.
virtual bool haveLoads() const;
//! \brief Writes the surface pressure for a given time step to VTF-file.
//! \param vtf The VTF-file object to receive the pressure vectors
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter
virtual bool writeGlvT(VTF* vtf, int iStep, int& nBlock) const;
//! \brief Returns whether there are any pressure values to write to VTF.
virtual bool hasTractionValues() const { return !presVal.empty(); }
//! \brief Returns a pointer to an Integrand for solution norm evaluation.
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
@ -177,6 +184,8 @@ protected:
LocalSystem* locSys; //!< Local coordinate system for result output
RealFunc* presFld; //!< Pointer to pressure field
mutable std::map<Vec3,Vec3> presVal; //!< Pressure field point values
// Work arrays declared as members to avoid frequent re-allocation
// within the numerical integration loop (for reduced overhead)
mutable Matrix Bmat; //!< Strain-displacement matrix

34
src/SIM/AdaptiveSIM.C
View File

@ -143,9 +143,39 @@ void AdaptiveSIM::printNorms (const Vector& norms, std::ostream& os)
os <<"Energy norm |u^h| = a(u^h,u^h)^0.5 : "<< norms(1)
<<"\nExternal energy ((h,u^h)+(t,u^h)^0.5 : "<< norms(2);
if (norms.size() > 2)
os <<"\nExact norm |u| = a(u,u)^0.5 : "<< norms(3);
os <<"\nExact norm |u| = a(u,u)^0.5 : "<< norms(3);
if (norms.size() > 3)
os <<"\nExact error a(e,e)^0.5, e=u-u^h : "<< norms(4)
os <<"\nExact error a(e,e)^0.5, e=u-u^h : "<< norms(4)
<<"\nExact relative error (%) : "<< 100.0*norms(4)/norms(3);
os << std::endl;
}
bool AdaptiveSIM::writeGlv (const char* infile, int format, const int* nViz,
int iStep, int& nBlock)
{
if (format < 0) return true;
// Write VTF-file with model geometry
if (!model->writeGlvG(nViz, nBlock, iStep == 1 ? infile : 0, format))
return false;
// Write boundary tractions, if any
if (!model->writeGlvT(iStep,nBlock))
return false;
// Write Dirichlet boundary conditions
if (!model->writeGlvBC(nViz,nBlock,iStep))
return false;
// Write solution fields
if (!model->writeGlvS(linsol,nViz,iStep,nBlock))
return false;
// Write element norms
if (!model->writeGlvN(eNorm,iStep,nBlock))
return false;
// Write state information
return model->writeGlvStep(iStep,iStep,1);
}

16
src/SIM/AdaptiveSIM.h
View File

@ -45,13 +45,17 @@ public:
//! \param[in] iStep Refinement step counter
bool adaptMesh(int iStep);
//! \brief Prints out the global norms to given stream
static void printNorms(const Vector& norms, std::ostream& os);
//! \brief Writes current mesh and results to the VTF-file.
//! \param[in] infile File name used to construct the VTF-file name from
//! \param[in] format Format of VTF-file (0=ASCII, 1=BINARY)
//! \param[in] nViz Number of visualization points over a knot-span
//! \param[in] iStep Refinement step identifier
//! \param nBlock Running result block counter
bool writeGlv(const char* infile, int format, const int* nViz,
int iStep, int& nBlock);
//! \brief Returns the current primary solution vector.
const Vector& getSolution() const { return linsol; }
//! \brief Returns the current element norms.
const Matrix& getElementNorms() const { return eNorm; }
//! \brief Prints out the global norms to given stream.
static void printNorms(const Vector& norms, std::ostream& os);
protected:

16
src/SIM/SIM1D.C
View File

@ -133,13 +133,9 @@ bool SIM1D::parse (char* keyWord, std::istream& is)
else
{
this->setPropertyType(code,Property::DIRICHLET_INHOM);
if ((cline = strtok(NULL," ")))
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,d));
else
{
std::cout << d;
myScalars[code] = new ConstFunc(d);
}
cline = strtok(NULL," ");
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,d));
}
std::cout << std::endl;
}
@ -300,10 +296,8 @@ bool SIM1D::parse (char* keyWord, std::istream& is)
if (!this->addConstraint(patch,pvert,0,bcode%1000,code))
return false;
if ((cline = strtok(NULL," ")))
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,pd));
else
myScalars[code] = new ConstFunc(pd);
cline = strtok(NULL," ");
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,pd));
}
std::cout << std::endl;
}

16
src/SIM/SIM3D.C
View File

@ -196,13 +196,9 @@ bool SIM3D::parse (char* keyWord, std::istream& is)
else
{
this->setPropertyType(code,Property::DIRICHLET_INHOM);
if ((cline = strtok(NULL," ")))
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,d));
else
{
std::cout << d;
myScalars[code] = new ConstFunc(d);
}
cline = strtok(NULL," ");
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,d));
}
std::cout << std::endl;
}
@ -430,10 +426,8 @@ bool SIM3D::parse (char* keyWord, std::istream& is)
if (!this->addConstraint(patch,abs(pface),ldim,bcode%1000,code))
return false;
if ((cline = strtok(NULL," ")))
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,pd));
else
myScalars[code] = new ConstFunc(pd);
cline = strtok(NULL," ");
myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(cline,pd));
}
if (pface < 10) std::cout << std::endl;
}

71
src/SIM/SIMbase.C
View File

@ -950,35 +950,37 @@ bool SIMbase::systemModes (std::vector<Mode>& solution,
bool SIMbase::writeGlv (const char* inpFile, const int* nViz, int format)
{
if (myVtf) return false;
#if HAS_VTFAPI == 2
const char* ext = ".vtfx";
#else
const char* ext = ".vtf";
#endif
// Open a new VTF-file
char* vtfName = new char[strlen(inpFile)+10];
strtok(strcpy(vtfName,inpFile),".");
if (nProc > 1)
sprintf(vtfName+strlen(vtfName),"_p%04d%s",myPid,ext);
else
strcat(vtfName,ext);
std::cout <<"\nWriting VTF-file "<< vtfName << std::endl;
myVtf = new VTF(vtfName,format);
delete[] vtfName;
// Dirty hack to be backwards compatible! (block 0 is usually unused)
int bar=0;
return writeGlvG(1,nViz,bar);
int nBlock = 0;
return this->writeGlvG(nViz,nBlock,inpFile,format);
}
bool SIMbase::writeGlvG(int iStep, const int* nViz, int& nBlock)
bool SIMbase::writeGlvG (const int* nViz, int& nBlock,
const char* inpFile, int format)
{
myVtf->clearGeometryBlocks();
if (inpFile)
{
if (myVtf) return false;
#if HAS_VTFAPI == 2
const char* ext = ".vtfx";
#else
const char* ext = ".vtf";
#endif
// Open a new VTF-file
char* vtfName = new char[strlen(inpFile)+10];
strtok(strcpy(vtfName,inpFile),".");
if (nProc > 1)
sprintf(vtfName+strlen(vtfName),"_p%04d%s",myPid,ext);
else
strcat(vtfName,ext);
std::cout <<"\nWriting VTF-file "<< vtfName << std::endl;
myVtf = new VTF(vtfName,format);
delete[] vtfName;
}
// Convert and write model geometry
char pname[16];
for (size_t i = 0; i < myModel.size(); i++)
@ -995,14 +997,15 @@ bool SIMbase::writeGlvG(int iStep, const int* nViz, int& nBlock)
sprintf(pname,"Patch %ld",i+1);
if (!myVtf->writeGrid(lvb,pname,++nBlock))
return false;
myVtf->writeGeometryBlocks(iStep);
}
// Do not write the geometry blocks to file yet, writeVectors might create
// an additional block for the point vectors, see method VTF::writeVectors
return true;
}
bool SIMbase::writeGlvBC (const int* nViz, int& nBlock) const
bool SIMbase::writeGlvBC (const int* nViz, int& nBlock, int iStep) const
{
if (!myVtf) return false;
@ -1058,7 +1061,7 @@ bool SIMbase::writeGlvBC (const int* nViz, int& nBlock) const
for (j = 0; j < 3; j++)
if (!dID[j].empty())
if (!myVtf->writeSblk(dID[j],label[j],1+j))
if (!myVtf->writeSblk(dID[j],label[j],1+j,iStep))
return false;
return true;
@ -1332,6 +1335,8 @@ bool SIMbase::writeGlvP (const Vector& ssol, const int* nViz,
bool SIMbase::writeGlvStep (int iStep, double value, int itype)
{
myVtf->writeGeometryBlocks(iStep);
if (itype == 0)
return myVtf->writeState(iStep,"Time %g",value,itype);
else
@ -1413,7 +1418,7 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
}
for (j = k = 0; j < field.rows() && j < 10; j++)
if (NormBase::hasElementContributions(j))
if (j != 1) // Norm #1 (external norm) does not have element contributions
if (!myVtf->writeEres(field.getRow(1+j),++nBlock,geomID))
return false;
else
@ -1421,13 +1426,9 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
}
int idBlock = 200;
for (j = k = 0; !sID[j].empty(); j++, k++)
{
while (!NormBase::hasElementContributions(k))
k++;
if (!myVtf->writeSblk(sID[j],NormBase::getName(k),++idBlock,iStep,true))
for (k = 0; !sID[k].empty(); k++)
if (!myVtf->writeSblk(sID[k],NormBase::getName(k),++idBlock,iStep,true))
return false;
}
return true;
}

11
src/SIM/SIMbase.h
View File

@ -324,16 +324,19 @@ public:
//! by the other \a writeGlv* methods.
bool writeGlv(const char* inpFile, const int* nViz, int format);
//! \brief Write current geometry to the VTF-file
//! \param[in] iStep Load/time step identifier
//! \brief Write current model geometry to the VTF-file.
//! \param[in] nViz Number of visualization points over each knot-span
//! \param nBlock Running result block counter
bool writeGlvG(int iStep, const int* nViz, int& nBlock);
//! \param[in] inpFile File name used to construct the VTF-file name from
//! \param[in] format Format of VTF-file (0=ASCII, 1=BINARY)
bool writeGlvG(const int* nViz, int& nBlock,
const char* inpFile = 0, int format = -1);
//! \brief Writes boundary conditions as scalar fields to the VTF-file.
//! \param[in] nViz Number of visualization points over each knot-span
//! \param nBlock Running result block counter
bool writeGlvBC(const int* nViz, int& nBlock) const;
//! \param[in] iStep Load/time step identifier
bool writeGlvBC(const int* nViz, int& nBlock, int iStep = 1) const;
//! \brief Writes boundary tractions for a given time step to the VTF-file.
//! \param[in] iStep Load/time step identifier

35
src/Utility/Functions.C
View File

@ -177,29 +177,54 @@ const RealFunc* utl::parseRealFunc (char* cline, real A)
switch (linear) {
case 1:
f = new LinearXFunc(A,atof(cline));
cline = strtok(NULL," ");
break;
case 2:
f = new LinearYFunc(A,atof(cline));
cline = strtok(NULL," ");
break;
case 3:
f = new LinearZFunc(A,atof(cline));
cline = strtok(NULL," ");
break;
case 4:
f = new LinearRotZFunc(true,A,atof(cline),atof(strtok(NULL," ")));
cline = strtok(NULL," ");
break;
case 5:
f = new LinearRotZFunc(false,A,atof(cline),atof(strtok(NULL," ")));
cline = strtok(NULL," ");
break;
case 6:
std::cout <<"StepX("<< cline <<"))";
f = new StepXFunc(A,atof(cline),atof(strtok(NULL," ")));
{
double x0 = atof(cline);
double x1 = atof(strtok(NULL," "));
std::cout <<"StepX("<< x0 <<","<< x1 <<"))";
f = new StepXFunc(A,x0,x1);
}
cline = strtok(NULL," ");
break;
case 7:
std::cout <<"StepXY("<< cline <<"))";
f = new StepXYFunc(A,atof(cline),atof(strtok(NULL," ")));
{
double x0 = atof(cline);
double y0 = atof(strtok(NULL," "));
cline = strtok(NULL," ");
if (cline && cline[0] == 't')
{
double x1 = atof(strtok(NULL," "));
double y1 = atof(strtok(NULL," "));
std::cout <<"StepXY(["<< x0<<","<<x1 <<"]x["<< y0<<","<<y1 <<"]))";
f = new StepXYFunc(A,x1,y1,x0,y0);
cline = strtok(NULL," ");
}
else
{
std::cout <<"StepXY([-inf,"<< x0 <<"]x[-inf,"<< y0 <<"]))";
f = new StepXYFunc(A,x0,y0);
}
}
break;
}
cline = strtok(NULL," ");
}
else if (quadratic > 0 && (cline = strtok(NULL," ")))
{

105
src/Utility/VTF.C
View File

@ -23,7 +23,7 @@
#endif
#include "ElementBlock.h"
#include <iostream>
#include <stdio.h>
#include <cstdio>
real VTF::vecOffset[3] = { 0.0, 0.0, 0.0 };
@ -32,14 +32,14 @@ real VTF::vecOffset[3] = { 0.0, 0.0, 0.0 };
VTF::VTF (const char* filename, int type)
{
myFile = 0;
geoBlock = new VTFAGeometryBlock;
myState = 0;
myGBlock = 0;
pointGeoID = 0;
if (!filename) return;
#if HAS_VTFAPI == 1
// Create the VTF file object
myFile = new VTFAFile();
// Enable debug info to stderr/console
myFile->SetOutputDebugError(1);
@ -74,10 +74,15 @@ VTF::~VTF ()
{
if (!myFile) return;
size_t i;
#if HAS_VTFAPI == 1
if (VTFA_FAILURE(myFile->WriteBlock(geoBlock)))
showError("Error writing Geometry Block");
if (myGBlock)
{
if (VTFA_FAILURE(myFile->WriteBlock(myGBlock)))
showError("Error writing Geometry Block");
delete myGBlock;
}
size_t i;
for (i = 0; i < myDBlock.size(); i++)
if (myDBlock[i])
{
@ -113,8 +118,14 @@ VTF::~VTF ()
delete myFile;
#elif HAS_VTFAPI == 2
if (VTFA_FAILURE(myDatabase->WriteBlock(geoBlock)))
showError("Error writing Geometry Block");
if (myGBlock)
{
if (VTFA_FAILURE(myDatabase->WriteBlock(myGBlock)))
showError("Error writing Geometry Block");
delete myGBlock;
}
size_t i;
for (i = 0; i < myDBlock.size(); i++)
if (myDBlock[i])
{
@ -125,16 +136,16 @@ VTF::~VTF ()
for (i = 0; i < myVBlock.size(); i++)
if (myVBlock[i])
{
if (VTFA_FAILURE(myDatabase->WriteBlock(myVBlock[i])))
showError("Error writing Vector Block");
delete myVBlock[i];
if (VTFA_FAILURE(myDatabase->WriteBlock(myVBlock[i])))
showError("Error writing Vector Block");
delete myVBlock[i];
}
for (i = 0; i < mySBlock.size(); i++)
if (mySBlock[i])
{
if (VTFA_FAILURE(myDatabase->WriteBlock(mySBlock[i])))
showError("Error writing Scalar Block");
showError("Error writing Scalar Block");
delete mySBlock[i];
}
@ -148,15 +159,14 @@ VTF::~VTF ()
VTFXACase* singleCase = new VTFXACase(myFile,"Single case",1,1);
VTFXACasePropertiesBlock frameGeneratorProps(VT_CT_FRAME_GENERATOR_SETTINGS);
frameGeneratorProps.AddInt(VT_PI_FG_FEM_MODEL_IDS, 1); // for VTFx just use always "1" here
frameGeneratorProps.AddInt(VT_PI_FG_FEM_MODEL_IDS,1);
singleCase->WritePropertiesBlock(&frameGeneratorProps);
for (i = 0; i < myBlocks.size(); i++) {
VTFXACasePropertiesBlock partAttr(VT_CT_PART_ATTRIBUTES);
partAttr.SetPartID(i+1);
// Turn on mesh
partAttr.AddBool(VT_PB_PA_MESH, VTFXA_FALSE);
partAttr.AddBool(VT_PB_PA_DISPLACEMENTS, VTFXA_FALSE);
partAttr.AddBool(VT_PB_PA_DISPLACEMENTS, VTFXA_FALSE);
singleCase->WritePropertiesBlock(&partAttr);
}
if (VTFA_FAILURE(myFile->CloseFile()))
@ -166,38 +176,36 @@ VTF::~VTF ()
delete myDatabase;
delete myFile;
#endif
delete geoBlock;
}
void VTF::writeGeometryBlocks(int iStep)
void VTF::writeGeometryBlocks (int iStep)
{
if (!myBlocks.size())
pointGeoID = 0;
if (myBlocks.empty())
return;
size_t i;
std::vector<int> geomID(myBlocks.size());
for (i = 0; i < myBlocks.size(); i++)
for (size_t i = 0; i < myBlocks.size(); i++)
{
geomID[i] = myBlocks[i].first;
delete myBlocks[i].second;
}
myBlocks.clear();
#ifdef HAS_VTFAPI
if (!myGBlock) myGBlock = new VTFAGeometryBlock();
#if HAS_VTFAPI == 1
geoBlock->SetGeometryElementBlocks(&geomID.front(),geomID.size(),iStep);
myGBlock->SetGeometryElementBlocks(&geomID.front(),geomID.size(),iStep);
#elif HAS_VTFAPI == 2
std::cout << "yo " << iStep << std::endl;
geoBlock->SetElementBlocksForState(&geomID.front(),geomID.size(),iStep);
myGBlock->SetElementBlocksForState(&geomID.front(),geomID.size(),iStep);
#endif
#endif
}
void VTF::clearGeometryBlocks()
{
for (size_t i=0;i<myBlocks.size();++i)
delete myBlocks[i].second;
myBlocks.clear();
}
bool VTF::writeGrid (const ElementBlock* block, const char* partname, int nBlock)
bool VTF::writeGrid (const ElementBlock* block, const char* partname,
int nBlock)
{
if (!myFile) return true;
@ -345,26 +353,25 @@ bool VTF::writeVectors (const std::map<Vec3,Vec3>& pntResult, int idBlock)
{
#if HAS_VTFAPI == 1
bool writePoints = false;
static int geomID = 0;
if (geomID == 0)
if (pointGeoID == 0)
{
// The big assumption here is that we have only one call to writeVectors
// per time step, and that all subsequent calls are with the same points
myBlocks.push_back(std::make_pair(myBlocks.size()+1,new ElementBlock()));
geomID = myBlocks.size();
pointGeoID = idBlock;
myBlocks.push_back(std::make_pair(pointGeoID,new ElementBlock()));
writePoints = true;
}
VTFANodeBlock nBlock(geomID,0);
VTFANodeBlock nBlock(pointGeoID,0);
VTFAResultBlock rBlock(idBlock,VTFA_DIM_VECTOR,VTFA_RESMAP_NODE,0);
size_t i = 0, np = pntResult.size();
if (writePoints && VTFA_FAILURE(nBlock.SetNumNodes(np)))
return showError("Error defining node block",geomID);
return showError("Error defining node block",pointGeoID);
else if (VTFA_FAILURE(rBlock.SetNumResults(np)))
return showError("Error defining result block",idBlock);
else
rBlock.SetMapToBlockID(geomID);
rBlock.SetMapToBlockID(pointGeoID);
int* mnpc = writePoints ? new int[np] : 0;
std::map<Vec3,Vec3>::const_iterator cit;
@ -372,7 +379,7 @@ bool VTF::writeVectors (const std::map<Vec3,Vec3>& pntResult, int idBlock)
if (writePoints && VTFA_FAILURE(nBlock.AddNode(vecOffset[0]+cit->first.x,
vecOffset[1]+cit->first.y,
vecOffset[2]+cit->first.z)))
return showError("Error adding node to block",geomID);
return showError("Error adding node to block",pointGeoID);
else if (VTFA_FAILURE(rBlock.AddResult(cit->second.x,
cit->second.y,
cit->second.z)))
@ -384,17 +391,17 @@ bool VTF::writeVectors (const std::map<Vec3,Vec3>& pntResult, int idBlock)
{
// We must define an element block (with point elements) also,
// otherwise GLview does not visualize the vectors
VTFAElementBlock eBlock(geomID,0,0);
eBlock.SetPartID(geomID);
eBlock.SetNodeBlockID(geomID);
VTFAElementBlock eBlock(pointGeoID,0,0);
eBlock.SetPartID(pointGeoID);
eBlock.SetNodeBlockID(pointGeoID);
if (VTFA_FAILURE(eBlock.AddElements(VTFA_POINTS,mnpc,np)))
return showError("Error defining element block",geomID);
return showError("Error defining element block",pointGeoID);
delete[] mnpc;
if (VTFA_FAILURE(myFile->WriteBlock(&nBlock)))
return showError("Error writing node block",geomID);
return showError("Error writing node block",pointGeoID);
else if (VTFA_FAILURE(myFile->WriteBlock(&eBlock)))
return showError("Error writing element block",geomID);
return showError("Error writing element block",pointGeoID);
}
if (VTFA_FAILURE(myFile->WriteBlock(&rBlock)))
return showError("Error writing result block",idBlock);
@ -492,7 +499,7 @@ bool VTF::writeVblk (const std::vector<int>& vBlockIDs, const char* resultName,
}
myVBlock[idBlock]->SetResultID(idBlock);
if (VTFA_FAILURE(myVBlock[idBlock]->SetResultValuesBlocks(&vBlockIDs.front(),
vBlockIDs.size(),
vBlockIDs.size(),
iStep)))
return showError("Error defining vector block",idBlock);
#endif
@ -534,7 +541,7 @@ bool VTF::writeSblk (int sBlockID, const char* resultName,
bool VTF::writeSblk (const std::vector<int>& sBlockIDs, const char* resultName,
int idBlock, int iStep, bool elementData)
int idBlock, int iStep, bool elementData)
{
if ((int)mySBlock.size() < idBlock) mySBlock.resize(idBlock,0);
@ -552,7 +559,7 @@ bool VTF::writeSblk (const std::vector<int>& sBlockIDs, const char* resultName,
{
int type = elementData?VTFXA_RESMAP_ELEMENT:VTFXA_RESMAP_NODE;
mySBlock[idBlock] = new VTFXAResultBlock(idBlock+1,
VTFXA_RESTYPE_SCALAR,type);
VTFXA_RESTYPE_SCALAR,type);
if (resultName) mySBlock[idBlock]->SetName(resultName);
}
mySBlock[idBlock]->SetResultID(idBlock);

18
src/Utility/VTF.h
View File

@ -61,7 +61,7 @@ public:
//! \param[in] partname Name of the geometry being written
//! \param[in] iStep Load/Time step identifier
virtual bool writeGrid(const ElementBlock* g, const char* partname,
int iStep=1);
int iStep = 1);
//! \brief Writes a block of scalar nodal results to the VTF-file.
//! \param[in] nodeResult Vector of nodal results,
@ -100,7 +100,7 @@ public:
//! \param[in] iStep Load/Time step identifier
//! \param[in] elementData false -> data per node, true -> data per element
bool writeSblk(int sBlockID, const char* resultName = 0, int idBlock = 1,
int iStep = 1, bool elementData=false);
int iStep = 1, bool elementData = false);
//! \brief Writes a scalar block definition to the VTF-file.
//! \param[in] sBlockIDs All result blocks that make up this scalar block
//! \param[in] resultName Name of the result quantity
@ -108,8 +108,8 @@ public:
//! \param[in] iStep Load/Time step identifier
//! \param[in] elementData false -> data per node, true -> data per element
virtual bool writeSblk(const std::vector<int>& sBlockIDs,
const char* resultName = 0,
int idBlock = 1, int iStep = 1, bool elementData=false);
const char* resultName = 0, int idBlock = 1,
int iStep = 1, bool elementData = false);
//! \brief Writes a vector block definition to the VTF-file.
//! \param[in] vBlockID The result block that makes up this vector block
//! \param[in] resultName Name of the result quantity
@ -145,12 +145,9 @@ public:
//! \brief Returns the pointer to a geometry block.
virtual const ElementBlock* getBlock(int geomID) const { return myBlocks[geomID-1].second; }
//! \brief Add the current FE geometry blocks to the description block
//! \brief Adds the current FE geometry blocks to the description block.
void writeGeometryBlocks(int iStep);
//! \brief Drop current FE geometry blocks
void clearGeometryBlocks();
private:
//! \brief Writes a node block to the VTF-file.
//! \details The coordinates of the last added element block are written.
@ -177,11 +174,12 @@ private:
VTFXADatabase* myDatabase; //!< Pointer to VTFx database object for this file
#endif
VTFAStateInfoBlock* myState; //!< The state info block for this file
VTFAGeometryBlock* myGBlock; //!< The geometry description block for this file
std::vector<VTFADisplacementBlock*> myDBlock; //!< Displacement blocks
std::vector<VTFAVectorBlock*> myVBlock; //!< Vector field blocks
std::vector<VTFAScalarBlock*> mySBlock; //!< Scalar field blocks
std::vector< std::pair<int,const ElementBlock*> > myBlocks; //!< The FE geometry
VTFAGeometryBlock* geoBlock;
std::vector< std::pair<int,const ElementBlock*> > myBlocks; //!< FE geometry
int pointGeoID; //!< ID of current point vector geometry block
};
#endif