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 GeomElement.hpp 19 \authors D. Martin, E. Lunéville, N. Kielbasiewicz 20 \since 13 apr 2012 21 \date 5 jul 2013 22 23 \brief Definition of the xlifepp::GeomElement class 24 25 A xlifepp::GeomElement object represents the geometric support of mesh elements and is of two kinds : 26 - plain mesh element : has always a xlifepp::MeshElement structure and not parent 27 - side of a mesh element : has always a parent and may have a MeshElement structure 28 29 A xlifepp::MeshElement object mainly carries 30 - Reference Element (Lagrange finite element) 31 - some numbering vectors for nodes, sides, sides of side and neighbor elements 32 */ 33 34 35 #ifndef GEOM_ELEMENT_HPP 36 #define GEOM_ELEMENT_HPP 37 38 #include "GeomMapData.hpp" 39 #include "finiteElements.h" 40 #include "utils.h" 41 #include "config.h" 42 43 namespace xlifepp 44 { 45 // forward class declarations 46 class Mesh; 47 class MeshDomain; 48 49 //------------------------------------------------------------------------------------ 50 /*! 51 \class MeshElement 52 object represents the geometric support of a mesh element. 53 It mainly stores global numberings of element 54 */ 55 //------------------------------------------------------------------------------------ 56 class MeshElement 57 { 58 public: 59 std::vector<Point*> nodes; //!< a subset of Mesh::nodes (direct access to nodes) 60 std::vector<number_t> nodeNumbers; //!< node numbers : a subset of Mesh::nodes (index starts from 1) 61 std::vector<number_t> vertexNumbers; //!< vertexNumbers : a subset of Mesh::nodes (may be equal to nodeNumbers, index starts from 1)) 62 std::vector<real_t> measures; //!< length or area or volume of element and of its sides 63 short int orientation; //!< element orientation (sign of jacobian determinant), 0 means unset orientation 64 Point centroid; //!< centroid of the element, set by computeMeasures 65 real_t size; //!< characteristic size of element, set by computeMeasures 66 mutable GeomMapData* geomMapData_p; //!< useful data for geometric map (jacobian,...) 67 bool linearMap; //!< true if geometric map from ref element is linear (for instance first order simplicial element) 68 69 //possibly additional useful information, only built if required 70 std::vector<number_t> sideNumbers; //!< element side numbers : a subset of Mesh::sides_ (built if required) 71 std::vector<number_t> sideOfSideNumbers; //!< element side of side numbers : : a subset of Mesh::sideOfSides_ (built if required) 72 const MeshElement* firstOrderParent_; //!< parent element when the current element is a first order element coming from an element of order>1 73 74 private: 75 const RefElement* refElt_p; //!< pointer to associated Reference Element (for geometric interpolation) 76 number_t index_; //!< element index (=number_ of GeomElement) 77 const dimen_t spaceDim_; //!< space dimension (node dimension) 78 79 public: 80 //constructor/destructor 81 MeshElement(); //!< default constructor 82 MeshElement(const RefElement*, dimen_t, number_t); //!< constructor of plain element from ReferenceElement and space dimension ~MeshElement()83 ~MeshElement() {if(geomMapData_p!=0) delete geomMapData_p;} //!< destructor 84 void setNodes(std::vector<Point>&); //!< update node pointers vector from Mesh::nodes 85 86 private: 87 void operator=(const MeshElement&); //!< no assignment operator= 88 89 public: 90 //accessors index() const91 number_t index() const //! return element index 92 {return index_;} index()93 number_t& index() //! return element index 94 {return index_;} spaceDim() const95 dimen_t spaceDim() const //! return space dimension 96 {return spaceDim_; } 97 98 //properties order() const99 number_t order() const //! return order of Lagrange interpolation 100 {return refElt_p->order();} isSimplex() const101 bool isSimplex() const //! return true for a segment or triangle or tetrahedron 102 {return refElt_p->isSimplex();} measure(const number_t sideNo=0) const103 real_t measure(const number_t sideNo = 0) const //! return measure of element if side_no = 0, measure of side number side_no > 0 104 {return measures[sideNo];} refElement(number_t i=0) const105 const RefElement* refElement(number_t i = 0) const //! return reference element if i=0 else refelement of side i 106 {return refElt_p->refElement(i);} geomRefElement(number_t i=0) const107 const GeomRefElement* geomRefElement(number_t i = 0) const //! return geometric reference element if i=0 else geomrefelement of side i 108 {return refElement(i)->geomRefElem_p;} elementDim() const109 dimen_t elementDim() const //! return space dimension of element 110 {return geomRefElement()->dim();} shapeType(number_t i=0) const111 ShapeType shapeType(number_t i = 0) const //! return the shape type of element (i=0) and its sides (_segment, _triangle, ...) 112 {return geomRefElement()->shapeType(i);} vertexNumber(number_t i) const113 number_t vertexNumber(number_t i) const //! return the global number of vertex i (i=1,...) 114 {return vertexNumbers[i - 1];} nodeNumber(number_t i) const115 number_t nodeNumber(number_t i) const //! return the global number of node i (i=1,...) 116 {return nodeNumbers[i - 1];} 117 std::vector<number_t> verticesNumbers(number_t s = 0) const; //!< return the global numbers of vertices of element (s=0), of side s>0 118 bool checkLinearMap() const; //!< check if linear map is available 119 vertex(number_t i) const120 Point& vertex(number_t i) const //! return the vertex i (i=1,...) 121 {return *nodes[i-1];} node(number_t i) const122 Point& node(number_t i) const //! return the node i (i=1,...) 123 {return *nodes[i-1];} vertexOnSide(number_t i,number_t s) const124 const Point& vertexOnSide(number_t i, number_t s) const //! return the vertex i (i=1,...) on side s 125 {return *nodes[geomRefElement()->sideVertexNumber(i, s)-1];} vertexOnSideOfSide(number_t i,number_t ss) const126 const Point& vertexOnSideOfSide(number_t i, number_t ss) const //! return the vertex i (i=1,...) on side of side ss 127 {return *nodes[geomRefElement()->sideOfSideVertexNumber(i, ss)-1];} vertexOppositeToSide(number_t s) const128 const Point& vertexOppositeToSide(number_t s) const //! return a vertex opposite to side s 129 {return *nodes[geomRefElement()->vertexOppositeSide(s)-1];} 130 numberOfNodes() const131 number_t numberOfNodes() const //! returns number of nodes 132 {return nodeNumbers.size();} numberOfVertices() const133 number_t numberOfVertices() const //! returns number of vertices 134 {return verticesNumbers().size();} numberOfSides() const135 number_t numberOfSides() const //! returns number of sides 136 {return geomRefElement()->nbSides();} numberOfSideOfSides() const137 number_t numberOfSideOfSides() const //! returns number of sides of sides 138 {return geomRefElement()->nbSideOfSides();} 139 140 // compute measures and orientation 141 void computeMeasures(); //!< compute measure of element and its side 142 void computeMeasure(); //!< compute measure of element 143 void computeMeasureOfSides(); //!< compute measure sides of element 144 void computeOrientation(); //!< compute orientation of element : sign(det(jacobian)) 145 Point center() const; //!< compute the center of element 146 real_t characteristicSize() const; //!< compute a characteristic size of element (diameter of the circumscribed ball or max of length of edges) 147 std::vector<real_t> normalVector(number_t s=0) const;//!< return outward normal vector on side s, if s=0 normal vector to element 148 std::vector<real_t> tangentVector(number_t e) const; //!< return tangent vector on edge 149 150 // various utilities 151 bool contains(const std::vector<real_t>&); //!< return true if mesh element contains a point given as Point or std::vector<real_t> 152 Point projection(const std::vector<real_t>&, real_t&) const; //!< projection of a point onto meshelement 153 std::vector<MeshElement*>splitP1() const; //!< split MeshElement in P1 (simplex and order 1) mesh elements 154 155 // GeomMapData management clearGeomMapData()156 void clearGeomMapData() //! clear geometric map data (i.e geoMapData object) 157 {if(geomMapData_p!=0)delete geomMapData_p; geomMapData_p = 0;} 158 159 void print(std::ostream&) const; //!< print utility print(PrintStream & os) const160 void print(PrintStream& os) const {print(os.currentStream());} 161 friend std::ostream& operator<< (std::ostream&, const MeshElement&); 162 163 }; // end of class MeshElement =================================================== 164 165 std::ostream& operator<< (std::ostream&, const MeshElement&); //!< prints characteristics and point numbers 166 real_t distance(const MeshElement&, const MeshElement&); //!< distance from two MeshElement 167 168 typedef std::pair<MeshElement*, number_t> EltNumPair; //!< useful alias of a pair of MeshElement* and number_t 169 170 //------------------------------------------------------------------------------------ 171 /*! 172 \class GeomElement 173 object represents the geometric support of mesh elements. 174 It either a plain element or a side element (side of a GeomElement). 175 - when it is a plain element, there is always a MeshElement associated to its but no parentSides 176 - when it is a side element, there is at least a parent (parentSides is not empty) and by default no MeshElement, 177 but if required, MeshElement may be constructed (buildSideMeshElement) 178 179 */ 180 //------------------------------------------------------------------------------------ 181 class GeomElement; 182 183 //! useful typedef for GeomElement class 184 typedef std::pair<GeomElement*, number_t> GeoNumPair; //!< alias of a pair of GeomElement* and a number 185 typedef real_t (*ColoringRule)(const GeomElement&, const std::vector<real_t>&); //!< alias of a function describing a GeomElement coloring rule 186 187 class GeomElement 188 { 189 protected: 190 const Mesh* mesh_p; //!< pointer to the const mesh 191 number_t number_; //!< unique id number of geometric element 192 mutable MeshElement* meshElement_p; //!< pointer to a MeshElement (may be 0 in case of side element) 193 std::vector<GeoNumPair> parentSides_; //!< if a side element, stores the pairs of parent element and side number 194 195 public: 196 //additionnal information 197 number_t materialId; //!< material id (default= 0) 198 number_t domainId; //!< domain id (default=0) 199 real_t theta; //!< angle to describe local curvature (theta in x-y plane), default 0 200 real_t phi; //!< angle to describe local curvature (phi in x-z plane), default 0 201 mutable GeomElement* twin_p; //!< twin geomelement in case of geomElement on an interface 202 mutable real_t color; //!< useful to mark an element (default = 0) 203 mutable real_t flag; //!< useful to mark an element in internal process (default = 0) 204 205 //constructor and destructor GeomElement(const Mesh * m=0,number_t n=0)206 GeomElement(const Mesh* m = 0, number_t n = 0) //! basic constructor 207 : mesh_p(m), number_(n), meshElement_p(0), materialId(0), domainId(0), theta(0), phi(0), twin_p(0), color(0), flag(0) {} 208 209 GeomElement(const Mesh*, const RefElement*, 210 dimen_t, number_t); //!< construct a plain element from Mesh, ReferenceElement and space dimension 211 GeomElement(GeomElement*, number_t, number_t); //!< construct a side element from GeomElement parent and side number 212 GeomElement(const GeomElement&); //!< copy constructor 213 GeomElement& operator=(const GeomElement&); //!< assign operator 214 ~GeomElement()215 ~GeomElement() //! destructor 216 {if(meshElement_p != 0) { delete meshElement_p; }} 217 218 //build utilities 219 MeshElement* buildSideMeshElement() const; //!< create a MeshElement object when element it a side element (not automatic) 220 void deleteMeshElement(); //!< destroy meshElement information 221 222 //simple accessors meshP() const223 const Mesh* meshP() const {return mesh_p;} //!< return mesh pointer meshP()224 const Mesh*& meshP() {return mesh_p;} //!< return mesh pointer number() const225 number_t number() const {return number_;} //!< returns element number number()226 number_t& number() {return number_;} //!< returns element number (writable) hasMeshElement() const227 bool hasMeshElement() const //! true if it handles a MeshElement structure 228 {return meshElement_p != 0;} isSideElement() const229 bool isSideElement() const //! true if it is a side element 230 {return parentSides_.size() > 0;} numberOfParents() const231 number_t numberOfParents() const //! number of parents of side element 232 {return parentSides_.size();} 233 234 //complex accessors depending of the GeomElement type (plain or side) 235 const MeshElement* meshElement() const; //!< return meshElement_p if not null else error 236 MeshElement* meshElement() ; //!< return meshElement_p if not null (non const) else error 237 const GeomElement* parent(number_t i = 0) const; //!< return i-th parent (no check) 238 const GeoNumPair parentSide(number_t i) const; //!< return i-th parent and its side (no check) parentSides()239 std::vector<GeoNumPair>& parentSides() //! access to parentSides_ 240 {return parentSides_;} 241 const GeomElement* parentInDomain(const MeshDomain*) const;//!< look for a parent in a domain (if exists return the first one) 242 dimen_t elementDim() const; //!< return dimension of the element 243 dimen_t spaceDim() const; //!< return space dimension, i.e nodes dimension 244 const RefElement* refElement(number_t i = 0) const; //!< return reference element if i=0 else ref element of side i geomRefElement(number_t i=0) const245 const GeomRefElement* geomRefElement(number_t i = 0) const //! return geomrefelement if i=0 else geomrefelement of side i 246 {return refElement(i)->geomRefElement();} 247 248 number_t numberOfNodes() const; //!< returns number of nodes of element 249 number_t numberOfVertices() const; //!< returns number of vertices of element 250 number_t numberOfSides() const; //!< returns number of sides (faces in 3D, edges in 2D) 251 number_t numberOfSideOfSides() const; //!< returns number of edges (edges in 3D) 252 ShapeType shapeType(number_t i = 0) const; //!< return the shape type of element (i=0) and its sides (_segment, _triangle, ...) 253 number_t vertexNumber(number_t i) const; //!< return the global number of vertex i (i=1,...) 254 number_t nodeNumber(number_t i) const; //!< return the global number of node i (i=1,...) 255 std::vector<number_t> vertexNumbers(number_t s = 0) const; //!< return the global numbers of vertices if s=0, of side s else 256 std::vector<number_t> nodeNumbers(number_t s = 0) const; //!< return the global numbers of nodes if s=0, of side s else 257 string_t encodeElement(number_t s = 0) const; //!< string key of element or side element from ascending vertex numbers 258 string_t encodeElementNodes(number_t s = 0) const; //!< string key of element or side element from ascending node numbers clearGeomMapData()259 void clearGeomMapData() //! clear geometric map data (i.e geoMapData object) 260 {if(meshElement_p!=0) meshElement_p->clearGeomMapData();} 261 262 //various stuff 263 GeoNumPair elementSharingSide(number_t s) const; //!< return the adjacent element by side s as a pair of (GeomElement*, number_t) 264 std::vector<GeoNumPair> elementsSharingVertex(number_t v, bool o = false) const; //!< return the list of adjacent elements by vertex v 265 std::vector<GeoNumPair> elementsSharingSideOfSide(number_t s, bool o = false) const; //!< return the list of adjacent elements by side of side s 266 std::vector<GeomElement*> splitP1() const; //!< split element in first order elements and return as a list 267 // std::vector<GeomElement*> splitO1() const; //!< split element in first order elements and return as a list (unused) 268 real_t measure(const number_t sideNo = 0) const; //! return measure of element if side_no = 0, measure of side number side_no > 0 269 real_t size() const; //!< characteristic size of element (computed if not available) 270 Point center() const; //!< compute the center of element 271 bool contains(const std::vector<real_t>&); //!< true if point p as vector belongs to element 272 Point projection(const std::vector<real_t>&, real_t&); //!< projection of a point onto element 273 std::vector<real_t> normalVector(number_t =0) const; //!< return outward normal vector on side s (not normalized) or normal to element when s=0 274 std::vector<real_t> tangentVector(number_t) const; //!< return tangent vector on an edge (not normalized) 275 276 //print utilities 277 void print(std::ostream&) const; //!< print geomelement characteristics print(PrintStream & os) const278 void print(PrintStream& os) const {print(os.currentStream());} 279 friend std::ostream& operator<< (std::ostream&, const GeomElement&); 280 }; 281 282 std::ostream& operator<< (std::ostream&, const GeomElement&); //!< outputs geomelement characteristics 283 284 //! create the side index map of a list of elements 285 void createSideIndex(const std::vector<GeomElement*>& elements, 286 std::map<string_t,std::vector<GeoNumPair> >& sideIndex); 287 //! create the side index map of side elements of a list of elements 288 void createSideEltIndex(const std::vector<GeomElement*>& elements, 289 std::map<string_t, GeomElement* >& sideEltIndex); 290 291 //Overloaded operators == and < (to sort Lists), and << for output 292 bool operator== (const GeomElement&, const GeomElement&); //!<overload operator == to sort elements 293 bool operator< (const GeomElement&, const GeomElement&); //!<operator < to sort elements 294 295 /*! 296 \class GeomElementLess 297 a very small class defining a functor to sort std::vector of GeomElement* 298 ascendingly according to their number 299 */ 300 class GeomElementLess 301 { 302 public: operator ()(const GeomElement * x,const GeomElement * y)303 bool operator()(const GeomElement* x, const GeomElement* y) //! operator < to sort elements 304 {return x->number() < y->number(); } 305 }; // end of class GeomElementLess =============================================== 306 307 //------------------------------------------------------------------------------------ 308 /* 309 some utilities 310 */ 311 //------------------------------------------------------------------------------------ 312 real_t triangleArea(const Point&, const Point&, const Point&); //!< area of a triangle 313 real_t tetrahedronVolume(const Point&, const Point&, const Point&, const Point&); //!< volume of a tetrahedron 314 real_t defaultColoringRule(const GeomElement& gelt, const std::vector<real_t>& val); //!< default Geomelement coloring rule 315 } // end of namespace xlifepp 316 317 #endif /* GEOM_ELEMENT_HPP */ 318 319 320 //=============================== oldstuff of GeomElment class ============================================================= 321 322 //real_t* coords(const number_t ) const; //!<returns pointer to coordinates of point of local number 323 324 //Global numbering 325 // void minmaxPointNumbers (std::pair<int, int>&) const; //!<compute min and max of point numbers in element 326 // void checkPointNumbers (const std::pair<int, int>&, 327 // std::vector<bool>&) const; 328 // void computeNumberOfVertices (number_t&, 329 // std::vector<number_t>&); //!<compute number of vertices in elements not seen in previous elements 330 // number_t vertices (std::vector<Point>&, 331 // std::vector<number_t>&); //!<create vertex global numbering 332 // number_t sides(std::vector<GeomSide>&, 333 // std::vector<number_t>&, 334 // std::vector<std::vector<number_t> >&); //!<create edge global numbering 335 // 336 // number_t sideofsides(std::vector<GeomSide>&, std::vector<number_t>&, std::vector<std::vector<number_t> >&); 337 // //!<create face global numbering 338 // 339 // void reorderVertices(number_t sideNo, number_t* veReordered); //!<reorder face vertices on face (side_no = 1, ...) using vertex global numbering 340 // 341 // number_t sharedVertices(const MeshElement&,std::vector< std::pair<number_t,number_t> >&) const; 342 // //!<returns the number of vertices shared with another element and the pairs of local numbers of vertices in elements 343 // 344 // void reorderSharedVertices(number_t,std::vector< std::pair<number_t,number_t> >&); 345 // //!<sort pairs of local numbers of shared vertices according to ascending global numbering 346 // 347 // number_t sharedEdges(const MeshElement&,std::vector< pair<number_t,number_t> >&) const; 348 // //!<returns the number of sides shared with another element and the pairs of local numbers of edges in elements 349 // 350 // number_t sharedSideOfSides(const MeshElement&,std::vector< pair<number_t,number_t> >&) const; 351 // //!<returns the number of faces shared with another element and the pairs of local numbers of faces in elements 352 353 // bool isNeighbor(const MeshElement& ge) const; //!<is element a neighbor of another ? 354 // 355 // // Elementary computations 356 // void computeMeasure(); 357 // void computeMeasureOfSides(); 358 // void orientation(MeshElementStorage&); //!<computes element orientation 359 // void orientation(const real_t); //!<set element orientation 360 // void normalize(Vector<real_t>&) const; //!<normalize normal std::vector into a outward unit normal std::vector 361 // void nullNormalVector(real_t) const; 362