gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
terms.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 // Contributor(s):
7 // Eric Bechet
8 //
9 
10 #include "terms.h"
11 
12 void BilinearTermToScalarTerm::get(MElement *ele, int npts, IntPt *GP,
13  double &val) const
14 {
15  fullMatrix<double> localMatrix;
16  bilterm.get(ele, npts, GP, localMatrix);
17  val = localMatrix(0, 0);
18 }
19 
20 void BilinearTermBase::get(MElement *ele, int npts, IntPt *GP,
21  fullMatrix<double> &m) const
22 {
23  std::vector<fullMatrix<double> > mv(npts);
24  get(ele, npts, GP, mv);
25  m.resize(mv[0].size1(), mv[0].size2());
26  m.setAll(0.);
27  double jac[3][3];
28  for(int k = 0; k < npts; k++) {
29  const double u = GP[k].pt[0];
30  const double v = GP[k].pt[1];
31  const double w = GP[k].pt[2];
32  const double weight = GP[k].weight;
33  const double detJ = ele->getJacobian(u, v, w, jac);
34  const double coeff = weight * detJ;
35  for(int i = 0; i < mv[k].size1(); ++i)
36  for(int j = 0; j < mv[k].size2(); ++j) m(i, j) += mv[k](i, j) * coeff;
37  }
38 }
39 
40 IsotropicElasticTerm::IsotropicElasticTerm(FunctionSpace<SVector3> &space1_,
41  FunctionSpace<SVector3> &space2_,
42  double E_, double nu_)
43  : BilinearTerm<SVector3, SVector3>(space1_, space2_), E(E_), nu(nu_), H(6, 6)
44 {
45  double FACT = E / (1 + nu);
46  double C11 = FACT * (1 - nu) / (1 - 2 * nu);
47  double C12 = FACT * nu / (1 - 2 * nu);
48  double C44 = (C11 - C12) / 2;
49  H.scale(0.);
50  for(int i = 0; i < 3; ++i) {
51  H(i, i) = C11;
52  H(i + 3, i + 3) = C44;
53  }
54  H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
55  sym = (&space1_ == &space2_);
56 }
57 IsotropicElasticTerm::IsotropicElasticTerm(FunctionSpace<SVector3> &space1_,
58  double E_, double nu_)
59  : BilinearTerm<SVector3, SVector3>(space1_, space1_), E(E_), nu(nu_), H(6, 6)
60 {
61  double FACT = E / (1 + nu);
62  double C11 = FACT * (1 - nu) / (1 - 2 * nu);
63  double C12 = FACT * nu / (1 - 2 * nu);
64  double C44 = (C11 - C12) / 2;
65  /* FACT = E / (1 - nu * nu); // plane stress (plates)
66  C11 = FACT;
67  C12 = nu * FACT;
68  C44 = (1. - nu) * .5 * FACT;*/
69  H.scale(0.);
70  for(int i = 0; i < 3; ++i) {
71  H(i, i) = C11;
72  H(i + 3, i + 3) = C44;
73  }
74  H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
75  sym = true;
76 }
77 
78 void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP,
79  fullMatrix<double> &m) const
80 {
81  if(ele->getParent()) ele = ele->getParent();
82  if(sym) {
83  int nbFF = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
84  double jac[3][3];
85  fullMatrix<double> B(6, nbFF);
86  fullMatrix<double> BTH(nbFF, 6);
87  fullMatrix<double> BT(nbFF, 6);
88  m.resize(nbFF, nbFF);
89  m.setAll(0.);
90  // std::cout << m.size1() << " " << m.size2() << std::endl;
91  for(int i = 0; i < npts; i++) {
92  const double u = GP[i].pt[0];
93  const double v = GP[i].pt[1];
94  const double w = GP[i].pt[2];
95  const double weight = GP[i].weight;
96  const double detJ = ele->getJacobian(u, v, w, jac);
97  std::vector<TensorialTraits<SVector3>::GradType> Grads;
98  BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w,
99  Grads); // a optimiser ??
100  for(int j = 0; j < nbFF; j++) {
101  BT(j, 0) = B(0, j) = Grads[j](0, 0);
102  BT(j, 1) = B(1, j) = Grads[j](1, 1);
103  BT(j, 2) = B(2, j) = Grads[j](2, 2);
104  BT(j, 3) = B(3, j) = Grads[j](0, 1) + Grads[j](1, 0);
105  BT(j, 4) = B(4, j) = Grads[j](1, 2) + Grads[j](2, 1);
106  BT(j, 5) = B(5, j) = Grads[j](0, 2) + Grads[j](2, 0);
107  }
108  BTH.setAll(0.);
109  BTH.gemm(BT, H);
110  m.gemm(BTH, B, weight * detJ, 1.); // m = m + w*detJ*BT*H*B
111  }
112  }
113  else {
114  int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
115  int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele);
116  double jac[3][3];
117  fullMatrix<double> B(6, nbFF2);
118  fullMatrix<double> BTH(nbFF2, 6);
119  fullMatrix<double> BT(nbFF1, 6);
120  m.resize(nbFF1, nbFF2);
121  m.setAll(0.);
122  // Sum on Gauss Points i
123  for(int i = 0; i < npts; i++) {
124  const double u = GP[i].pt[0];
125  const double v = GP[i].pt[1];
126  const double w = GP[i].pt[2];
127  const double weight = GP[i].weight;
128  const double detJ = ele->getJacobian(u, v, w, jac);
129  std::vector<TensorialTraits<SVector3>::GradType>
130  Grads; // tableau de matrices...
131  std::vector<TensorialTraits<SVector3>::GradType>
132  GradsT; // tableau de matrices...
133  BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads);
134  BilinearTerm<SVector3, SVector3>::space2.gradf(ele, u, v, w, GradsT);
135  for(int j = 0; j < nbFF1; j++) {
136  BT(j, 0) = Grads[j](0, 0);
137  BT(j, 1) = Grads[j](1, 1);
138  BT(j, 2) = Grads[j](2, 2);
139  BT(j, 3) = Grads[j](0, 1) + Grads[j](1, 0);
140  BT(j, 4) = Grads[j](1, 2) + Grads[j](2, 1);
141  BT(j, 5) = Grads[j](0, 2) + Grads[j](2, 0);
142  }
143  for(int j = 0; j < nbFF2; j++) {
144  B(0, j) = GradsT[j](0, 0);
145  B(1, j) = GradsT[j](1, 1);
146  B(2, j) = GradsT[j](2, 2);
147  B(3, j) = GradsT[j](0, 1) + GradsT[j](1, 0);
148  B(4, j) = GradsT[j](1, 2) + GradsT[j](2, 1);
149  B(5, j) = GradsT[j](0, 2) + GradsT[j](2, 0);
150  }
151  BTH.setAll(0.);
152  BTH.gemm(BT, H);
153  // gemm add the product to m so there is a sum on gauss' points here
154  m.gemm(BTH, B, weight * detJ, 1.);
155  }
156  }
157 }
158 
159 void LagMultTerm::get(MElement *ele, int npts, IntPt *GP,
160  fullMatrix<double> &m) const
161 {
162  int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(
163  ele); // nbVertices*nbcomp of parent
164  int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(
165  ele); // nbVertices of boundary
166  double jac[3][3];
167  m.resize(nbFF1, nbFF2);
168  m.setAll(0.);
169  for(int i = 0; i < npts; i++) {
170  double u = GP[i].pt[0];
171  double v = GP[i].pt[1];
172  double w = GP[i].pt[2];
173  const double weight = GP[i].weight;
174  const double detJ = ele->getJacobian(u, v, w, jac);
175  std::vector<TensorialTraits<SVector3>::ValType> Vals;
176  std::vector<TensorialTraits<SVector3>::ValType> ValsT;
177  BilinearTerm<SVector3, SVector3>::space1.f(ele, u, v, w, Vals);
178  BilinearTerm<SVector3, SVector3>::space2.f(ele, u, v, w, ValsT);
179  for(int j = 0; j < nbFF1; j++) {
180  for(int k = 0; k < nbFF2; k++) {
181  m(j, k) += _eqfac * dot(Vals[j], ValsT[k]) * weight * detJ;
182  }
183  }
184  }
185 }
terms.h
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
IsotropicElasticTerm::E
double E
Definition: terms.h:291
fullMatrix::setAll
void setAll(const scalar &m)
Definition: fullMatrix.h:422
MElement::getParent
virtual MElement * getParent() const
Definition: MElement.h:231
fullMatrix::scale
void scale(const scalar s)
Definition: fullMatrix.h:447
SVector3
Definition: SVector3.h:16
LagMultTerm::get
virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix< double > &m) const
Definition: terms.cpp:159
IntPt::pt
double pt[3]
Definition: GaussIntegration.h:13
fullMatrix< double >
IntPt::weight
double weight
Definition: GaussIntegration.h:14
LagMultTerm::_eqfac
double _eqfac
Definition: terms.h:363
IsotropicElasticTerm::sym
bool sym
Definition: terms.h:292
IsotropicElasticTerm::nu
double nu
Definition: terms.h:291
BilinearTermToScalarTerm::get
virtual void get(MElement *ele, int npts, IntPt *GP, double &val) const
Definition: terms.cpp:12
MElement
Definition: MElement.h:30
IsotropicElasticTerm::IsotropicElasticTerm
IsotropicElasticTerm(FunctionSpace< SVector3 > &space1_, FunctionSpace< SVector3 > &space2_, double E_, double nu_)
Definition: terms.cpp:40
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
IsotropicElasticTerm::get
virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix< double > &m) const
Definition: terms.cpp:78
IsotropicElasticTerm::H
fullMatrix< double > H
Definition: terms.h:293
BilinearTermToScalarTerm::bilterm
BilinearTermBase & bilterm
Definition: terms.h:79
IntPt
Definition: GaussIntegration.h:12
BilinearTermBase::get
virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix< double > &m) const
Definition: terms.cpp:20
BilinearTerm
Definition: terms.h:159
FunctionSpace< SVector3 >
MElement::getJacobian
virtual double getJacobian(const fullMatrix< double > &gsf, double jac[3][3]) const
Definition: MElement.cpp:868
fullMatrix::resize
bool resize(int r, int c, bool resetValue=true)
Definition: fullMatrix.h:307