gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
elasticityTerm.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 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 #include "elasticityTerm.h"
7 #include "Numeric.h"
8 
9 // The SElement (Solver element) that has been sent to the function
10 // contains 2 enrichments, that can enrich both shape and test functions
11 
12 void elasticityTerm::createData(MElement *e) const
13 {
14  if(_data.find(e->getTypeForMSH()) != _data.end()) return;
15  elasticityDataAtGaussPoint d;
16  int nbSF = (int)e->getNumShapeFunctions();
17  int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
18  int npts;
19  IntPt *GP;
20  double gs[100][3];
21  e->getIntegrationPoints(integrationOrder, &npts, &GP);
22  for(int i = 0; i < npts; i++) {
23  fullMatrix<double> a(nbSF, 3);
24  const double u = GP[i].pt[0];
25  const double v = GP[i].pt[1];
26  const double w = GP[i].pt[2];
27  const double weight = GP[i].weight;
28  e->getGradShapeFunctions(u, v, w, gs);
29  for(int j = 0; j < nbSF; j++) {
30  a(j, 0) = gs[j][0];
31  a(j, 1) = gs[j][1];
32  a(j, 2) = gs[j][2];
33  }
34  d.gradSF.push_back(a);
35  d.u.push_back(u);
36  d.v.push_back(v);
37  d.w.push_back(w);
38  d.weight.push_back(weight);
39  }
40  _data[e->getTypeForMSH()] = d;
41 }
42 
43 void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
44 {
45  MElement *e = se->getMeshElement();
46  createData(e);
47 
48  int nbSF = (int)e->getNumShapeFunctions();
49  elasticityDataAtGaussPoint &d = _data[e->getTypeForMSH()];
50  int npts = d.u.size();
51  m.setAll(0.);
52 
53  double FACT = _e / (1 + _nu);
54  double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
55  double C12 = FACT * _nu / (1 - 2 * _nu);
56  double C44 = (C11 - C12) / 2;
57  const double C[6][6] = {{C11, C12, C12, 0, 0, 0}, {C12, C11, C12, 0, 0, 0},
58  {C12, C12, C11, 0, 0, 0}, {0, 0, 0, C44, 0, 0},
59  {0, 0, 0, 0, C44, 0}, {0, 0, 0, 0, 0, C44}};
60 
61  fullMatrix<double> H(6, 6);
62  fullMatrix<double> B(6, 3 * nbSF);
63  fullMatrix<double> BTH(3 * nbSF, 6);
64  fullMatrix<double> BT(3 * nbSF, 6);
65  for(int i = 0; i < 6; i++)
66  for(int j = 0; j < 6; j++) H(i, j) = C[i][j];
67 
68  double jac[3][3], invjac[3][3], Grads[100][3];
69  for(int i = 0; i < npts; i++) {
70  const double weight = d.weight[i];
71  const fullMatrix<double> &grads = d.gradSF[i];
72  const double detJ = e->getJacobian(grads, jac);
73  inv3x3(jac, invjac);
74 
75  for(int j = 0; j < nbSF; j++) {
76  Grads[j][0] = invjac[0][0] * grads(j, 0) + invjac[0][1] * grads(j, 1) +
77  invjac[0][2] * grads(j, 2);
78  Grads[j][1] = invjac[1][0] * grads(j, 0) + invjac[1][1] * grads(j, 1) +
79  invjac[1][2] * grads(j, 2);
80  Grads[j][2] = invjac[2][0] * grads(j, 0) + invjac[2][1] * grads(j, 1) +
81  invjac[2][2] * grads(j, 2);
82  }
83 
84  B.setAll(0.);
85  BT.setAll(0.);
86 
87  if(se->getShapeEnrichement() == se->getTestEnrichement()) {
88  for(int j = 0; j < nbSF; j++) {
89  BT(j, 0) = B(0, j) = Grads[j][0];
90  BT(j, 3) = B(3, j) = Grads[j][1];
91  BT(j, 5) = B(5, j) = Grads[j][2];
92 
93  BT(j + nbSF, 1) = B(1, j + nbSF) = Grads[j][1];
94  BT(j + nbSF, 3) = B(3, j + nbSF) = Grads[j][0];
95  BT(j + nbSF, 4) = B(4, j + nbSF) = Grads[j][2];
96 
97  BT(j + 2 * nbSF, 2) = B(2, j + 2 * nbSF) = Grads[j][2];
98  BT(j + 2 * nbSF, 4) = B(4, j + 2 * nbSF) = Grads[j][1];
99  BT(j + 2 * nbSF, 5) = B(5, j + 2 * nbSF) = Grads[j][0];
100  }
101  }
102  else {
103  /*
104  se->gradNodalTestFunctions (u, v, w, invjac,GradsT);
105  for (int j = 0; j < nbSF; j++){
106  BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
107  BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1];
108  BT(j, 5) = Grads[j][2]; B(5, j) = GradsT[j][2];
109 
110  BT(j + nbSF, 1) = Grads[j][1]; B(1, j + nbSF) = GradsT[j][1];
111  BT(j + nbSF, 3) = Grads[j][0]; B(3, j + nbSF) = GradsT[j][0];
112  BT(j + nbSF, 4) = Grads[j][2]; B(4, j + nbSF) = GradsT[j][2];
113 
114  BT(j + 2 * nbSF, 2) = Grads[j][2]; B(2, j + 2 * nbSF) = GradsT[j][2];
115  BT(j + 2 * nbSF, 4) = Grads[j][1]; B(4, j + 2 * nbSF) = GradsT[j][1];
116  BT(j + 2 * nbSF, 5) = Grads[j][0]; B(5, j + 2 * nbSF) = GradsT[j][0];
117  }
118  */
119  }
120 
121  BTH.setAll(0.);
122  BTH.gemm(BT, H);
123  m.gemm(BTH, B, weight * detJ, 1.);
124  }
125 }
126 
127 void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
128 {
129  MElement *e = se->getMeshElement();
130  int nbSF = (int)e->getNumShapeFunctions();
131  int integrationOrder = 2 * e->getPolynomialOrder();
132  int npts;
133  IntPt *GP;
134  double jac[3][3];
135  double ff[256];
136  e->getIntegrationPoints(integrationOrder, &npts, &GP);
137 
138  m.scale(0.);
139 
140  for(int i = 0; i < npts; i++) {
141  const double u = GP[i].pt[0];
142  const double v = GP[i].pt[1];
143  const double w = GP[i].pt[2];
144  const double weight = GP[i].weight;
145  const double detJ = e->getJacobian(u, v, w, jac);
146  se->nodalTestFunctions(u, v, w, ff);
147  for(int j = 0; j < nbSF; j++) {
148  m(j) += _volumeForce.x() * ff[j] * weight * detJ * .5;
149  m(j + nbSF) += _volumeForce.y() * ff[j] * weight * detJ * .5;
150  m(j + 2 * nbSF) += _volumeForce.z() * ff[j] * weight * detJ * .5;
151  }
152  }
153 }
154 
155 /*
156 
157 B_d is the deviatoric part of the FE strains
158 
159 K_{uu} = \int_\Omega B^T_d H B_d dv
160 
161 */
162 
163 void elasticityMixedTerm::elementMatrix(SElement *se,
164  fullMatrix<double> &m) const
165 {
166  setPolynomialBasis(se);
167  MElement *e = se->getMeshElement();
168  int integrationOrder = 4 * (e->getPolynomialOrder()) + 2;
169  int npts;
170  IntPt *GP;
171  double jac[3][3];
172  double invjac[3][3];
173  double Grads[100][3];
174  double SF[100];
175  e->getIntegrationPoints(integrationOrder, &npts, &GP);
176 
177  m.setAll(0.);
178  fullMatrix<double> KUU(3 * _sizeN, 3 * _sizeN);
179  fullMatrix<double> KUP(3 * _sizeN, _sizeM);
180  fullMatrix<double> KUG(3 * _sizeN, _sizeM);
181  fullMatrix<double> KPG(_sizeM, _sizeM);
182  fullMatrix<double> KGG(_sizeM, _sizeM);
183  fullMatrix<double> KPP(_sizeM, _sizeM); // stabilization
184 
185  double FACT = _e / ((1 + _nu) * (1. - 2 * _nu));
186  double C11 = FACT * (1 - _nu);
187  double C12 = FACT * _nu;
188  double C44 = FACT * (1 - 2. * _nu) * .5;
189  // double _k = 3*_e / (1 - 2 * _nu);
190  // FIXME : PLANE STRESS !!!
191  // FACT = _e / (1-_nu*_nu);
192  // C11 = FACT;
193  // C12 = _nu * FACT;
194  // C44 = (1.-_nu)*.5*FACT;
195 
196  const double C[6][6] = {{C11, C12, C12, 0, 0, 0}, {C12, C11, C12, 0, 0, 0},
197  {C12, C12, C11, 0, 0, 0}, {0, 0, 0, C44, 0, 0},
198  {0, 0, 0, 0, C44, 0}, {0, 0, 0, 0, 0, C44}};
199 
200  fullMatrix<double> H(6, 6);
201  fullMatrix<double> B(6, 3 * _sizeN);
202  fullMatrix<double> BDEV(6, 3 * _sizeN);
203  fullMatrix<double> BDIL(1, 3 * _sizeN);
204  fullMatrix<double> BDILT(3 * _sizeN, 1);
205  fullMatrix<double> BTH(3 * _sizeN, 6);
206  fullMatrix<double> BT(3 * _sizeN, 6);
207  fullMatrix<double> trBTH(3 * _sizeN, 1);
208  fullMatrix<double> MT(_sizeM, 1);
209  fullMatrix<double> M(1, _sizeM);
210  for(int i = 0; i < 6; i++)
211  for(int j = 0; j < 6; j++) H(i, j) = C[i][j];
212 
213  double _dimension = 2.0;
214 
215  for(int i = 0; i < npts; i++) {
216  const double u = GP[i].pt[0];
217  const double v = GP[i].pt[1];
218  const double w = GP[i].pt[2];
219  const double weight = GP[i].weight;
220  const double detJ = e->getJacobian(u, v, w, jac);
221  inv3x3(jac, invjac);
222  se->gradNodalShapeFunctions(u, v, w, invjac, Grads);
223  e->getShapeFunctions(u, v, w, SF, _polyOrderM);
224 
225  B.setAll(0.);
226  BT.setAll(0.);
227 
228  const double K =
229  .0000002 * weight * detJ *
230  pow(detJ, 2 / _dimension); // the second detJ is for stabilization
231 
232  for(int j = 0; j < _sizeM; j++) {
233  for(int k = 0; k < _sizeM; k++) {
234  KPP(j, k) += (Grads[j][0] * Grads[k][0] + Grads[j][1] * Grads[k][1] +
235  Grads[j][2] * Grads[k][2]) *
236  K;
237  }
238  }
239 
240  const double K2 = weight * detJ;
241 
242  for(int j = 0; j < _sizeN; j++) {
243  for(int k = 0; k < _sizeM; k++) {
244  KUP(j + 0 * _sizeN, k) += Grads[j][0] * SF[k] * K2;
245  KUP(j + 1 * _sizeN, k) += Grads[j][1] * SF[k] * K2;
246  KUP(j + 2 * _sizeN, k) += Grads[j][2] * SF[k] * K2;
247  }
248  }
249 
250  /*
251  const double K3 = weight * detJ * _e/(2*(1+_nu));
252  for (int j = 0; j < _sizeN; j++){
253  for (int k = 0; k < _sizeN; k++){
254  const double fa = (Grads[j][0] * Grads[k][0] +
255  Grads[j][1] * Grads[k][1] +
256  Grads[j][2] * Grads[k][2]) * K3;
257  KUU(j+0*_sizeN, k+0*_sizeN) += fa;
258  KUU(j+1*_sizeN, k+1*_sizeN) += fa;
259  KUU(j+2*_sizeN, k+2*_sizeN) += fa;
260  }
261  }
262  */
263 
264  for(int j = 0; j < _sizeN; j++) {
265  // printf("%g %g %g\n",Grads[j][0],Grads[j][1],Grads[j][2]);
266 
267  BDILT(j, 0) = BDIL(0, j) = Grads[j][0] / _dimension;
268 
269  BT(j, 0) = B(0, j) = Grads[j][0]; //-Grads[j][0]/_dimension;
270  // BT(j, 1) = B(1, j) = -Grads[j][1]/_dimension;
271  // BT(j, 2) = B(2, j) = -Grads[j][2]/_dimension;
272 
273  BT(j, 3) = B(3, j) = Grads[j][1];
274  BT(j, 4) = B(4, j) = Grads[j][2];
275 
276  BDILT(j + _sizeN, 0) = BDIL(0, j + _sizeN) = Grads[j][1] / _dimension;
277 
278  // BT(j + _sizeN, 0) = B(0, j + _sizeN) = -Grads[j][0]/_dimension;
279  BT(j + _sizeN, 1) = B(1, j + _sizeN) =
280  Grads[j][1]; //-Grads[j][1]/_dimension;
281  // BT(j + _sizeN, 2) = B(2, j + _sizeN) = -Grads[j][2]/_dimension;
282 
283  BT(j + _sizeN, 3) = B(3, j + _sizeN) = Grads[j][0];
284  BT(j + _sizeN, 5) = B(5, j + _sizeN) = Grads[j][2];
285 
286  BDILT(j + 2 * _sizeN, 0) = BDIL(0, j + 2 * _sizeN) =
287  Grads[j][2] / _dimension;
288 
289  // BT(j + 2 * _sizeN, 0) = B(0, j + 2 * _sizeN) =
290  // -Grads[j][0]/_dimension; BT(j + 2 * _sizeN, 1) = B(1, j + 2 *
291  // _sizeN) = -Grads[j][1]/_dimension;
292  BT(j + 2 * _sizeN, 2) = B(2, j + 2 * _sizeN) =
293  Grads[j][2]; //-Grads[j][2]/_dimension;
294 
295  BT(j + 2 * _sizeN, 4) = B(4, j + 2 * _sizeN) = Grads[j][1];
296  BT(j + 2 * _sizeN, 5) = B(5, j + 2 * _sizeN) = Grads[j][0];
297  }
298 
299  for(int j = 0; j < _sizeM; j++) {
300  MT(j, 0) = M(0, j) = SF[j];
301  }
302 
303  // printf("BT (%d %d) x H (%d,%d) = BTH(%d,%d)\n",BT.size1(),BT.size2(),
304  // H.size1(),H.size2(),BTH.size1(),BTH.size2());
305 
306  BTH.setAll(0.);
307  BTH.gemm(BT, H);
308 
309  for(int j = 0; j < 3 * _sizeN; j++) {
310  trBTH(j, 0) = BTH(j, 0) + BTH(j, 1) + BTH(j, 2);
311  }
312 
313  KUU.gemm(BTH, B, weight * detJ);
314  // KUP.gemm(BDILT, M, weight * detJ);
315  // KPG.gemm(MT, M, -weight * detJ);
316  // KUG.gemm(trBTH, M, weight * detJ/_dimension);
317  // KGG.gemm(MT, M, weight * detJ * (_k)/(_dimension*_dimension));
318  }
319  /*
320  3*_sizeN _sizeM _sizeM
321 
322  KUU KUP KUG
323  m = KUP^t KPP KPG
324  KUG^T KPG^t KGG
325 
326  */
327 
328  for(int i = 0; i < 3 * _sizeN; i++) {
329  for(int j = 0; j < 3 * _sizeN; j++) m(i, j) = KUU(i, j); // assemble KUU
330 
331  for(int j = 0; j < _sizeM; j++) {
332  m(i, j + 3 * _sizeN) = KUP(i, j); // assemble KUP
333  m(j + 3 * _sizeN, i) = KUP(i, j); // assemble KUP^T
334  }
335 
336  for(int j = 0; j < _sizeM; j++) {
337  m(i, j + 3 * _sizeN + _sizeM) = KUG(i, j); // assemble KUG
338  m(j + 3 * _sizeN + _sizeM, i) = KUG(i, j); // assemble KUG^t
339  }
340  }
341  for(int i = 0; i < _sizeM; i++) {
342  for(int j = 0; j < _sizeM; j++) {
343  m(i + 3 * _sizeN, j + 3 * _sizeN + _sizeM) = KPG(i, j); // assemble KPG
344  m(j + 3 * _sizeN + _sizeM, i + 3 * _sizeN) = KPG(i, j); // assemble KPG^t
345 
346  m(i + 3 * _sizeN + _sizeM, j + 3 * _sizeN + _sizeM) =
347  KGG(i, j); // assemble KGG
348 
349  m(i + 3 * _sizeN, j + 3 * _sizeN) = KPP(i, j); // assemble KPP
350  }
351  }
352  // m.print("Mixed");
353  return;
354 }
elasticityMixedTerm::elementMatrix
void elementMatrix(SElement *se, fullMatrix< double > &m) const
Definition: elasticityTerm.cpp:163
inv3x3
double inv3x3(double mat[3][3], double inv[3][3])
Definition: Numeric.cpp:232
elasticityTerm::createData
void createData(MElement *) const
Definition: elasticityTerm.cpp:12
elasticityMixedTerm::_sizeM
int _sizeM
Definition: elasticityTerm.h:88
elasticityMixedTerm::_polyOrderM
int _polyOrderM
Definition: elasticityTerm.h:88
fullVector< double >
MElement::getTypeForMSH
virtual int getTypeForMSH() const
Definition: MElement.h:488
MElement::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MElement.h:436
fullMatrix::setAll
void setAll(const scalar &m)
Definition: fullMatrix.h:422
fullVector::scale
void scale(const scalar s)
Definition: fullMatrix.h:144
SElement::nodalTestFunctions
void nodalTestFunctions(double u, double v, double w, double s[])
Definition: SElement.cpp:96
IntPt::pt
double pt[3]
Definition: GaussIntegration.h:13
elasticityTerm::_nu
double _nu
Definition: elasticityTerm.h:23
elasticityTerm::_volumeForce
SVector3 _volumeForce
Definition: elasticityTerm.h:25
elasticityDataAtGaussPoint::u
std::vector< double > u
Definition: elasticityTerm.h:18
elasticityDataAtGaussPoint::weight
std::vector< double > weight
Definition: elasticityTerm.h:18
fullMatrix< double >
IntPt::weight
double weight
Definition: GaussIntegration.h:14
elasticityTerm::elementMatrix
void elementMatrix(SElement *se, fullMatrix< double > &m) const
Definition: elasticityTerm.cpp:43
elasticityDataAtGaussPoint::v
std::vector< double > v
Definition: elasticityTerm.h:18
elasticityDataAtGaussPoint::w
std::vector< double > w
Definition: elasticityTerm.h:18
Numeric.h
elasticityMixedTerm::_sizeN
int _sizeN
Definition: elasticityTerm.h:88
MElement::getNumShapeFunctions
virtual std::size_t getNumShapeFunctions() const
Definition: MElement.h:387
elasticityTerm::elementVector
void elementVector(SElement *se, fullVector< double > &m) const
Definition: elasticityTerm.cpp:127
elasticityTerm::_e
double _e
Definition: elasticityTerm.h:23
elasticityMixedTerm::setPolynomialBasis
void setPolynomialBasis(SElement *se) const
Definition: elasticityTerm.h:90
SElement::getShapeEnrichement
static const simpleFunction< double > * getShapeEnrichement()
Definition: SElement.h:44
SElement
Definition: SElement.h:18
elasticityDataAtGaussPoint
Definition: elasticityTerm.h:16
SElement::getTestEnrichement
static const simpleFunction< double > * getTestEnrichement()
Definition: SElement.h:48
MElement
Definition: MElement.h:30
SVector3::x
double x(void) const
Definition: SVector3.h:30
fullMatrix::gemm
void gemm(const fullMatrix< scalar > &a, const fullMatrix< scalar > &b, scalar alpha=1., scalar beta=1., bool transposeA=false, bool transposeB=false)
Definition: fullMatrix.h:580
SElement::getMeshElement
MElement * getMeshElement() const
Definition: SElement.h:35
elasticityMixedTerm::_e
double _e
Definition: elasticityTerm.h:87
elasticityMixedTerm::_nu
double _nu
Definition: elasticityTerm.h:87
IntPt
Definition: GaussIntegration.h:12
SVector3::y
double y(void) const
Definition: SVector3.h:31
SVector3::z
double z(void) const
Definition: SVector3.h:32
MElement::getShapeFunctions
virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const
Definition: MElement.cpp:458
elasticityTerm::_data
std::map< int, elasticityDataAtGaussPoint > _data
Definition: elasticityTerm.h:26
MElement::getJacobian
virtual double getJacobian(const fullMatrix< double > &gsf, double jac[3][3]) const
Definition: MElement.cpp:868
elasticityTerm.h
SElement::gradNodalShapeFunctions
void gradNodalShapeFunctions(double u, double v, double w, double invjac[3][3], double grad[][3])
Definition: SElement.cpp:79
MElement::getPolynomialOrder
virtual int getPolynomialOrder() const
Definition: MElement.h:78
elasticityDataAtGaussPoint::gradSF
std::vector< fullMatrix< double > > gradSF
Definition: elasticityTerm.h:17
MElement::getGradShapeFunctions
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1) const
Definition: MElement.cpp:468