diff --git a/Apps/LinearElasticity/main_LinEl3D.C b/Apps/LinearElasticity/main_LinEl3D.C index d1f40ed3..15af3b22 100644 --- a/Apps/LinearElasticity/main_LinEl3D.C +++ b/Apps/LinearElasticity/main_LinEl3D.C @@ -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] <<" [-dense|-spr|-superlu[]|-samg|-petsc]\n " - <<" [-free] [-lag] [-spec] [-2D[pstrain|axisymm]|-KL]" - <<" [-nGauss ]\n [-vtf [-nviz ]" + <<" [-free] [-lag] [-spec] [-LR] [-2D[pstrain|axisymm]|-KL]" + <<" [-adap] [-nGauss ]\n [-vtf [-nviz ]" <<" [-nu ] [-nv ] [-nw ]] [-hdf5]\n" <<" [-eig [-nev ] [-ncv ]]\n" <<" [-ignore ...] [-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 modes; std::vector::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; } diff --git a/src/Integrands/KirchhoffLovePlate.C b/src/Integrands/KirchhoffLovePlate.C index 8bf85316..6be5aa70 100644 --- a/src/Integrands/KirchhoffLovePlate.C +++ b/src/Integrands/KirchhoffLovePlate.C @@ -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 + } } diff --git a/src/Integrands/KirchhoffLovePlate.h b/src/Integrands/KirchhoffLovePlate.h index 9f4a1923..1dc7b221 100644 --- a/src/Integrands/KirchhoffLovePlate.h +++ b/src/Integrands/KirchhoffLovePlate.h @@ -15,6 +15,8 @@ #define _KIRCHHOFF_LOVE_PLATE_H #include "IntegrandBase.h" +#include "Vec3.h" +#include 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& 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 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 diff --git a/src/SIM/AdaptiveSIM.C b/src/SIM/AdaptiveSIM.C index 9a5ea9f6..e7390d40 100644 --- a/src/SIM/AdaptiveSIM.C +++ b/src/SIM/AdaptiveSIM.C @@ -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); +} diff --git a/src/SIM/AdaptiveSIM.h b/src/SIM/AdaptiveSIM.h index 06acd790..17e944e8 100644 --- a/src/SIM/AdaptiveSIM.h +++ b/src/SIM/AdaptiveSIM.h @@ -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: diff --git a/src/SIM/SIM1D.C b/src/SIM/SIM1D.C index 395d7633..a4342e3a 100644 --- a/src/SIM/SIM1D.C +++ b/src/SIM/SIM1D.C @@ -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(utl::parseRealFunc(cline,d)); - else - { - std::cout << d; - myScalars[code] = new ConstFunc(d); - } + + cline = strtok(NULL," "); + myScalars[code] = const_cast(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(utl::parseRealFunc(cline,pd)); - else - myScalars[code] = new ConstFunc(pd); + cline = strtok(NULL," "); + myScalars[code] = const_cast(utl::parseRealFunc(cline,pd)); } std::cout << std::endl; } diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C index 2a77bf63..267abe93 100644 --- a/src/SIM/SIM3D.C +++ b/src/SIM/SIM3D.C @@ -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(utl::parseRealFunc(cline,d)); - else - { - std::cout << d; - myScalars[code] = new ConstFunc(d); - } + + cline = strtok(NULL," "); + myScalars[code] = const_cast(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(utl::parseRealFunc(cline,pd)); - else - myScalars[code] = new ConstFunc(pd); + cline = strtok(NULL," "); + myScalars[code] = const_cast(utl::parseRealFunc(cline,pd)); } if (pface < 10) std::cout << std::endl; } diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 6cf00a51..cd232494 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -950,35 +950,37 @@ bool SIMbase::systemModes (std::vector& 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; } diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index df41d6e0..500af0d3 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -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 diff --git a/src/Utility/Functions.C b/src/Utility/Functions.C index d4cc326e..5df9d50b 100644 --- a/src/Utility/Functions.C +++ b/src/Utility/Functions.C @@ -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<<","< 0 && (cline = strtok(NULL," "))) { diff --git a/src/Utility/VTF.C b/src/Utility/VTF.C index b280f9a1..607d4cc3 100644 --- a/src/Utility/VTF.C +++ b/src/Utility/VTF.C @@ -23,7 +23,7 @@ #endif #include "ElementBlock.h" #include -#include +#include 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 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& 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::const_iterator cit; @@ -372,7 +379,7 @@ bool VTF::writeVectors (const std::map& 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& 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& 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& 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& 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); diff --git a/src/Utility/VTF.h b/src/Utility/VTF.h index a23078ea..37a0eb75 100644 --- a/src/Utility/VTF.h +++ b/src/Utility/VTF.h @@ -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& 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 myDBlock; //!< Displacement blocks std::vector myVBlock; //!< Vector field blocks std::vector mySBlock; //!< Scalar field blocks - std::vector< std::pair > myBlocks; //!< The FE geometry - VTFAGeometryBlock* geoBlock; + std::vector< std::pair > myBlocks; //!< FE geometry + int pointGeoID; //!< ID of current point vector geometry block }; #endif