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