NGSolve  4.9
fem/elementtransformation.hpp
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