MFEM  v3.4
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 
17 void NonlinearFormIntegrator::AssembleElementVector(
18  const FiniteElement &el, ElementTransformation &Tr,
19  const Vector &elfun, Vector &elvect)
20 {
21  mfem_error("NonlinearFormIntegrator::AssembleElementVector"
22  " is not overloaded!");
23 }
24 
25 void NonlinearFormIntegrator::AssembleFaceVector(
26  const FiniteElement &el1, const FiniteElement &el2,
27  FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
28 {
29  mfem_error("NonlinearFormIntegrator::AssembleFaceVector"
30  " is not overloaded!");
31 }
32 
33 void NonlinearFormIntegrator::AssembleElementGrad(
34  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
35  DenseMatrix &elmat)
36 {
37  mfem_error("NonlinearFormIntegrator::AssembleElementGrad"
38  " is not overloaded!");
39 }
40 
41 void NonlinearFormIntegrator::AssembleFaceGrad(
42  const FiniteElement &el1, const FiniteElement &el2,
43  FaceElementTransformations &Tr, const Vector &elfun,
44  DenseMatrix &elmat)
45 {
46  mfem_error("NonlinearFormIntegrator::AssembleFaceGrad"
47  " is not overloaded!");
48 }
49 
50 double NonlinearFormIntegrator::GetElementEnergy(
51  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
52 {
53  mfem_error("NonlinearFormIntegrator::GetElementEnergy"
54  " is not overloaded!");
55  return 0.0;
56 }
57 
58 
59 void BlockNonlinearFormIntegrator::AssembleElementVector(
60  const Array<const FiniteElement *> &el,
61  ElementTransformation &Tr,
62  const Array<const Vector *> &elfun,
63  const Array<Vector *> &elvec)
64 {
65  mfem_error("BlockNonlinearFormIntegrator::AssembleElementVector"
66  " is not overloaded!");
67 }
68 
69 void BlockNonlinearFormIntegrator::AssembleFaceVector(
70  const Array<const FiniteElement *> &el1,
71  const Array<const FiniteElement *> &el2,
72  FaceElementTransformations &Tr,
73  const Array<const Vector *> &elfun,
74  const Array<Vector *> &elvect)
75 {
76  mfem_error("BlockNonlinearFormIntegrator::AssembleFaceVector"
77  " is not overloaded!");
78 }
79 
80 void BlockNonlinearFormIntegrator::AssembleElementGrad(
81  const Array<const FiniteElement*> &el,
82  ElementTransformation &Tr,
83  const Array<const Vector *> &elfun,
84  const Array2D<DenseMatrix *> &elmats)
85 {
86  mfem_error("BlockNonlinearFormIntegrator::AssembleElementGrad"
87  " is not overloaded!");
88 }
89 
90 void BlockNonlinearFormIntegrator::AssembleFaceGrad(
91  const Array<const FiniteElement *>&el1,
92  const Array<const FiniteElement *>&el2,
93  FaceElementTransformations &Tr,
94  const Array<const Vector *> &elfun,
95  const Array2D<DenseMatrix *> &elmats)
96 {
97  mfem_error("BlockNonlinearFormIntegrator::AssembleFaceGrad"
98  " is not overloaded!");
99 }
100 
101 double BlockNonlinearFormIntegrator::GetElementEnergy(
102  const Array<const FiniteElement *>&el,
103  ElementTransformation &Tr,
104  const Array<const Vector *>&elfun)
105 {
106  mfem_error("BlockNonlinearFormIntegrator::GetElementEnergy"
107  " is not overloaded!");
108  return 0.0;
109 }
110 
111 
112 double InverseHarmonicModel::EvalW(const DenseMatrix &J) const
113 {
114  Z.SetSize(J.Width());
115  CalcAdjugateTranspose(J, Z);
116  return 0.5*(Z*Z)/J.Det();
117 }
118 
119 void InverseHarmonicModel::EvalP(const DenseMatrix &J, DenseMatrix &P) const
120 {
121  int dim = J.Width();
122  double t;
123 
124  Z.SetSize(dim);
125  S.SetSize(dim);
126  CalcAdjugateTranspose(J, Z);
127  MultAAt(Z, S);
128  t = 0.5*S.Trace();
129  for (int i = 0; i < dim; i++)
130  {
131  S(i,i) -= t;
132  }
133  t = J.Det();
134  S *= -1.0/(t*t);
135  Mult(S, Z, P);
136 }
137 
138 void InverseHarmonicModel::AssembleH(
139  const DenseMatrix &J, const DenseMatrix &DS, const double weight,
140  DenseMatrix &A) const
141 {
142  int dof = DS.Height(), dim = DS.Width();
143  double t;
144 
145  Z.SetSize(dim);
146  S.SetSize(dim);
147  G.SetSize(dof, dim);
148  C.SetSize(dof, dim);
149 
150  CalcAdjugateTranspose(J, Z);
151  MultAAt(Z, S);
152 
153  t = 1.0/J.Det();
154  Z *= t; // Z = J^{-t}
155  S *= t; // S = |J| (J.J^t)^{-1}
156  t = 0.5*S.Trace();
157 
158  MultABt(DS, Z, G); // G = DS.J^{-1}
159  Mult(G, S, C);
160 
161  // 1.
162  for (int i = 0; i < dof; i++)
163  for (int j = 0; j <= i; j++)
164  {
165  double a = 0.0;
166  for (int d = 0; d < dim; d++)
167  {
168  a += G(i,d)*G(j,d);
169  }
170  a *= weight;
171  for (int k = 0; k < dim; k++)
172  for (int l = 0; l <= k; l++)
173  {
174  double b = a*S(k,l);
175  A(i+k*dof,j+l*dof) += b;
176  if (i != j)
177  {
178  A(j+k*dof,i+l*dof) += b;
179  }
180  if (k != l)
181  {
182  A(i+l*dof,j+k*dof) += b;
183  if (i != j)
184  {
185  A(j+l*dof,i+k*dof) += b;
186  }
187  }
188  }
189  }
190 
191  // 2.
192  for (int i = 1; i < dof; i++)
193  for (int j = 0; j < i; j++)
194  {
195  for (int k = 1; k < dim; k++)
196  for (int l = 0; l < k; l++)
197  {
198  double a =
199  weight*(C(i,l)*G(j,k) - C(i,k)*G(j,l) +
200  C(j,k)*G(i,l) - C(j,l)*G(i,k) +
201  t*(G(i,k)*G(j,l) - G(i,l)*G(j,k)));
202 
203  A(i+k*dof,j+l*dof) += a;
204  A(j+l*dof,i+k*dof) += a;
205 
206  A(i+l*dof,j+k*dof) -= a;
207  A(j+k*dof,i+l*dof) -= a;
208  }
209  }
210 }
211 
212 
213 inline void NeoHookeanModel::EvalCoeffs() const
214 {
215  mu = c_mu->Eval(*Ttr, Ttr->GetIntPoint());
216  K = c_K->Eval(*Ttr, Ttr->GetIntPoint());
217  if (c_g)
218  {
219  g = c_g->Eval(*Ttr, Ttr->GetIntPoint());
220  }
221 }
222 
223 double NeoHookeanModel::EvalW(const DenseMatrix &J) const
224 {
225  int dim = J.Width();
226 
227  if (have_coeffs)
228  {
229  EvalCoeffs();
230  }
231 
232  double dJ = J.Det();
233  double sJ = dJ/g;
234  double bI1 = pow(dJ, -2.0/dim)*(J*J); // \bar{I}_1
235 
236  return 0.5*(mu*(bI1 - dim) + K*(sJ - 1.0)*(sJ - 1.0));
237 }
238 
239 void NeoHookeanModel::EvalP(const DenseMatrix &J, DenseMatrix &P) const
240 {
241  int dim = J.Width();
242 
243  if (have_coeffs)
244  {
245  EvalCoeffs();
246  }
247 
248  Z.SetSize(dim);
249  CalcAdjugateTranspose(J, Z);
250 
251  double dJ = J.Det();
252  double a = mu*pow(dJ, -2.0/dim);
253  double b = K*(dJ/g - 1.0)/g - a*(J*J)/(dim*dJ);
254 
255  P = 0.0;
256  P.Add(a, J);
257  P.Add(b, Z);
258 }
259 
260 void NeoHookeanModel::AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
261  const double weight, DenseMatrix &A) const
262 {
263  int dof = DS.Height(), dim = DS.Width();
264 
265  if (have_coeffs)
266  {
267  EvalCoeffs();
268  }
269 
270  Z.SetSize(dim);
271  G.SetSize(dof, dim);
272  C.SetSize(dof, dim);
273 
274  double dJ = J.Det();
275  double sJ = dJ/g;
276  double a = mu*pow(dJ, -2.0/dim);
277  double bc = a*(J*J)/dim;
278  double b = bc - K*sJ*(sJ - 1.0);
279  double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0);
280 
281  CalcAdjugateTranspose(J, Z);
282  Z *= (1.0/dJ); // Z = J^{-t}
283 
284  MultABt(DS, J, C); // C = DS J^t
285  MultABt(DS, Z, G); // G = DS J^{-1}
286 
287  a *= weight;
288  b *= weight;
289  c *= weight;
290 
291  // 1.
292  for (int i = 0; i < dof; i++)
293  for (int k = 0; k <= i; k++)
294  {
295  double s = 0.0;
296  for (int d = 0; d < dim; d++)
297  {
298  s += DS(i,d)*DS(k,d);
299  }
300  s *= a;
301 
302  for (int d = 0; d < dim; d++)
303  {
304  A(i+d*dof,k+d*dof) += s;
305  }
306 
307  if (k != i)
308  for (int d = 0; d < dim; d++)
309  {
310  A(k+d*dof,i+d*dof) += s;
311  }
312  }
313 
314  a *= (-2.0/dim);
315 
316  // 2.
317  for (int i = 0; i < dof; i++)
318  for (int j = 0; j < dim; j++)
319  for (int k = 0; k < dof; k++)
320  for (int l = 0; l < dim; l++)
321  {
322  A(i+j*dof,k+l*dof) +=
323  a*(C(i,j)*G(k,l) + G(i,j)*C(k,l)) +
324  b*G(i,l)*G(k,j) + c*G(i,j)*G(k,l);
325  }
326 }
327 
328 
329 double HyperelasticNLFIntegrator::GetElementEnergy(const FiniteElement &el,
330  ElementTransformation &Ttr,
331  const Vector &elfun)
332 {
333  int dof = el.GetDof(), dim = el.GetDim();
334  double energy;
335 
336  DSh.SetSize(dof, dim);
337  Jrt.SetSize(dim);
338  Jpr.SetSize(dim);
339  Jpt.SetSize(dim);
340  PMatI.UseExternalData(elfun.GetData(), dof, dim);
341 
342  const IntegrationRule *ir = IntRule;
343  if (!ir)
344  {
345  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
346  }
347 
348  energy = 0.0;
349  model->SetTransformation(Ttr);
350  for (int i = 0; i < ir->GetNPoints(); i++)
351  {
352  const IntegrationPoint &ip = ir->IntPoint(i);
353  Ttr.SetIntPoint(&ip);
354  CalcInverse(Ttr.Jacobian(), Jrt);
355 
356  el.CalcDShape(ip, DSh);
357  MultAtB(PMatI, DSh, Jpr);
358  Mult(Jpr, Jrt, Jpt);
359 
360  energy += ip.weight * Ttr.Weight() * model->EvalW(Jpt);
361  }
362 
363  return energy;
364 }
365 
366 void HyperelasticNLFIntegrator::AssembleElementVector(
367  const FiniteElement &el, ElementTransformation &Ttr,
368  const Vector &elfun, Vector &elvect)
369 {
370  int dof = el.GetDof(), dim = el.GetDim();
371 
372  DSh.SetSize(dof, dim);
373  DS.SetSize(dof, dim);
374  Jrt.SetSize(dim);
375  Jpt.SetSize(dim);
376  P.SetSize(dim);
377  PMatI.UseExternalData(elfun.GetData(), dof, dim);
378  elvect.SetSize(dof*dim);
379  PMatO.UseExternalData(elvect.GetData(), dof, dim);
380 
381  const IntegrationRule *ir = IntRule;
382  if (!ir)
383  {
384  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
385  }
386 
387  elvect = 0.0;
388  model->SetTransformation(Ttr);
389  for (int i = 0; i < ir->GetNPoints(); i++)
390  {
391  const IntegrationPoint &ip = ir->IntPoint(i);
392  Ttr.SetIntPoint(&ip);
393  CalcInverse(Ttr.Jacobian(), Jrt);
394 
395  el.CalcDShape(ip, DSh);
396  Mult(DSh, Jrt, DS);
397  MultAtB(PMatI, DS, Jpt);
398 
399  model->EvalP(Jpt, P);
400 
401  P *= ip.weight * Ttr.Weight();
402  AddMultABt(DS, P, PMatO);
403  }
404 }
405 
406 void HyperelasticNLFIntegrator::AssembleElementGrad(const FiniteElement &el,
407  ElementTransformation &Ttr,
408  const Vector &elfun,
409  DenseMatrix &elmat)
410 {
411  int dof = el.GetDof(), dim = el.GetDim();
412 
413  DSh.SetSize(dof, dim);
414  DS.SetSize(dof, dim);
415  Jrt.SetSize(dim);
416  Jpt.SetSize(dim);
417  PMatI.UseExternalData(elfun.GetData(), dof, dim);
418  elmat.SetSize(dof*dim);
419 
420  const IntegrationRule *ir = IntRule;
421  if (!ir)
422  {
423  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
424  }
425 
426  elmat = 0.0;
427  model->SetTransformation(Ttr);
428  for (int i = 0; i < ir->GetNPoints(); i++)
429  {
430  const IntegrationPoint &ip = ir->IntPoint(i);
431  Ttr.SetIntPoint(&ip);
432  CalcInverse(Ttr.Jacobian(), Jrt);
433 
434  el.CalcDShape(ip, DSh);
435  Mult(DSh, Jrt, DS);
436  MultAtB(PMatI, DS, Jpt);
437 
438  model->AssembleH(Jpt, DS, ip.weight * Ttr.Weight(), elmat);
439  }
440 }
441 
442 double IncompressibleNeoHookeanIntegrator::GetElementEnergy(
443  const Array<const FiniteElement *>&el,
444  ElementTransformation &Tr,
445  const Array<const Vector *>&elfun)
446 {
447  if (el.Size() != 2)
448  {
449  mfem_error("IncompressibleNeoHookeanIntegrator::GetElementEnergy"
450  " has incorrect block finite element space size!");
451  }
452 
453  int dof_u = el[0]->GetDof();
454  int dim = el[0]->GetDim();
455 
456  DSh_u.SetSize(dof_u, dim);
457  J0i.SetSize(dim);
458  J1.SetSize(dim);
459  J.SetSize(dim);
460  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
461 
462  int intorder = 2*el[0]->GetOrder() + 3; // <---
463  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
464 
465  double energy = 0.0;
466  double mu = 0.0;
467 
468  for (int i = 0; i < ir.GetNPoints(); ++i)
469  {
470  const IntegrationPoint &ip = ir.IntPoint(i);
471  Tr.SetIntPoint(&ip);
472  CalcInverse(Tr.Jacobian(), J0i);
473 
474  el[0]->CalcDShape(ip, DSh_u);
475  MultAtB(PMatI_u, DSh_u, J1);
476  Mult(J1, J0i, J);
477 
478  mu = c_mu->Eval(Tr, ip);
479 
480  energy += ip.weight*Tr.Weight()*(mu/2.0)*(J*J - 3);
481  }
482 
483  return energy;
484 }
485 
486 void IncompressibleNeoHookeanIntegrator::AssembleElementVector(
487  const Array<const FiniteElement *> &el,
488  ElementTransformation &Tr,
489  const Array<const Vector *> &elfun,
490  const Array<Vector *> &elvec)
491 {
492  if (el.Size() != 2)
493  {
494  mfem_error("IncompressibleNeoHookeanIntegrator::AssembleElementVector"
495  " has finite element space of incorrect block number");
496  }
497 
498  int dof_u = el[0]->GetDof();
499  int dof_p = el[1]->GetDof();
500 
501  int dim = el[0]->GetDim();
502  int spaceDim = Tr.GetSpaceDim();
503 
504  if (dim != spaceDim)
505  {
506  mfem_error("IncompressibleNeoHookeanIntegrator::AssembleElementVector"
507  " is not defined on manifold meshes");
508  }
509 
510 
511  DSh_u.SetSize(dof_u, dim);
512  DS_u.SetSize(dof_u, dim);
513  J0i.SetSize(dim);
514  F.SetSize(dim);
515  FinvT.SetSize(dim);
516  P.SetSize(dim);
517  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
518  elvec[0]->SetSize(dof_u*dim);
519  PMatO_u.UseExternalData(elvec[0]->GetData(), dof_u, dim);
520 
521  Sh_p.SetSize(dof_p);
522  elvec[1]->SetSize(dof_p);
523 
524  int intorder = 2*el[0]->GetOrder() + 3; // <---
525  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
526 
527  *elvec[0] = 0.0;
528  *elvec[1] = 0.0;
529 
530  for (int i = 0; i < ir.GetNPoints(); ++i)
531  {
532  const IntegrationPoint &ip = ir.IntPoint(i);
533  Tr.SetIntPoint(&ip);
534  CalcInverse(Tr.Jacobian(), J0i);
535 
536  el[0]->CalcDShape(ip, DSh_u);
537  Mult(DSh_u, J0i, DS_u);
538  MultAtB(PMatI_u, DS_u, F);
539 
540  el[1]->CalcShape(ip, Sh_p);
541 
542  double pres = Sh_p * *elfun[1];
543  double mu = c_mu->Eval(Tr, ip);
544  double dJ = F.Det();
545 
546  CalcInverseTranspose(F, FinvT);
547 
548  P = 0.0;
549  P.Add(mu * dJ, F);
550  P.Add(-1.0 * pres * dJ, FinvT);
551  P *= ip.weight*Tr.Weight();
552 
553  AddMultABt(DS_u, P, PMatO_u);
554 
555  elvec[1]->Add(ip.weight * Tr.Weight() * (dJ - 1.0), Sh_p);
556  }
557 
558 }
559 
560 void IncompressibleNeoHookeanIntegrator::AssembleElementGrad(
561  const Array<const FiniteElement*> &el,
562  ElementTransformation &Tr,
563  const Array<const Vector *> &elfun,
564  const Array2D<DenseMatrix *> &elmats)
565 {
566  int dof_u = el[0]->GetDof();
567  int dof_p = el[1]->GetDof();
568 
569  int dim = el[0]->GetDim();
570 
571  elmats(0,0)->SetSize(dof_u*dim, dof_u*dim);
572  elmats(0,1)->SetSize(dof_u*dim, dof_p);
573  elmats(1,0)->SetSize(dof_p, dof_u*dim);
574  elmats(1,1)->SetSize(dof_p, dof_p);
575 
576  *elmats(0,0) = 0.0;
577  *elmats(0,1) = 0.0;
578  *elmats(1,0) = 0.0;
579  *elmats(1,1) = 0.0;
580 
581  DSh_u.SetSize(dof_u, dim);
582  DS_u.SetSize(dof_u, dim);
583  J0i.SetSize(dim);
584  F.SetSize(dim);
585  FinvT.SetSize(dim);
586  Finv.SetSize(dim);
587  P.SetSize(dim);
588  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
589  Sh_p.SetSize(dof_p);
590 
591  int intorder = 2*el[0]->GetOrder() + 3; // <---
592  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
593 
594  for (int i = 0; i < ir.GetNPoints(); ++i)
595  {
596  const IntegrationPoint &ip = ir.IntPoint(i);
597  Tr.SetIntPoint(&ip);
598  CalcInverse(Tr.Jacobian(), J0i);
599 
600  el[0]->CalcDShape(ip, DSh_u);
601  Mult(DSh_u, J0i, DS_u);
602  MultAtB(PMatI_u, DS_u, F);
603 
604  el[1]->CalcShape(ip, Sh_p);
605  double pres = Sh_p * *elfun[1];
606  double mu = c_mu->Eval(Tr, ip);
607  double dJ = F.Det();
608  double dJ_FinvT_DS;
609 
610  CalcInverseTranspose(F, FinvT);
611 
612  // u,u block
613  for (int i_u = 0; i_u < dof_u; ++i_u)
614  {
615  for (int i_dim = 0; i_dim < dim; ++i_dim)
616  {
617  for (int j_u = 0; j_u < dof_u; ++j_u)
618  {
619  for (int j_dim = 0; j_dim < dim; ++j_dim)
620  {
621 
622  // m = j_dim;
623  // k = i_dim;
624 
625  for (int n=0; n<dim; ++n)
626  {
627  for (int l=0; l<dim; ++l)
628  {
629  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
630  dJ * (mu * F(i_dim, l) - pres * FinvT(i_dim,l)) *
631  FinvT(j_dim,n) * DS_u(i_u,l) * DS_u(j_u, n) *
632  ip.weight * Tr.Weight();
633 
634  if (j_dim == i_dim && n==l)
635  {
636  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
637  dJ * mu * DS_u(i_u, l) * DS_u(j_u,n) *
638  ip.weight * Tr.Weight();
639  }
640 
641  // a = n;
642  // b = m;
643  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
644  dJ * pres * FinvT(i_dim, n) *
645  FinvT(j_dim,l) * DS_u(i_u,l) * DS_u(j_u,n) *
646  ip.weight * Tr.Weight();
647  }
648  }
649  }
650  }
651  }
652  }
653 
654  // u,p and p,u blocks
655  for (int i_p = 0; i_p < dof_p; ++i_p)
656  {
657  for (int j_u = 0; j_u < dof_u; ++j_u)
658  {
659  for (int dim_u = 0; dim_u < dim; ++dim_u)
660  {
661  for (int l=0; l<dim; ++l)
662  {
663  dJ_FinvT_DS = dJ * FinvT(dim_u,l) * DS_u(j_u, l) * Sh_p(i_p) *
664  ip.weight * Tr.Weight();
665  (*elmats(1,0))(i_p, j_u + dof_u * dim_u) += dJ_FinvT_DS;
666  (*elmats(0,1))(j_u + dof_u * dim_u, i_p) -= dJ_FinvT_DS;
667 
668  }
669  }
670  }
671  }
672  }
673 
674 }
675 
676 }
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:232
mfem::FiniteElement
Abstract class for Finite Elements.
Definition: fe.hpp:140
mfem::MultABt
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:3356
mfem::Array::Size
int Size() const
Logical size of the array.
Definition: array.hpp:133
mfem::NeoHookeanModel::EvalP
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: nonlininteg.cpp:239
mfem::NeoHookeanModel::EvalW
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: nonlininteg.cpp:223
mfem::FiniteElement::GetDim
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:211
mfem::IncompressibleNeoHookeanIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:560
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
mfem::BlockNonlinearFormIntegrator::AssembleFaceVector
virtual void AssembleFaceVector(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvect)
Definition: nonlininteg.cpp:69
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
mfem::InverseHarmonicModel::AssembleH
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: nonlininteg.cpp:138
mfem::DenseMatrix::Det
double Det() const
Definition: densemat.cpp:460
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
mfem::FiniteElement::GetOrder
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:221
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
mfem::HyperelasticModel::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
mfem::InverseHarmonicModel::EvalP
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: nonlininteg.cpp:119
mfem::ElementTransformation::GetIntPoint
const IntegrationPoint & GetIntPoint()
Definition: eltrans.hpp:54
mfem::Vector::GetData
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
mfem::BlockNonlinearFormIntegrator::AssembleFaceGrad
virtual void AssembleFaceGrad(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Definition: nonlininteg.cpp:90
mfem::HyperelasticModel::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
mfem::BlockNonlinearFormIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:80
mfem::Array2D::SetSize
void SetSize(int m, int n)
Definition: array.hpp:308
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:235
mfem::CalcAdjugateTranspose
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:3114
dim
int dim
Definition: ex3.cpp:47
mfem::NeoHookeanModel::EvalCoeffs
void EvalCoeffs() const
Definition: nonlininteg.cpp:213
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
mfem::DenseMatrix::Add
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:553
mfem::mfem_error
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:146
mfem::ElementTransformation::Jacobian
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:67
mfem::ElementTransformation::GetSpaceDim
int GetSpaceDim() const
Get the dimension of the target (physical) space.
Definition: eltrans.hpp:93
mfem::NonlinearFormIntegrator::AssembleFaceGrad
virtual void AssembleFaceGrad(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local action of the gradient of the NonlinearFormIntegrator resulting from a face integr...
Definition: nonlininteg.cpp:41
mfem::IncompressibleNeoHookeanIntegrator::GetElementEnergy
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
Definition: nonlininteg.cpp:442
mfem::HyperelasticNLFIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:406
mfem::NonlinearFormIntegrator::AssembleElementVector
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: nonlininteg.cpp:17
mfem::HyperelasticNLFIntegrator::GetElementEnergy
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun)
Computes the integral of W(Jacobian(Trt)) over a target zone.
Definition: nonlininteg.cpp:329
mfem::Array2D
Dynamic 2D array using row-major layout.
Definition: array.hpp:289
mfem::HyperelasticModel::SetTransformation
void SetTransformation(ElementTransformation &_Ttr)
Definition: nonlininteg.hpp:128
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3150
mfem::FiniteElement::GetGeomType
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:214
mfem::NonlinearFormIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:33
mfem::DenseMatrix::Trace
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:436
mfem::AddMultABt
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:3479
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:217
mfem::Array::SetSize
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
mfem::MultAAt
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:3298
mfem::NonlinearFormIntegrator::AssembleFaceVector
virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator resulting from a face integral term...
Definition: nonlininteg.cpp:25
mfem::CalcInverseTranspose
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:3233
mfem::HyperelasticModel::Ttr
ElementTransformation * Ttr
Definition: nonlininteg.hpp:117
mfem::NonlinearFormIntegrator::IntRule
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:29
mfem::NeoHookeanModel::AssembleH
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: nonlininteg.cpp:260
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::InverseHarmonicModel::EvalW
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: nonlininteg.cpp:112
mfem::HyperelasticModel::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
mfem::IncompressibleNeoHookeanIntegrator::AssembleElementVector
virtual void AssembleElementVector(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvec)
Perform the local action of the NonlinearFormIntegrator.
Definition: nonlininteg.cpp:486
mfem::Coefficient::Eval
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
mfem::HyperelasticNLFIntegrator::AssembleElementVector
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: nonlininteg.cpp:366
mfem::BlockNonlinearFormIntegrator::GetElementEnergy
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
Definition: nonlininteg.cpp:101
mfem::MultAtB
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:3631
mfem::Vector
Vector data type.
Definition: vector.hpp:48
mfem::NonlinearFormIntegrator::GetElementEnergy
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
Definition: nonlininteg.cpp:50
mfem::DenseMatrix::UseExternalData
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:64
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:86
mfem::FiniteElement::CalcDShape
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...
mfem::BlockNonlinearFormIntegrator::AssembleElementVector
virtual void AssembleElementVector(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvec)
Perform the local action of the BlockNonlinearFormIntegrator.
Definition: nonlininteg.cpp:59
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:353
Generated on Fri Feb 18 2022 19:03:35 for MFEM by   doxygen 1.8.5