00001
00024 #ifndef MESHINIT_H_
00025 #define MESHINIT_H_
00026
00027 #include <vector>
00028 #include <string>
00029 #include <algorithm>
00030 #include <exception>
00031
00032 #include <cassert>
00033 #include <cmath>
00034 #include <cstdio>
00035
00036 #include <boost/noncopyable.hpp>
00037
00038 #include "AuxDefs.h"
00039 #include "StandardAVI.h"
00040 #include "Element.h"
00041 #include "Femap.h"
00042 #include "Material.h"
00043 #include "CoordConn.h"
00044 #include "TriLinearCoordConn.h"
00045 #include "TetLinearCoordConn.h"
00046
00047 class MeshInit : boost::noncopyable {
00048
00049 public:
00050 typedef StandardAVI::BCFunc BCFunc;
00051 typedef StandardAVI::BCImposedType BCImposedType;
00052
00053 private:
00054 static const double RHO = 1.0;
00055 static const double MU = 0.5;
00056 static const double LAMBDA = 0.0;
00057 static const int PID = 0;
00058
00059 static const double DELTA = 0.1;
00060 static const double T_INIT = 0.0;
00061
00062
00063 static const size_t MAX_FNAME = 1024;
00064
00065
00066
00067 private:
00068 double simEndTime;
00069 bool wave;
00070
00071
00072
00073
00075 LocalToGlobalMap* l2gMap;
00076 CoordConn* cc;
00077 SimpleMaterial* ile;
00078
00080 std::vector<ElementGeometry*> geomVec;
00081 std::vector<Element*> elemVec;
00082 std::vector<Residue*> massResidueVec;
00083 std::vector<DResidue*> operationsVec;
00084 std::vector<AVI*> aviVec;
00085
00086
00087 std::vector<size_t> fieldsUsed;
00088
00089 double writeInc;
00090 int writeInterval;
00091 std::vector<size_t> aviWriteInterval;
00092 FILE* syncFileWriter;
00093
00094 private:
00095 void stretchInternal (VecDouble& dispOrVel, bool isVel) const ;
00096 void getBCs (const Element& e, std::vector< std::vector<BCImposedType> >& itypeMat,
00097 std::vector<BCFunc>& bcfuncVec) const;
00098
00099 static bool computeDiffAVI (std::vector<AVI*> listA, std::vector<AVI*> listB, bool printDiff);
00100
00101 template <typename T>
00102 static void destroyVecOfPtr (std::vector<T*>& vec) {
00103 for (typename std::vector<T*>::iterator i = vec.begin (), ei = vec.end ();
00104 i != ei; ++i) {
00105
00106 delete *i;
00107 *i = NULL;
00108 }
00109 }
00110
00111 void destroy () {
00112
00113 destroyVecOfPtr (geomVec);
00114 destroyVecOfPtr (elemVec);
00115 destroyVecOfPtr (massResidueVec);
00116 destroyVecOfPtr (operationsVec);
00117 destroyVecOfPtr (aviVec);
00118
00119 delete l2gMap;
00120 delete cc;
00121 delete ile;
00122
00123 }
00124
00125 public:
00126
00132 MeshInit (double simEndTime, bool wave):
00133 simEndTime (simEndTime), wave(wave) {
00134
00135
00136
00137 if (wave) {
00138 this->writeInc = 0.005;
00139 }
00140 else {
00141 this->writeInc = 0.1;
00142 }
00143
00144 writeInterval = 0;
00145 syncFileWriter = NULL;
00146
00147 }
00148
00149 virtual ~MeshInit () {
00150 destroy ();
00151 }
00152
00161 void initializeMesh (const std::string& fileName, int ndiv);
00162
00163 virtual size_t getSpatialDim () const = 0;
00164
00165 virtual size_t getNodesPerElem () const = 0;
00166
00167 bool isWave () const { return wave; }
00168
00169 double getSimEndTime () const { return simEndTime; }
00170
00172 int getNumElements () const { return cc->getNumElements (); }
00173
00175 int getNumNodes () const { return cc->getNumNodes (); }
00176
00178 unsigned int getTotalNumDof () const { return l2gMap->getTotalNumDof (); }
00179
00180 const std::vector<AVI*>& getAVIVec () const { return aviVec; }
00181
00184
00185 const LocalToGlobalMap& getLocalToGlobalMap () const { return *l2gMap; }
00186
00193 void setupDisplacements (VecDouble& disp) const { stretchInternal (disp, false); }
00194
00201 void setupVelocities (VecDouble& vel) const { stretchInternal (vel, true); }
00202
00212 void writeSync (const AVI& avi, const VecDouble& Qval, const VecDouble& Vbval, const VecDouble& Tval) ;
00213
00214
00221 bool cmpState (const MeshInit& that) const { return computeDiffAVI (this->aviVec, that.aviVec, false); }
00222
00230 void printDiff (const MeshInit& that) const { computeDiffAVI (this->aviVec, that.aviVec, true); }
00231
00232
00233
00234 protected:
00235
00238 virtual BCFunc getBCFunc (const double* coord) const = 0;
00239
00241 virtual CoordConn* makeCoordConn () const = 0;
00242
00244 virtual const double* getParam () const = 0;
00245
00249 virtual double initDisp (const double* coord, int f) const = 0;
00250
00254 virtual double initVel (const double* coord, int f) const = 0;
00255
00257 virtual size_t numFields () const = 0;
00258
00259 };
00260
00261 class TriMeshInit: public MeshInit {
00262 private:
00263
00264 static double topBC (int f, int a, double t) {
00265 if (f == 0) {
00266 return (0.1 * cos (t));
00267 } else {
00268 return (0.0);
00269 }
00270 }
00271
00272 static double botBC (int f, int a, double t) {
00273 return (0.0);
00274 }
00275
00276 public:
00277 static const double PARAM[];
00278
00279 TriMeshInit (double simEndTime, bool wave=false): MeshInit(simEndTime, wave) {
00280 }
00281
00282
00283 virtual size_t getSpatialDim () const {
00284 return TriLinearTraits::SPD;
00285 }
00286 virtual size_t getNodesPerElem () const {
00287 return TriLinearTraits::NODES_PER_ELEM;
00288 }
00289
00290
00291 protected:
00292
00293
00294 virtual CoordConn* makeCoordConn () const { return new TriLinearCoordConn(); }
00295
00296 virtual const double* getParam () const { return PARAM; }
00297
00298 virtual size_t numFields () const { return TriLinearTraits::NFIELDS;}
00299
00300 virtual BCFunc getBCFunc (const double* coord) const {
00301 if (coord[0] == 0.0) {
00302 return botBC;
00303 } else if (coord[0] == 10.0) {
00304 return topBC;
00305 } else {
00306 return NULL;
00307 }
00308 }
00309
00310 virtual double initDisp (const double* coord, int f) const {
00311 double stretch;
00312 if (f == 0) {
00313
00314 stretch = coord[0] * 0.2 - 1.0;
00315 stretch = coord[0] * 0.01 - 0.05;
00316 } else {
00317 stretch = coord[1] * 0.2 - 1.0;
00318 stretch = coord[1] * 0.01 - 0.05;
00319 }
00320
00321 return stretch;
00322 }
00323
00324 virtual double initVel (const double* coord, int f) const {
00325 if (coord[0] == 0.0) {
00326 return 0.1;
00327 } else {
00328 return 0.0;
00329 }
00330 }
00331
00332 };
00333
00334 class TetMeshInit: public MeshInit {
00335 private:
00336
00337 static double topBC (int f, int a, double t) {
00338 if (f == 2) {
00339 return (0.1 * sin (t));
00340 } else {
00341 return (0.0);
00342 }
00343 }
00344
00345 static double botBC (int f, int a, double t) {
00346 return (0.0);
00347 }
00348
00349
00350 public:
00351 static const double PARAM[];
00352
00353 TetMeshInit (double simEndTime, bool wave=false): MeshInit(simEndTime, wave) {
00354 }
00355
00356
00357 virtual size_t getSpatialDim () const {
00358 return TetLinearTraits::SPD;
00359 }
00360 virtual size_t getNodesPerElem () const {
00361 return TetLinearTraits::NODES_PER_ELEM;
00362 }
00363
00364
00365 protected:
00366
00367
00368 virtual CoordConn* makeCoordConn () const { return new TetLinearCoordConn(); }
00369
00370 virtual const double* getParam () const { return PARAM; }
00371
00372 virtual size_t numFields () const { return TetLinearTraits::NFIELDS;}
00373
00374 virtual BCFunc getBCFunc (const double* coord) const {
00375 if (fabs (coord[0]) < 0.01) {
00376 return botBC;
00377 } else if (fabs (coord[0] - 5.0) < 0.01) {
00378 return topBC;
00379 } else {
00380 return NULL;
00381 }
00382 }
00383
00384 virtual double initDisp (const double* coord, int f) const {
00385 double stretch = coord[f] * 0.10 - 0.25;
00386 return stretch;
00387 }
00388
00389 virtual double initVel (const double* coord, int f) const {
00390 return 0.0;
00391 }
00392
00393 };
00394
00395
00396 #endif
