NGSolve
4.9
|
00001 #ifndef FILE_ELEMENTTRANSFORMATION 00002 #define FILE_ELEMENTTRANSFORMATION 00003 00004 /*********************************************************************/ 00005 /* File: elementtransformation.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 25. Mar. 2000 */ 00008 /*********************************************************************/ 00009 00010 00011 00012 /* 00013 Transformation from reference element to actual element 00014 */ 00015 00016 00017 00018 namespace ngfem 00019 { 00020 00024 class NGS_DLL_HEADER ElementTransformation 00025 { 00026 protected: 00028 int elnr; 00030 int elindex; 00032 // int dim; 00034 // bool boundary; 00036 ELEMENT_TYPE eltype; 00037 00038 bool higher_integration_order; 00039 00041 bool iscurved; 00042 00043 public: 00045 // ElementTransformation * specific; 00047 ElementTransformation () { higher_integration_order = false; } // specific = NULL; } 00049 virtual ~ElementTransformation() { ; } // delete specific; } 00051 virtual void SetElement (bool aboundary, int aelnr, int aelindex) 00052 { 00053 // boundary = Boundary(); // aboundary; 00054 elnr = aelnr; 00055 elindex = aelindex; 00056 // if (specific) specific -> SetElement (aboundary, aelnr, aelindex); 00057 // iscurved = false; 00058 } 00060 void SetElementType (ELEMENT_TYPE aet) { eltype = aet; } 00062 int GetElementNr () const { return elnr; } 00064 int GetElementIndex () const { return elindex; } 00066 ELEMENT_TYPE GetElementType () const { return eltype; } 00068 virtual void CalcJacobian (const IntegrationPoint & ip, 00069 FlatMatrix<> dxdxi) const = 0; 00070 00072 virtual void CalcPoint (const IntegrationPoint & ip, 00073 FlatVector<> point) const = 0; 00074 00076 virtual void CalcPointJacobian (const IntegrationPoint & ip, 00077 FlatVector<> point, FlatMatrix<> dxdxi) const = 0; 00078 00080 virtual void CalcMultiPointJacobian (const IntegrationRule & ir, 00081 BaseMappedIntegrationRule & mir) const = 0; 00082 00084 void CalcNormalVector (const IntegrationPoint & ip, 00085 FlatVector<> nv, 00086 LocalHeap & lh) const 00087 { 00088 if (Boundary()) 00089 { 00090 if (SpaceDim() == 2) 00091 { 00092 Mat<2,1> dxdxi; 00093 CalcJacobian (ip, dxdxi); 00094 // Ng_GetSurfaceElementTransformation (elnr+1, &ip(0), 0, &dxdxi(0)); 00095 double len = sqrt (sqr (dxdxi(0,0)) + sqr(dxdxi(1,0))); 00096 nv(0) = -dxdxi(1,0) / len; //SZ 00097 nv(1) = dxdxi(0,0) / len; 00098 } 00099 else 00100 { 00101 Mat<3,2> dxdxi; 00102 CalcJacobian (ip, dxdxi); 00103 // Ng_GetSurfaceElementTransformation (elnr+1, &ip(0), 0, &dxdxi(0)); 00104 nv(0) = dxdxi(1,0) * dxdxi(2,1) - dxdxi(2,0) * dxdxi(1,1); 00105 nv(1) = dxdxi(2,0) * dxdxi(0,1) - dxdxi(0,0) * dxdxi(2,1); 00106 nv(2) = dxdxi(0,0) * dxdxi(1,1) - dxdxi(1,0) * dxdxi(0,1); 00107 nv /= L2Norm (nv); 00108 } 00109 } 00110 } 00111 00112 00114 virtual int SpaceDim () const = 0; 00115 00117 virtual bool Boundary() const = 0; 00118 00119 void SetHigherIntegrationOrder(void) {higher_integration_order = true;} 00120 void UnSetHigherIntegrationOrder(void) {higher_integration_order = false;} 00121 bool HigherIntegrationOrderSet(void) const 00122 { 00123 return higher_integration_order; 00124 } 00125 00127 virtual bool IsCurvedElement() const 00128 { 00129 return false; 00130 } 00131 00132 virtual void GetSort (FlatArray<int> sort) const 00133 { ; } 00134 00136 virtual BaseMappedIntegrationPoint & operator() (const IntegrationPoint & ip, LocalHeap & lh) const = 0; 00137 00139 virtual BaseMappedIntegrationRule & operator() (const IntegrationRule & ir, LocalHeap & lh) const = 0; 00140 00141 private: 00142 ElementTransformation (const ElementTransformation & eltrans2) { ; } 00143 ElementTransformation & operator= (const ElementTransformation & eltrans2) 00144 { return *this; } 00145 }; 00146 00147 00148 00149 00150 00151 00152 00153 00154 00155 00156 00157 /* 00158 Transformation from reference element to physical element. 00159 Uses finite element fel to describe mapping 00160 */ 00161 template <int DIMS, int DIMR> 00162 class NGS_DLL_HEADER FE_ElementTransformation : public ElementTransformation 00163 { 00165 const ScalarFiniteElement<DIMS> * fel; 00166 00168 Matrix<> pointmat; 00170 // bool pointmat_ownmem; 00171 00173 FlatMatrix<> nvmat; 00174 public: 00176 FE_ElementTransformation (); 00177 00179 ~FE_ElementTransformation (); 00180 00182 virtual void SetElement (const FiniteElement * afel, int aelnr, int aelindex) 00183 { 00184 fel = static_cast<const ScalarFiniteElement<DIMS>*> (afel); 00185 elnr = aelnr; 00186 elindex = aelindex; 00187 SetElementType (fel->ElementType()); 00188 pointmat.SetSize (DIMR, fel->GetNDof()); 00189 } 00190 00192 const FiniteElement & GetElement () const { return *fel; } 00193 00194 00195 ELEMENT_TYPE GetElementType () const { return fel->ElementType(); } 00196 00198 virtual void CalcJacobian (const IntegrationPoint & ip, 00199 FlatMatrix<> dxdxi) const; 00200 00202 virtual void CalcPoint (const IntegrationPoint & ip, 00203 FlatVector<> point) const; 00204 00206 virtual void CalcPointJacobian (const IntegrationPoint & ip, 00207 FlatVector<> point, 00208 FlatMatrix<> dxdxi) const; 00209 00210 00211 00212 00213 virtual void CalcMultiPointJacobian (const IntegrationRule & ir, 00214 BaseMappedIntegrationRule & bmir) const; 00215 00217 // const FlatMatrix<> & PointMatrix () const { return pointmat; } 00219 FlatMatrix<> PointMatrix () const { return pointmat; } 00221 /* 00222 void AllocPointMatrix (int spacedim, int vertices) 00223 { 00224 if (pointmat_ownmem) delete [] &pointmat(0,0); 00225 pointmat.AssignMemory (spacedim, vertices, new double[spacedim*vertices]); 00226 pointmat_ownmem = 1; 00227 } 00228 */ 00230 const FlatMatrix<> & NVMatrix () const { return nvmat; } 00231 00232 template <typename T> 00233 void CalcNormalVector (const IntegrationPoint & ip, 00234 T & nv, 00235 LocalHeap & lh) const 00236 { 00237 for (int i = 0; i < nvmat.Height(); i++) 00238 nv(i) = nvmat(i,0) ; 00239 } 00240 00242 int SpaceDim () const 00243 { 00244 return pointmat.Height(); 00245 } 00246 00247 bool Boundary(void) const 00248 { 00249 return pointmat.Height() != ElementTopology::GetSpaceDim (fel->ElementType()); 00250 } 00251 00252 void GetSort (FlatArray<int> sort) const 00253 { ; } 00254 00255 00256 virtual BaseMappedIntegrationPoint & operator() (const IntegrationPoint & ip, LocalHeap & lh) const 00257 { 00258 return *new (lh) MappedIntegrationPoint<DIMS,DIMR> (ip, *this); 00259 } 00260 00261 virtual BaseMappedIntegrationRule & operator() (const IntegrationRule & ir, LocalHeap & lh) const 00262 { 00263 return *new (lh) MappedIntegrationRule<DIMS,DIMR> (ir, *this, lh); 00264 } 00265 00266 }; 00267 00268 00269 00270 00271 00272 00273 } 00274 00275 00276 00277 00278 #endif 00279 00280 00281 00282 00283 00284