MFEM  v3.3
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nonlininteg.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "fem.hpp"
13 
14 namespace mfem
15 {
16 
18  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
19  DenseMatrix &elmat)
20 {
21  mfem_error("NonlinearFormIntegrator::AssembleElementGrad"
22  " is not overloaded!");
23 }
24 
26  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
27 {
28  mfem_error("NonlinearFormIntegrator::GetElementEnergy"
29  " is not overloaded!");
30  return 0.0;
31 }
32 
33 
35 {
36  Z.SetSize(J.Width());
38  return 0.5*(Z*Z)/J.Det();
39 }
40 
42 {
43  int dim = J.Width();
44  double t;
45 
46  Z.SetSize(dim);
47  S.SetSize(dim);
49  MultAAt(Z, S);
50  t = 0.5*S.Trace();
51  for (int i = 0; i < dim; i++)
52  {
53  S(i,i) -= t;
54  }
55  t = J.Det();
56  S *= -1.0/(t*t);
57  Mult(S, Z, P);
58 }
59 
61  const DenseMatrix &J, const DenseMatrix &DS, const double weight,
62  DenseMatrix &A) const
63 {
64  int dof = DS.Height(), dim = DS.Width();
65  double t;
66 
67  Z.SetSize(dim);
68  S.SetSize(dim);
69  G.SetSize(dof, dim);
70  C.SetSize(dof, dim);
71 
73  MultAAt(Z, S);
74 
75  t = 1.0/J.Det();
76  Z *= t; // Z = J^{-t}
77  S *= t; // S = |J| (J.J^t)^{-1}
78  t = 0.5*S.Trace();
79 
80  MultABt(DS, Z, G); // G = DS.J^{-1}
81  Mult(G, S, C);
82 
83  // 1.
84  for (int i = 0; i < dof; i++)
85  for (int j = 0; j <= i; j++)
86  {
87  double a = 0.0;
88  for (int d = 0; d < dim; d++)
89  {
90  a += G(i,d)*G(j,d);
91  }
92  a *= weight;
93  for (int k = 0; k < dim; k++)
94  for (int l = 0; l <= k; l++)
95  {
96  double b = a*S(k,l);
97  A(i+k*dof,j+l*dof) += b;
98  if (i != j)
99  {
100  A(j+k*dof,i+l*dof) += b;
101  }
102  if (k != l)
103  {
104  A(i+l*dof,j+k*dof) += b;
105  if (i != j)
106  {
107  A(j+l*dof,i+k*dof) += b;
108  }
109  }
110  }
111  }
112 
113  // 2.
114  for (int i = 0; i < dof; i++)
115  for (int j = 0; j < i; j++)
116  {
117  for (int k = 0; k < dim; k++)
118  for (int l = 0; l < k; l++)
119  {
120  double a =
121  weight*(C(i,l)*G(j,k) - C(i,k)*G(j,l) +
122  C(j,k)*G(i,l) - C(j,l)*G(i,k) +
123  t*(G(i,k)*G(j,l) - G(i,l)*G(j,k)));
124 
125  A(i+k*dof,j+l*dof) += a;
126  A(j+l*dof,i+k*dof) += a;
127 
128  A(i+l*dof,j+k*dof) -= a;
129  A(j+k*dof,i+l*dof) -= a;
130  }
131  }
132 }
133 
134 
135 inline void NeoHookeanModel::EvalCoeffs() const
136 {
137  mu = c_mu->Eval(*T, T->GetIntPoint());
138  K = c_K->Eval(*T, T->GetIntPoint());
139  if (c_g)
140  {
141  g = c_g->Eval(*T, T->GetIntPoint());
142  }
143 }
144 
145 double NeoHookeanModel::EvalW(const DenseMatrix &J) const
146 {
147  int dim = J.Width();
148 
149  if (have_coeffs)
150  {
151  EvalCoeffs();
152  }
153 
154  double dJ = J.Det();
155  double sJ = dJ/g;
156  double bI1 = pow(dJ, -2.0/dim)*(J*J); // \bar{I}_1
157 
158  return 0.5*(mu*(bI1 - dim) + K*(sJ - 1.0)*(sJ - 1.0));
159 }
160 
162 {
163  int dim = J.Width();
164 
165  if (have_coeffs)
166  {
167  EvalCoeffs();
168  }
169 
170  Z.SetSize(dim);
172 
173  double dJ = J.Det();
174  double a = mu*pow(dJ, -2.0/dim);
175  double b = K*(dJ/g - 1.0)/g - a*(J*J)/(dim*dJ);
176 
177  P = 0.0;
178  P.Add(a, J);
179  P.Add(b, Z);
180 }
181 
183  const double weight, DenseMatrix &A) const
184 {
185  int dof = DS.Height(), dim = DS.Width();
186 
187  if (have_coeffs)
188  {
189  EvalCoeffs();
190  }
191 
192  Z.SetSize(dim);
193  G.SetSize(dof, dim);
194  C.SetSize(dof, dim);
195 
196  double dJ = J.Det();
197  double sJ = dJ/g;
198  double a = mu*pow(dJ, -2.0/dim);
199  double bc = a*(J*J)/dim;
200  double b = bc - K*sJ*(sJ - 1.0);
201  double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0);
202 
204  Z *= (1.0/dJ); // Z = J^{-t}
205 
206  MultABt(DS, J, C); // C = DS J^t
207  MultABt(DS, Z, G); // G = DS J^{-1}
208 
209  a *= weight;
210  b *= weight;
211  c *= weight;
212 
213  // 1.
214  for (int i = 0; i < dof; i++)
215  for (int k = 0; k <= i; k++)
216  {
217  double s = 0.0;
218  for (int d = 0; d < dim; d++)
219  {
220  s += DS(i,d)*DS(k,d);
221  }
222  s *= a;
223 
224  for (int d = 0; d < dim; d++)
225  {
226  A(i+d*dof,k+d*dof) += s;
227  }
228 
229  if (k != i)
230  for (int d = 0; d < dim; d++)
231  {
232  A(k+d*dof,i+d*dof) += s;
233  }
234  }
235 
236  a *= (-2.0/dim);
237 
238  // 2.
239  for (int i = 0; i < dof; i++)
240  for (int j = 0; j < dim; j++)
241  for (int k = 0; k < dof; k++)
242  for (int l = 0; l < dim; l++)
243  {
244  A(i+j*dof,k+l*dof) +=
245  a*(C(i,j)*G(k,l) + G(i,j)*C(k,l)) +
246  b*G(i,l)*G(k,j) + c*G(i,j)*G(k,l);
247  }
248 }
249 
250 
253  const Vector &elfun)
254 {
255  int dof = el.GetDof(), dim = el.GetDim();
256  double energy;
257 
258  DSh.SetSize(dof, dim);
259  J0i.SetSize(dim);
260  J1.SetSize(dim);
261  J.SetSize(dim);
262  PMatI.UseExternalData(elfun.GetData(), dof, dim);
263 
264  int intorder = 2*el.GetOrder() + 3; // <---
265  const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
266 
267  energy = 0.0;
268  model->SetTransformation(Tr);
269  for (int i = 0; i < ir.GetNPoints(); i++)
270  {
271  const IntegrationPoint &ip = ir.IntPoint(i);
272  Tr.SetIntPoint(&ip);
273  CalcInverse(Tr.Jacobian(), J0i);
274 
275  el.CalcDShape(ip, DSh);
276  MultAtB(PMatI, DSh, J1);
277  Mult(J1, J0i, J);
278 
279  energy += ip.weight*Tr.Weight()*model->EvalW(J);
280  }
281 
282  return energy;
283 }
284 
286  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
287  Vector &elvect)
288 {
289  int dof = el.GetDof(), dim = el.GetDim();
290 
291  DSh.SetSize(dof, dim);
292  DS.SetSize(dof, dim);
293  J0i.SetSize(dim);
294  J.SetSize(dim);
295  P.SetSize(dim);
296  PMatI.UseExternalData(elfun.GetData(), dof, dim);
297  elvect.SetSize(dof*dim);
298  PMatO.UseExternalData(elvect.GetData(), dof, dim);
299 
300  int intorder = 2*el.GetOrder() + 3; // <---
301  const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
302 
303  elvect = 0.0;
304  model->SetTransformation(Tr);
305  for (int i = 0; i < ir.GetNPoints(); i++)
306  {
307  const IntegrationPoint &ip = ir.IntPoint(i);
308  Tr.SetIntPoint(&ip);
309  CalcInverse(Tr.Jacobian(), J0i);
310 
311  el.CalcDShape(ip, DSh);
312  Mult(DSh, J0i, DS);
313  MultAtB(PMatI, DS, J);
314 
315  model->EvalP(J, P);
316 
317  P *= ip.weight*Tr.Weight();
318  AddMultABt(DS, P, PMatO);
319  }
320 }
321 
323  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
324  DenseMatrix &elmat)
325 {
326  int dof = el.GetDof(), dim = el.GetDim();
327 
328  DSh.SetSize(dof, dim);
329  DS.SetSize(dof, dim);
330  J0i.SetSize(dim);
331  J.SetSize(dim);
332  PMatI.UseExternalData(elfun.GetData(), dof, dim);
333  elmat.SetSize(dof*dim);
334 
335  int intorder = 2*el.GetOrder() + 3; // <---
336  const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
337 
338  elmat = 0.0;
339  model->SetTransformation(Tr);
340  for (int i = 0; i < ir.GetNPoints(); i++)
341  {
342  const IntegrationPoint &ip = ir.IntPoint(i);
343  Tr.SetIntPoint(&ip);
344  CalcInverse(Tr.Jacobian(), J0i);
345 
346  el.CalcDShape(ip, DSh);
347  Mult(DSh, J0i, DS);
348  MultAtB(PMatI, DS, J);
349 
350  model->AssembleH(J, DS, ip.weight*Tr.Weight(), elmat);
351  }
352 }
353 
355 {
356  PMatI.ClearExternalData();
357  PMatO.ClearExternalData();
358 }
359 
360 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
Abstract class for Finite Elements.
Definition: fe.hpp:46
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
Definition: densemat.cpp:3216
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P=P(J).
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W=W(J).
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:115
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:310
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Definition: nonlininteg.cpp:60
double Det() const
Calculates the determinant of the matrix (for 2x2 or 3x3 matrices)
Definition: densemat.cpp:416
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:51
int GetOrder() const
Returns the order of the finite element.
Definition: fe.hpp:124
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P=P(J).
Definition: nonlininteg.cpp:41
const IntegrationPoint & GetIntPoint()
Definition: eltrans.hpp:53
double * GetData() const
Definition: vector.hpp:114
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2974
int dim
Definition: ex3.cpp:47
void EvalCoeffs() const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:470
const DenseMatrix & Jacobian()
Definition: eltrans.hpp:64
void SetTransformation(ElementTransformation &_T)
An element transformation that can be used to evaluate coefficients.
Definition: nonlininteg.hpp:58
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3010
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:118
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:17
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:392
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:3339
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:121
void mfem_error(const char *msg)
Definition: error.cpp:106
void ClearExternalData()
Definition: densemat.hpp:72
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:3158
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const =0
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Class for integration point with weight.
Definition: intrules.hpp:25
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W=W(J).
Definition: nonlininteg.cpp:34
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P=P(J).
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:3491
Vector data type.
Definition: vector.hpp:36
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
Definition: nonlininteg.cpp:25
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:60
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:82
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:343
virtual double EvalW(const DenseMatrix &J) const =0
Evaluate the strain energy density function, W=W(J).
ElementTransformation * T
Definition: nonlininteg.hpp:52