← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Solver/elasticityTerm.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-27 17:36:40 UTC
  • mfrom: (1.4.2 upstream)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20090927173640-meoeywl4f5hq5qas
Tags: 2.4.2.dfsg-1
[Christophe Prud'homme]
* New upstream release
  + solver code refactoring
  + better IDE integration.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Solver/elasticityTerm.cpp
 
1
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
 
2
//
 
3
// See the LICENSE.txt file for license information. Please report all
 
4
// bugs and problems to <gmsh@geuz.org>.
 
5
 
 
6
#include "elasticityTerm.h"
 
7
#include "Numeric.h"
 
8
 
 
9
void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
 
10
{
 
11
  MElement *e = se->getMeshElement();
 
12
  int nbNodes = e->getNumVertices();
 
13
  int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
 
14
  int npts;
 
15
  IntPt *GP;
 
16
  double jac[3][3];
 
17
  double invjac[3][3];
 
18
  double Grads[256][3], grads[100][3];
 
19
  e->getIntegrationPoints(integrationOrder, &npts, &GP);
 
20
  
 
21
  m.set_all(0.);
 
22
  
 
23
  double FACT = _E / (1 + _nu);
 
24
  double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
 
25
  double C12 = FACT * _nu / (1 - 2 * _nu);
 
26
  double C44 = (C11 - C12) / 2;
 
27
  const double C[6][6] =
 
28
    { {C11, C12, C12,    0,   0,   0}, 
 
29
      {C12, C11, C12,    0,   0,   0}, 
 
30
      {C12, C12, C11,    0,   0,   0}, 
 
31
      {  0,   0,   0,  C44,   0,   0},
 
32
      {  0,   0,   0,    0, C44,   0}, 
 
33
      {  0,   0,   0,    0,   0, C44} };
 
34
  
 
35
  fullMatrix<double> H(6, 6);
 
36
  fullMatrix<double> B(6, 3 * nbNodes);
 
37
  fullMatrix<double> BTH(3 * nbNodes, 6);
 
38
  fullMatrix<double> BT(3 * nbNodes, 6);
 
39
  for (int i = 0; i < 6; i++)
 
40
    for (int j = 0; j < 6; j++)
 
41
      H(i, j) = C[i][j];
 
42
  
 
43
  for (int i = 0; i < npts; i++){
 
44
    const double u = GP[i].pt[0];
 
45
    const double v = GP[i].pt[1];
 
46
    const double w = GP[i].pt[2];
 
47
    const double weight = GP[i].weight;
 
48
    const double detJ = e->getJacobian(u, v, w, jac);
 
49
    e->getGradShapeFunctions(u, v, w, grads);
 
50
    inv3x3(jac, invjac);
 
51
    B.set_all(0.);
 
52
    BT.set_all(0.);
 
53
    for (int j = 0; j < nbNodes; j++){
 
54
      Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] + 
 
55
        invjac[0][2] * grads[j][2];
 
56
      Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] + 
 
57
        invjac[1][2] * grads[j][2];
 
58
      Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] + 
 
59
        invjac[2][2] * grads[j][2];
 
60
      
 
61
      BT(j, 0) = B(0, j) = Grads[j][0];
 
62
      BT(j, 3) = B(3, j) = Grads[j][1];
 
63
      BT(j, 4) = B(4, j) = Grads[j][2];
 
64
      
 
65
      BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
 
66
      BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
 
67
      BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
 
68
      
 
69
      BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
 
70
      BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
 
71
      BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
 
72
    }
 
73
    BTH.set_all(0.);
 
74
    BTH.gemm(BT, H); 
 
75
    m.gemm(BTH, B, weight * detJ, 1.);
 
76
  } 
 
77
}
 
78
 
 
79
void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
 
80
{
 
81
  MElement *e = se->getMeshElement();
 
82
  int nbNodes = e->getNumVertices();
 
83
  int integrationOrder = 2 * e->getPolynomialOrder();
 
84
  int npts;
 
85
  IntPt *GP;
 
86
  double jac[3][3];
 
87
  double ff[256];
 
88
  e->getIntegrationPoints(integrationOrder, &npts, &GP);
 
89
  
 
90
  m.scale(0.);
 
91
  
 
92
  for (int i = 0; i < npts; i++){
 
93
    const double u = GP[i].pt[0];
 
94
    const double v = GP[i].pt[1];
 
95
    const double w = GP[i].pt[2];
 
96
    const double weight = GP[i].weight;
 
97
    const double detJ = e->getJacobian(u, v, w, jac);
 
98
    e->getShapeFunctions(u, v, w, ff);
 
99
    for (int j = 0; j < nbNodes; j++){
 
100
      m(j) += _volumeForce.x() *ff[j] * weight * detJ *.5;
 
101
      m(j + nbNodes) += _volumeForce.y() *ff[j] * weight * detJ *.5;
 
102
      m(j + 2 * nbNodes) += _volumeForce.z() *ff[j] * weight * detJ *.5;
 
103
    }
 
104
  } 
 
105
}