1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle 2 // 3 // See the LICENSE.txt file in the Gmsh root directory for license information. 4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues. 5 6 #ifndef ELASTICITY_TERM_H 7 #define ELASTICITY_TERM_H 8 9 #include "femTerm.h" 10 #include "GmshGlobal.h" 11 #include "GModel.h" 12 #include "polynomialBasis.h" 13 #include "SElement.h" 14 #include "fullMatrix.h" 15 16 struct elasticityDataAtGaussPoint { 17 std::vector<fullMatrix<double> > gradSF; 18 std::vector<double> u, v, w, weight; 19 }; 20 21 class elasticityTerm : public femTerm<double> { 22 protected: 23 double _e, _nu; 24 int _iFieldR, _iFieldC; 25 SVector3 _volumeForce; 26 mutable std::map<int, elasticityDataAtGaussPoint> _data; 27 void createData(MElement *) const; 28 29 public: setFieldC(int i)30 void setFieldC(int i) { _iFieldC = i; } setFieldR(int i)31 void setFieldR(int i) { _iFieldR = i; } 32 // element matrix size : 3 dofs per vertex sizeOfR(SElement * se)33 virtual int sizeOfR(SElement *se) const 34 { 35 return 3 * se->getMeshElement()->getNumShapeFunctions(); 36 } sizeOfC(SElement * se)37 virtual int sizeOfC(SElement *se) const 38 { 39 return 3 * se->getMeshElement()->getNumShapeFunctions(); 40 } 41 // order dofs in the local matrix : 42 // dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn getLocalDofR(SElement * se,int iRow)43 Dof getLocalDofR(SElement *se, int iRow) const 44 { 45 MElement *e = se->getMeshElement(); 46 int iCompR = iRow / (int)e->getNumShapeFunctions(); 47 int ithLocalVertex = iRow % (int)e->getNumShapeFunctions(); 48 return Dof(e->getShapeFunctionNode(ithLocalVertex)->getNum(), 49 Dof::createTypeWithTwoInts(iCompR, _iFieldR)); 50 } getLocalDofC(SElement * se,int iCol)51 Dof getLocalDofC(SElement *se, int iCol) const 52 { 53 MElement *e = se->getMeshElement(); 54 int iCompC = iCol / (int)e->getNumShapeFunctions(); 55 int ithLocalVertex = iCol % (int)e->getNumShapeFunctions(); 56 return Dof(e->getShapeFunctionNode(ithLocalVertex)->getNum(), 57 Dof::createTypeWithTwoInts(iCompC, _iFieldC)); 58 } 59 60 public: elasticityTerm(GModel * gm,double E,double nu,int fieldr,int fieldc)61 elasticityTerm(GModel *gm, double E, double nu, int fieldr, int fieldc) 62 : femTerm<double>(gm), _e(E), _nu(nu), _iFieldR(fieldr), _iFieldC(fieldc) 63 { 64 } elasticityTerm(GModel * gm,double E,double nu,int fieldr)65 elasticityTerm(GModel *gm, double E, double nu, int fieldr) 66 : femTerm<double>(gm), _e(E), _nu(nu), _iFieldR(fieldr), _iFieldC(fieldr) 67 { 68 } setVector(const SVector3 & f)69 void setVector(const SVector3 &f) { _volumeForce = f; } 70 void elementMatrix(SElement *se, fullMatrix<double> &m) const; 71 void elementVector(SElement *se, fullVector<double> &m) const; 72 }; 73 74 /* 75 Formulation of elasticity with 3 fields 76 -) displacement U 77 -) lagrange multipliers p and g 78 79 g = trace (epsilon) 80 p = trace (epsilon) 81 82 83 */ 84 85 class elasticityMixedTerm : public femTerm<double> { 86 protected: 87 double _e, _nu; 88 mutable int _iField, _polyOrderN, _polyOrderM, _sizeN, _sizeM; 89 mutable polynomialBasis *_pN, *_pM; setPolynomialBasis(SElement * se)90 void setPolynomialBasis(SElement *se) const 91 { 92 _polyOrderN = se->getMeshElement()->getPolynomialOrder(); 93 _polyOrderM = se->getMeshElement()->getPolynomialOrder(); 94 _pN = 95 (polynomialBasis *)se->getMeshElement()->getFunctionSpace(_polyOrderN); 96 _pM = 97 (polynomialBasis *)se->getMeshElement()->getFunctionSpace(_polyOrderM); 98 _sizeN = _pN->coefficients.size1(); 99 _sizeM = _pM->coefficients.size1(); 100 } 101 102 public: setField(int i)103 void setField(int i) { _iField = i; } 104 // element matrix size : 3 dofs per vertex sizeOfR(SElement * se)105 virtual int sizeOfR(SElement *se) const 106 { 107 setPolynomialBasis(se); 108 return 3 * _sizeN + 2 * _sizeM; 109 } sizeOfC(SElement * se)110 virtual int sizeOfC(SElement *se) const { return sizeOfR(se); } 111 // order dofs in the local matrix : 112 // dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn , 113 // p1,p2,...,pm, g1,g2,...,gm getLocalDofR(SElement * se,int iRow)114 Dof getLocalDofR(SElement *se, int iRow) const 115 { 116 setPolynomialBasis(se); 117 MElement *e = se->getMeshElement(); 118 int iComp; 119 int ithLocalVertex; 120 if(iRow < 3 * _sizeN) { 121 iComp = iRow / _sizeN; 122 ithLocalVertex = iRow % _sizeN; 123 } 124 else { 125 iRow -= 3 * _sizeN; 126 iComp = 3 + iRow / _sizeM; 127 ithLocalVertex = iRow % _sizeM; 128 } 129 return Dof(e->getShapeFunctionNode(ithLocalVertex)->getNum(), 130 Dof::createTypeWithTwoInts(iComp, _iField)); 131 } getLocalDofC(SElement * se,int iCol)132 Dof getLocalDofC(SElement *se, int iCol) const 133 { 134 return getLocalDofR(se, iCol); 135 } 136 137 public: elasticityMixedTerm(GModel * gm,double E,double nu,int field)138 elasticityMixedTerm(GModel *gm, double E, double nu, int field) 139 : femTerm<double>(gm), _e(E), _nu(nu), _iField(field) 140 { 141 } 142 void elementMatrix(SElement *se, fullMatrix<double> &m) const; elementVector(SElement * se,fullVector<double> & m)143 void elementVector(SElement *se, fullVector<double> &m) const {} setYoung(double E)144 void setYoung(double E) { _e = E; } 145 }; 146 147 #endif 148