1 /* 2 XLiFE++ is an extended library of finite elements written in C++ 3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4 5 This program is free software: you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation, either version 3 of the License, or 8 (at your option) any later version. 9 This program is distributed in the hope that it will be useful, 10 but WITHOUT ANY WARRANTY; without even the implied warranty of 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 GNU General Public License for more details. 13 You should have received a copy of the GNU General Public License 14 along with this program. If not, see <http://www.gnu.org/licenses/>. 15 */ 16 17 /*! 18 \file RefElement.hpp 19 \authors D. Martin, N. Kielbasiewicz 20 \since 27 nov 2002 21 \date 6 aug 2012 22 23 \brief Definition of the xlifepp::RefElement class 24 25 A xlifepp::RefElement object represents the standard FE reference element. 26 It is mainly defined by a xlifepp::GeomRefElement and an interpolation which 27 implies a predefined set of degrees of freedom (D.o.F) 28 Only required xlifepp::RefElement objects are instantiated at run time. 29 30 Inheriting classes : 31 - xlifepp::RefSegment (1D) 32 - xlifepp::RefTriangle, xlifepp::RefQuadrangle (2D) 33 - xlifepp::RefTetrahedron, xlifepp::RefPrism, xlifepp::RefHexahedron (3D) 34 35 IMPORTANT : A finite element involves a lot of arrays of numbering such as 36 local numbers in an element, on a side ... 37 38 These numbers always start with first item bearing number 1, though this value is 39 stored in item 0 of an array or vector. 40 41 IMPORTANT : Implicit local numbering of D.o.F 42 43 D.o.F in an element are stored and referred to according to what is called 44 "the local numbering of D.o.F's" which never appears explicitly in the code. 45 Those local numbers of D.o.F can only be "seen" in the arrays of local numbering 46 on edges or faces "numbering_of_DOF". 47 Sketches of local numbering for 1d an 2d elements can be found in implementation 48 of classes xlifepp::RefSegment, xlifepp::RefTriangle or xlifepp::RefQuadrangle. 49 50 That "IMPLICIT" numbering of D.o.F's should obey the following rules which 51 are used to create global numbering of D.o.F. in a Finite Element interpolation space : 52 53 -# D.o.F's supported by vertices of element should be numbered first in the order 54 of (local) numbering of vertices, that is D.o.F supported by first vertex 55 in element should bear first (or smallest) numbers, and so on. 56 This set of D.o.F's can be empty, in which case nbDofsOnVertices_ = 0, 57 otherwise nbDofsOnVertices_ is the number of such D.o.F's. 58 -# each first D.o.F's supported by a node on an edge (excluding vertices) or 59 supported by an edge should be numbered ascendingly according to local numbering 60 of edges, then continue with second D.o.F's and so on 61 that is first D.o.F with support in edge #1 should bear first number following 62 the local numbers of D.o.F's on vertices, then come first D.o.F's on edge #2 and so on. 63 This set of D.o.F's can be empty, in which case nbDofsInSideOfSides_ = 0, 64 otherwise nbDofsInSideOfSides_ is the number of such D.o.F.s 65 -# idem for D.o.F's on faces for 3d elements (replace nbDofsInSideOfSides_ with nbDofsInSides_) 66 -# then come all other D.o.F's (for instance, internal or non-sharable D.o.F's) 67 68 As an example, nodes of Pk-Lagrange triangles are numbered as follows 69 70 \verbatim 71 2 2 2 2 2 72 | \ | \ | \ | \ | \ 73 3---1 5 4 5 7 5 10 5 13 74 k=1 | \ | \ | \ | \ 75 3---6---1 8 10 4 8 14 7 8 17 10 76 k=2 | \ | \ | \ 77 3---6---9---1 11 15 13 4 11 20 19 7 78 k=3 | \ | \ 79 3---6---9--12---1 14 18 21 16 4 80 k=4 | \ 81 3---6---9--12--15--1 82 k=5 83 \endverbatim 84 */ 85 #ifndef REF_ELEMENT_HPP 86 #define REF_ELEMENT_HPP 87 88 #include "config.h" 89 #include "utils.h" 90 #include "Interpolation.hpp" 91 #include "ShapeValues.hpp" 92 #include "RefDof.hpp" 93 #include "GeomRefElement.hpp" 94 95 96 namespace xlifepp 97 { 98 class Quadrature; //forward class 99 100 /*! types of map applied to shape functions of reference element 101 standard map : J{-t} Lagrange FE 102 contravariantPiolaMap : J*J{-t}/|J| div conforming FE 103 covariantPiolaMap : J{-t}*J{-t} curl conforming FE 104 */ 105 enum FEMapType {_standardMap,_contravariantPiolaMap,_covariantPiolaMap}; 106 107 /*! 108 types of compatibility rule to applied to side dofs 109 _noDofCompatibility : no rule to apply (case of Lagrange dofs) 110 _signDofCompatibility : sign correction due to normal/tangent vector orientation (case of Hcurl, Hdiv element) 111 */ 112 enum DofCompatibility {_noDofCompatibility,_signDofCompatibility}; 113 114 /*! 115 \class RefElement 116 A RefElement object represents the standard FE reference element. 117 */ 118 typedef std::vector<std::pair<ShapeType, std::vector<number_t> > > splitvec_t; 119 class RefElement 120 { 121 public: 122 // pointers to GeomRefElement and Interpolation 123 GeomRefElement* geomRefElem_p; //!< geometric data is gathered from geometric reference element 124 const Interpolation* interpolation_p; //!< interpolation parameters 125 std::vector<RefDof*> refDofs; //!< vector of local reference Degrees of Freedom 126 FEMapType mapType; //!< type of map applied to shape functions of reference element 127 DofCompatibility dofCompatibility; //!< compatibility rule to applied to side dofs (_noDofCompatibility,_signDofCompatibility) 128 dimen_t dimShapeFunction; //!< dimension of shape function 129 bool rotateDof; //!< if true, apply rotation (given by childs) to shape values in order to insure correct matching on shared side 130 number_t maxDegree; //!< maximum degree of shape functions 131 132 protected: 133 string_t name_; //!< reference element name (for documentation or debug purposes) 134 number_t nbDofs_; //!< number of Degrees Of Freedom 135 number_t nbPts_; //!< number of Points (for local coordinates of points) (NEEDED only as dimension of local coordinates array ?) 136 // number_t nbNodes_; //!< number of Nodes (i.e. D.o.F support, might not be useful) 137 138 139 // These are only needed to create global numbering of D.o.F's 140 number_t nbDofsOnVertices_; //!< sharable D.o.F on vertices are first D.o.F in element 141 number_t nbDofsInSideOfSides_; //!< sharable D.o.F on side of sides (edges), excluding the previous ones, are next 142 number_t nbDofsInSides_; //!< sharable D.o.F on sides (faces or edges), excluding the previous ones, are next 143 number_t nbInternalDofs_; //!< finally come non-sharable D.o.F's 144 145 /* 146 These are only used (and should only be relied upon) for FE computation on sides 147 They need not be defined for all kinds of elements, for instance for elements 148 each side of which does not "naturally" define a finite element in lower dimension 149 (e.g. Hcurl or Hdiv edge elements; also see actualRefElement() function) 150 */ 151 std::vector<RefElement*> sideRefElems_; //!< pointers to side reference elements (in 2d sideRefElems_ is set to sideOfSideRefElems_) 152 std::vector<RefElement*> sideOfSideRefElems_; //!< pointers to side of side (1D) reference elements 153 154 // vector<ReferenceSideInterpolation*> side_interpolations_; 155 // interpolation data on sides (in 2d side_interpolations_ is set to edge_interpolations_) 156 // vector<ReferenceSideInterpolation*> edge_interpolations_; 157 // interpolation data on edges 158 // Matrix<real_t> shapeFunCoeffs; //!< matrix representing shape functions coefficients in polynoms basis 159 160 public: 161 std::vector<std::vector<number_t> > sideDofNumbers_; //!< list of dof numbers for each side 162 std::vector<std::vector<number_t> > sideOfSideDofNumbers_; //!< list of dof numbers for each side of side (edge in 3D) 163 //internal side dofs mapping when vertices of side are permuted (used when build global dofs) dofsMap(const number_t & i,const number_t & j,const number_t & k=0) const164 virtual std::vector<number_t> dofsMap(const number_t& i, const number_t& j, const number_t& k=0) const //!< dofs permutation 165 {return std::vector<number_t>();} //no permutation sideDofsMap(const number_t & n,const number_t & i,const number_t & j,const number_t & k=0) const166 virtual number_t sideDofsMap(const number_t & n, const number_t & i, const number_t & j, const number_t & k=0) const //face in 3d or edge in 2d 167 {return n;} sideofsideDofsMap(const number_t & n,const number_t & i,const number_t & j=0) const168 virtual number_t sideofsideDofsMap(const number_t & n, const number_t & i, const number_t & j=0) const //edge in 3d 169 {return n;} 170 number_t sideOf(number_t) const; //return side number of a dof if it is located on a side (0 if not), 2D and 3D 171 number_t sideOfSideOf(number_t) const; //return side of side number of a dof if it is located on a side of side (0 if not), 3D only 172 173 //data and tools to manage formal representation of shape values, say shape values as polynomials (optional) 174 public: 175 PolynomialsBasis Wk; //!< shape functions basis 176 std::vector<PolynomialsBasis> dWk; //!< derivatives of shape functions basis 177 178 public: 179 bool hasShapeValues; //!< tell that ref element has shapevalues, generally true 180 mutable ShapeValues shapeValues; //!< to store shape values at a point 181 mutable std::vector<ShapeValues> shvs_; //!< temporary structure to store shape values at some points 182 mutable std::map<Quadrature*,std::vector<ShapeValues> > qshvs_; //!< temporary structure to store shape values at some points of quadrature 183 mutable std::map<Quadrature*,std::vector<ShapeValues> > qshvs_aux; //!< temporary structure to store auxiliary shape values at some points of quadrature 184 185 private: 186 static const bool isSharable_ = true; //!< returns true is sharable (means that it can be a side of another element) 187 188 public: 189 static std::vector<RefElement*> theRefElements; //!< vector to store run-time RefElement pointers 190 static void clearGlobalVector(); //!< delete all RefElement objects 191 192 //-------------------------------------------------------------------------------- 193 // Constructors, Destructor 194 // Warning* Constructor should not be invoked directly: 195 // use RefElementFind function to create or find such an object 196 //-------------------------------------------------------------------------------- 197 public: 198 RefElement(); //!< default constructor 199 RefElement(ShapeType, const Interpolation* ); //!< constructor by shape & interpolation 200 virtual ~RefElement(); //!< destructor 201 202 private: // copy of a RefElement object is forbidden 203 RefElement(const RefElement&); //!< no copy constructor 204 void operator=(const RefElement&); //!< no assignment operator= 205 206 public: 207 //-------------------------------------------------------------------------------- 208 // Public accessors 209 //-------------------------------------------------------------------------------- 210 211 //! returns associated Geometric Reference Element interpolation() const212 const Interpolation& interpolation() const { return *interpolation_p; } 213 //! returns element numtype (order for Lagrange element) order() const214 number_t order() const {return interpolation_p->numtype;} isSimplex() const215 bool isSimplex() const //! returns true for a segment or triangle or tetrahedron 216 {return geomRefElem_p->isSimplex();} 217 //! returns element dimension dim() const218 dimen_t dim() const { return geomRefElem_p->dim(); } 219 //! returns (protected) reference element name name() const220 string_t name() const { return name_; } 221 //! returns (protected) number of points nbPts() const222 number_t nbPts() const { return nbPts_; } 223 //! returns (protected) number of D.o.F supported by vertices nbDofsOnVertices() const224 number_t nbDofsOnVertices() const { return nbDofsOnVertices_; } 225 //! returns (protected) number of D.o.F supported by edges excluding D.o.F support by vertices nbDofsInSideOfSides() const226 number_t nbDofsInSideOfSides() const { return nbDofsInSideOfSides_; } 227 //! returns (protected) number of D.o.F supported by faces excluding D.o.F supported by vertices or edges nbDofsInSides() const228 number_t nbDofsInSides() const { return nbDofsInSides_; } 229 //! returns this if sideNum = 0 or reference element of sideDim side sideNum(sideNum = 1, ...) 230 // Eric : I am not sure that it is longer true actualRefElement(const number_t sideNum=0,const dimen_t sideDim=0) const231 const RefElement* actualRefElement(const number_t sideNum = 0, const dimen_t sideDim = 0) const 232 { 233 if ( sideNum == 0 ) { return this; } 234 else if ( sideDim == 1 ) { return sideOfSideRefElems_[sideNum - 1]; } 235 else { return sideRefElems_[sideNum - 1];} 236 } 237 sideRefElems() const238 const std::vector<RefElement*>& sideRefElems() const 239 {return sideRefElems_;} 240 /*! 241 returns this if sideNum = 0 or reference element of side sideNum (sideNum = 1, ...) 242 non const version, to avoid a forbidden const to non const cast in Space class 243 */ refElement(number_t sideNum=0)244 RefElement* refElement(number_t sideNum = 0) 245 { 246 if ( sideNum == 0 ) { return this; } 247 return sideRefElems_[sideNum - 1]; 248 } 249 //! returns this if sideNum = 0 or reference element of side sideNum (sideNum = 1, ...) refElement(number_t sideNum=0) const250 const RefElement* refElement(number_t sideNum = 0) const //! 251 { 252 if ( sideNum == 0 ) { return this; } 253 return sideRefElems_[sideNum - 1]; 254 } 255 //! return associated Geometric Reference Element (i=0) or GeomRefElement of its side i geomRefElement(number_t sideNum=0) const256 const GeomRefElement* geomRefElement(number_t sideNum = 0) const 257 { 258 if ( sideNum == 0 ) { return geomRefElem_p; } 259 return sideRefElems_[sideNum - 1]->geomRefElem_p; 260 } 261 // ReferenceSideInterpolation* sideInterpolation(const number_t sideNum, const dimen_t sideDim = 0) const 262 // //! return interpolation pointer if sideNum = 0 263 // //! interpolation pointer of sideDim side sideNum (sideNum = 1, ...) 264 // { 265 // if ( sideDim == 1 ) return edge_interpolations_[sideNum-1]; 266 // else return side_interpolations_[sideNum-1]; 267 // } 268 //! returns number of element D.o.F's (sideNum = 0) number of D.o.F's on a given sideDim side sideNum (sideNum = 1, ...) 269 number_t nbDofs(const number_t sideNum = 0, const dimen_t sideDim = 0) const; 270 //! returns number of element internal D.o.F (sideNum = 0 ) or number of internal D.o.F's on sideDim side sideNum (sideNum = 1, ...) 271 number_t nbInternalDofs(const number_t sideNum = 0, const dimen_t sideDim = 0) const; 272 273 ShapeType shapeType() const; //!<returns shape of element 274 275 // const std::vector<number_t>& vertexDofNumbers()const //! dof numbers of dofs on vertices 276 // {return vertexDofNumbers_;} 277 278 //-------------------------------------------------------------------------------- 279 // Public Member Utility functions (declarations) 280 //-------------------------------------------------------------------------------- 281 void buildPolynomialTree(); //!< build tree representation of polynomial shape functions computeShapeFunctions()282 virtual void computeShapeFunctions() //! compute shape functions 283 { error("not_handled","computeShapeFunctions()"); } 284 285 size_t shapeValueSize() const; //!< length of shape function storage computed from set of D.o.Fs 286 // computes shape values 287 virtual void computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv = true) const = 0; computeShapeValues(std::vector<real_t>::const_iterator it_pt,const bool withDeriv=true) const288 void computeShapeValues(std::vector<real_t>::const_iterator it_pt, const bool withDeriv = true) const 289 {computeShapeValues(it_pt,shapeValues,withDeriv);} 290 std::vector<ShapeValues> getShapeValues(const std::vector< std::vector<real_t> >& xs, const bool withDeriv = true) const; 291 void computeShapeValuesFromShapeFunctions(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const; computeShapeValuesFromShapeFunctions(std::vector<real_t>::const_iterator it_pt,const bool withDeriv) const292 void computeShapeValuesFromShapeFunctions(std::vector<real_t>::const_iterator it_pt, const bool withDeriv) const 293 {computeShapeValuesFromShapeFunctions(it_pt,shapeValues,withDeriv);} rotateDofs(const std::vector<number_t> &,ShapeValues & shv,bool withDerivative=false) const294 virtual void rotateDofs(const std::vector<number_t>&, ShapeValues& shv, bool withDerivative=false) const {} //!< apply rotation to some shapevalues (rotation given by children) 295 296 //split stuff 297 virtual std::vector<std::vector<number_t> > splitP1() const; //!< return nodes numbers of first order elements (P1) when splitting current element 298 //! return nodes numbers of first order elements of same shape when splitting current element 299 // (Lagrange elements only) 300 virtual splitvec_t splitO1() const; 301 //! returns reference to splitting scheme computed by splitO1 302 virtual const splitvec_t& getO1splitting() const = 0; 303 virtual std::vector<std::pair<number_t,number_t> > splitP1Side(number_t s) const; //!< return side numbers of first order elements (P1) when splitting current element 304 305 //print utility 306 virtual void print(std::ostream&, bool withDerivative=false) const; //!< print RefElement print(PrintStream & os,bool withDerivative=false) const307 virtual void print(PrintStream& os, bool withDerivative=false) const //! print RefElement 308 {print(os.currentStream(), withDerivative);} print(CoutStream & os,bool withDerivative=false) const309 virtual void print(CoutStream& os, bool withDerivative=false) const //! print RefElement 310 { print(std::cout,withDerivative); 311 if(os.traceOnFile) print(os.printStream->currentStream(), withDerivative); 312 } 313 friend std::ostream& operator<<(std::ostream&, const RefElement&); 314 void printShapeValues(std::ostream&, std::vector<real_t>::const_iterator it_pt) const; //!< print RefElement shape functions printShapeValues(PrintStream & os,std::vector<real_t>::const_iterator it_pt) const315 void printShapeValues(PrintStream& os, std::vector<real_t>::const_iterator it_pt) const //! print RefElement shape functions 316 {printShapeValues(os.currentStream(), it_pt);} 317 void printShapeValues(std::vector<real_t>::const_iterator it_pt) const; //!< print RefElement shape functions 318 //! print reference element as P1 reference 319 virtual void outputAsP1(std::ofstream& os, const int, const std::vector<number_t>& ranks) const = 0; 320 321 protected: 322 void sideOfSideRefElement(); //!< reference element constructors on element edges 323 void sideRefElement(); //!< reference element constructors on element faces 324 void lagrangeRefDofs(dimen_t dimE, number_t nbV, number_t nbDI=0, 325 number_t nbE=0, number_t nbDE=0, 326 number_t nbF=0, number_t nbDF=0); //!< creates Reference D.O.F in Lagrange standard reference element 327 void lagrangeRefDofs(dimen_t dimE, number_t nbV, number_t nbDI, 328 number_t nbE, number_t nbDE, 329 std::vector<number_t> nbDF); //!< creates Reference D.O.F in Lagrange standard reference element 330 //-------------------------------------------------------------------------------- 331 // errors handlers 332 //-------------------------------------------------------------------------------- 333 void noSuchFunction(const string_t& s) const; //!< function not yet implemented in a child class 334 335 friend class Element; 336 337 // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // 338 // #ifdef FIG4TEX_OUTPUT 339 // string fig4TexFileName(); 340 // void fig4TeXVertexPt(ofstream&, vector<real_t>::const_iterator, number_t&); 341 // void fig4TeXEdgePt(ofstream&, const number_t, number_t&, const string_t&, const string_t& = ""); 342 // void fig4TeXFacePt(ofstream&, const number_t, const number_t, number_t& pts, const string_t& posit); 343 // void fig4TeXInternalPoints2d(ofstream&, number_t&); 344 // void fig4TeXSmallFace(ofstream&, const number_t, const number_t, const number_t); 345 // #endif /* FIG4TEX_OUTPUT */ 346 // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // 347 348 public: 349 // splitting tools 350 std::vector<std::vector<number_t> > splitLagrange1DToP1() const; //!< creates the list of nodes numbers of first order segment elements when splitting Lagrange element 351 std::vector<std::vector<number_t> > splitLagrange2DToP1() const; //!< creates the list of nodes numbers of first order triangle elements when splitting Lagrange element 352 std::vector<std::vector<number_t> > splitLagrange3DToP1() const; //!< creates the list of nodes numbers of first order tetrahedron elements when splitting Lagrange element 353 // splitvec_t splitLagrange2DToQ1() const; //!< creates the list of nodes numbers of first order quadrangle elements when splitting Lagrange element (unused) 354 splitvec_t splitLagrange3DToQ1() const; //!< creates the list of nodes numbers of first order hexahedron elements when splitting Lagrange element 355 splitvec_t splitLagrange3DToPrismO1() const; //!< creates the list of nodes numbers of first order prism elements when splitting Lagrange element 356 357 }; // end of class RefElement 358 359 //================================================================================= 360 // Extern class related functions 361 //================================================================================= 362 std::ostream& operator<<(std::ostream&, const RefElement&); //!< output operator 363 void simplexVertexOutput(std::ofstream& os, const int, const int v1, const int v2, const int v3 = 0, const int v4 = 0); 364 bool tetrahedraIntersect(const Point& A1,const Point& B1, const Point& C1, const Point& D1, 365 const Point& A2,const Point& B2, const Point& C2, const Point& D2); //utility 366 367 //================================================================================= 368 // Splitting utility functions 369 //================================================================================= 370 371 /*! 372 \struct SortPointsByXAndY 373 sort nodes according to coordinates 374 */ 375 struct SortPointsByXAndY 376 { 377 bool operator() (const std::pair<number_t,Point>& p, const std::pair<number_t,Point>& q); 378 }; 379 380 int_t whichSide(const Point& p, const Point& n, const std::vector<Point>& pts); //!< determines if a plane splits a set of points 381 bool tetrahedraOverlap(const std::vector<Point>& t1,const std::vector<Point>& t2); //!< determines if 2 tetrahedra overlap 382 bool prismValid(const std::set<std::pair<number_t,Point>, SortPointsByXAndY>& face1, const std::set<std::pair<number_t,Point>,SortPointsByXAndY>& face2); //! determines if a prism is valid 383 bool areNeighbors2D(const Point& p, const Point& q, real_t tol); //!< determines if 2D points are neighbors 384 bool areNeighbors3D(const Point& p, const Point& q, real_t tol); //!< determines if 3D points are neighbors 385 386 /*! 387 RefElementFind 388 definition of a Reference Element by shape number and interpolation 389 use existing Reference Element object if one already exists 390 otherwise create a new one, in both cases returns pointer to the object 391 */ 392 RefElement* findRefElement(ShapeType, const Interpolation*); 393 RefElement* selectRefSegment(const Interpolation*); 394 RefElement* selectRefTriangle(const Interpolation*); 395 RefElement* selectRefQuadrangle(const Interpolation*); 396 RefElement* selectRefTetrahedron(const Interpolation*); 397 RefElement* selectRefPrism(const Interpolation*); 398 RefElement* selectRefHexahedron(const Interpolation*); 399 RefElement* selectRefPyramid(const Interpolation*); 400 401 } // end of namespace xlifepp 402 403 #endif /* REF_ELEMENT_HPP */ 404