00001
00029 #ifndef ELEMENTGEOMETRY
00030 #define ELEMENTGEOMETRY
00031
00032 #include "AuxDefs.h"
00033 #include <string>
00034 #include <vector>
00035 #include <iostream>
00036
00037 #include <cmath>
00056 class ElementGeometry
00057 {
00058 public:
00059 inline ElementGeometry(){}
00060
00061 inline virtual ~ElementGeometry(){}
00062
00063 inline ElementGeometry(const ElementGeometry &){}
00064
00065 virtual ElementGeometry * clone() const = 0;
00066
00068 virtual const size_t getNumVertices() const = 0;
00069
00071 virtual const std::vector <GlobalNodalIndex> &getConnectivity() const = 0;
00072
00074 virtual const std::string getPolytopeName() const = 0;
00075
00077 virtual size_t getSpatialDimension () const = 0;
00078
00080 virtual const size_t getParametricDimension() const = 0;
00081
00083 virtual const size_t getEmbeddingDimension() const = 0;
00084
00088 virtual void map(const double * X, double *Y) const = 0;
00089
00097 virtual void dMap(const double * X, double *DY, double &Jac) const = 0;
00098
00103 virtual bool consistencyTest(const double * X, const double Pert) const = 0;
00104
00106 virtual size_t getNumFaces() const = 0;
00107
00115 virtual ElementGeometry * getFaceGeometry(size_t e) const = 0;
00116
00120 virtual const double getInRadius() const = 0;
00121
00125 virtual const double getOutRadius() const = 0;
00126
00131 virtual void computeNormal (size_t e, std::vector<double>& vNormal) const = 0;
00132 };
00133
00137 template <size_t SPD>
00138 class AbstractGeom : public ElementGeometry {
00139 private:
00140 const std::vector<double>& globalCoordVec;
00141 std::vector<GlobalNodalIndex> connectivity;
00142
00143 public:
00150 AbstractGeom (const std::vector<double>& globalCoordVec, const std::vector<GlobalNodalIndex>& connectivity)
00151 :ElementGeometry (), globalCoordVec(globalCoordVec), connectivity (connectivity) {
00152
00153 }
00154
00155 AbstractGeom (const AbstractGeom<SPD>& that)
00156 : ElementGeometry (that), globalCoordVec (that.globalCoordVec), connectivity (that.connectivity) {
00157
00158 }
00159
00160 virtual size_t getSpatialDimension () const {
00161 return SPD;
00162 }
00163
00164 virtual const std::vector<GlobalNodalIndex>& getConnectivity () const {
00165 return connectivity;
00166 }
00167
00168 virtual bool consistencyTest (const double* X, const double Pert) const {
00169 double *DYNum = new double[getParametricDimension() * getEmbeddingDimension()];
00170 double *DY = new double[getParametricDimension() * getEmbeddingDimension()];
00171 double *Xpert = new double[getParametricDimension()];
00172 double *Yplus = new double[getEmbeddingDimension()];
00173 double *Yminus = new double[getEmbeddingDimension()];
00174 double Jac;
00175
00176 if (Pert <= 0)
00177 std::cerr << "ElementGeometry::ConsistencyTest - Pert cannot be less or equal than zero\n";
00178
00179 for (size_t a = 0; a < getParametricDimension(); a++)
00180 Xpert[a] = X[a];
00181
00182 dMap(X, DY, Jac);
00183
00184 for (size_t a = 0; a < getParametricDimension(); a++) {
00185 Xpert[a] = X[a] + Pert;
00186 map(Xpert, Yplus);
00187
00188 Xpert[a] = X[a] - Pert;
00189 map(Xpert, Yminus);
00190
00191 Xpert[a] = X[a];
00192
00193 for (size_t i = 0; i < getEmbeddingDimension(); i++)
00194 DYNum[a * getEmbeddingDimension() + i] = (Yplus[i] - Yminus[i]) / (2 * Pert);
00195 }
00196
00197 double error = 0;
00198 double normX = 0;
00199 double normDYNum = 0;
00200 double normDY = 0;
00201
00202 for (size_t a = 0; a < getParametricDimension(); a++) {
00203 normX += X[a] * X[a];
00204
00205 for (size_t i = 0; i < getEmbeddingDimension(); i++) {
00206 error += pow(DY[a * getEmbeddingDimension() + i] - DYNum[a * getEmbeddingDimension() + i], 2.);
00207 normDY += pow(DY[a * getEmbeddingDimension() + i], 2.);
00208 normDYNum += pow(DYNum[a * getEmbeddingDimension() + i], 2.);
00209 }
00210 }
00211 error = sqrt(error);
00212 normX = sqrt(normX);
00213 normDY = sqrt(normDY);
00214 normDYNum = sqrt(normDYNum);
00215
00216 delete[] Yplus;
00217 delete[] Yminus;
00218 delete[] Xpert;
00219 delete[] DYNum;
00220 delete[] DY;
00221
00222 if (error * (normX + Pert) < (normDY < normDYNum ? normDYNum : normDY) * Pert * 10)
00223 return true;
00224 else
00225 return false;
00226 }
00227
00228
00229
00230 protected:
00231
00237 double getCoordinate (size_t a, size_t i) const {
00238
00239 size_t index = getConnectivity ()[a] * getSpatialDimension() + i;
00240 return globalCoordVec[index];
00241 }
00242
00246 const std::vector<double>& getGlobalCoordVec () const { return globalCoordVec; }
00247
00248
00249 };
00250
00251 #endif
