MFEM  v4.2.0
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-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "fem.hpp"
13 #include "../general/forall.hpp"
14 
15 namespace mfem
16 {
17 
19 {
20  mfem_error ("NonlinearFormIntegrator::AssemblePA(...)\n"
21  " is not implemented for this class.");
22 }
23 
25  const FiniteElementSpace &)
26 {
27  mfem_error ("NonlinearFormIntegrator::AssemblePA(...)\n"
28  " is not implemented for this class.");
29 }
30 
32 {
33  mfem_error ("NonlinearFormIntegrator::AddMultPA(...)\n"
34  " is not implemented for this class.");
35 }
36 
38  const FiniteElement &el, ElementTransformation &Tr,
39  const Vector &elfun, Vector &elvect)
40 {
41  mfem_error("NonlinearFormIntegrator::AssembleElementVector"
42  " is not overloaded!");
43 }
44 
46  const FiniteElement &el1, const FiniteElement &el2,
47  FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
48 {
49  mfem_error("NonlinearFormIntegrator::AssembleFaceVector"
50  " is not overloaded!");
51 }
52 
54  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
55  DenseMatrix &elmat)
56 {
57  mfem_error("NonlinearFormIntegrator::AssembleElementGrad"
58  " is not overloaded!");
59 }
60 
62  const FiniteElement &el1, const FiniteElement &el2,
63  FaceElementTransformations &Tr, const Vector &elfun,
64  DenseMatrix &elmat)
65 {
66  mfem_error("NonlinearFormIntegrator::AssembleFaceGrad"
67  " is not overloaded!");
68 }
69 
71  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
72 {
73  mfem_error("NonlinearFormIntegrator::GetElementEnergy"
74  " is not overloaded!");
75  return 0.0;
76 }
77 
78 
82  const Array<const Vector *> &elfun,
83  const Array<Vector *> &elvec)
84 {
85  mfem_error("BlockNonlinearFormIntegrator::AssembleElementVector"
86  " is not overloaded!");
87 }
88 
93  const Array<const Vector *> &elfun,
94  const Array<Vector *> &elvect)
95 {
96  mfem_error("BlockNonlinearFormIntegrator::AssembleFaceVector"
97  " is not overloaded!");
98 }
99 
101  const Array<const FiniteElement*> &el,
103  const Array<const Vector *> &elfun,
104  const Array2D<DenseMatrix *> &elmats)
105 {
106  mfem_error("BlockNonlinearFormIntegrator::AssembleElementGrad"
107  " is not overloaded!");
108 }
109 
114  const Array<const Vector *> &elfun,
115  const Array2D<DenseMatrix *> &elmats)
116 {
117  mfem_error("BlockNonlinearFormIntegrator::AssembleFaceGrad"
118  " is not overloaded!");
119 }
120 
124  const Array<const Vector *>&elfun)
125 {
126  mfem_error("BlockNonlinearFormIntegrator::GetElementEnergy"
127  " is not overloaded!");
128  return 0.0;
129 }
130 
131 
133 {
134  Z.SetSize(J.Width());
136  return 0.5*(Z*Z)/J.Det();
137 }
138 
140 {
141  int dim = J.Width();
142  double t;
143 
144  Z.SetSize(dim);
145  S.SetSize(dim);
147  MultAAt(Z, S);
148  t = 0.5*S.Trace();
149  for (int i = 0; i < dim; i++)
150  {
151  S(i,i) -= t;
152  }
153  t = J.Det();
154  S *= -1.0/(t*t);
155  Mult(S, Z, P);
156 }
157 
159  const DenseMatrix &J, const DenseMatrix &DS, const double weight,
160  DenseMatrix &A) const
161 {
162  int dof = DS.Height(), dim = DS.Width();
163  double t;
164 
165  Z.SetSize(dim);
166  S.SetSize(dim);
167  G.SetSize(dof, dim);
168  C.SetSize(dof, dim);
169 
171  MultAAt(Z, S);
172 
173  t = 1.0/J.Det();
174  Z *= t; // Z = J^{-t}
175  S *= t; // S = |J| (J.J^t)^{-1}
176  t = 0.5*S.Trace();
177 
178  MultABt(DS, Z, G); // G = DS.J^{-1}
179  Mult(G, S, C);
180 
181  // 1.
182  for (int i = 0; i < dof; i++)
183  for (int j = 0; j <= i; j++)
184  {
185  double a = 0.0;
186  for (int d = 0; d < dim; d++)
187  {
188  a += G(i,d)*G(j,d);
189  }
190  a *= weight;
191  for (int k = 0; k < dim; k++)
192  for (int l = 0; l <= k; l++)
193  {
194  double b = a*S(k,l);
195  A(i+k*dof,j+l*dof) += b;
196  if (i != j)
197  {
198  A(j+k*dof,i+l*dof) += b;
199  }
200  if (k != l)
201  {
202  A(i+l*dof,j+k*dof) += b;
203  if (i != j)
204  {
205  A(j+l*dof,i+k*dof) += b;
206  }
207  }
208  }
209  }
210 
211  // 2.
212  for (int i = 1; i < dof; i++)
213  for (int j = 0; j < i; j++)
214  {
215  for (int k = 1; k < dim; k++)
216  for (int l = 0; l < k; l++)
217  {
218  double a =
219  weight*(C(i,l)*G(j,k) - C(i,k)*G(j,l) +
220  C(j,k)*G(i,l) - C(j,l)*G(i,k) +
221  t*(G(i,k)*G(j,l) - G(i,l)*G(j,k)));
222 
223  A(i+k*dof,j+l*dof) += a;
224  A(j+l*dof,i+k*dof) += a;
225 
226  A(i+l*dof,j+k*dof) -= a;
227  A(j+k*dof,i+l*dof) -= a;
228  }
229  }
230 }
231 
232 
233 inline void NeoHookeanModel::EvalCoeffs() const
234 {
235  mu = c_mu->Eval(*Ttr, Ttr->GetIntPoint());
236  K = c_K->Eval(*Ttr, Ttr->GetIntPoint());
237  if (c_g)
238  {
239  g = c_g->Eval(*Ttr, Ttr->GetIntPoint());
240  }
241 }
242 
243 double NeoHookeanModel::EvalW(const DenseMatrix &J) const
244 {
245  int dim = J.Width();
246 
247  if (have_coeffs)
248  {
249  EvalCoeffs();
250  }
251 
252  double dJ = J.Det();
253  double sJ = dJ/g;
254  double bI1 = pow(dJ, -2.0/dim)*(J*J); // \bar{I}_1
255 
256  return 0.5*(mu*(bI1 - dim) + K*(sJ - 1.0)*(sJ - 1.0));
257 }
258 
260 {
261  int dim = J.Width();
262 
263  if (have_coeffs)
264  {
265  EvalCoeffs();
266  }
267 
268  Z.SetSize(dim);
270 
271  double dJ = J.Det();
272  double a = mu*pow(dJ, -2.0/dim);
273  double b = K*(dJ/g - 1.0)/g - a*(J*J)/(dim*dJ);
274 
275  P = 0.0;
276  P.Add(a, J);
277  P.Add(b, Z);
278 }
279 
281  const double weight, DenseMatrix &A) const
282 {
283  int dof = DS.Height(), dim = DS.Width();
284 
285  if (have_coeffs)
286  {
287  EvalCoeffs();
288  }
289 
290  Z.SetSize(dim);
291  G.SetSize(dof, dim);
292  C.SetSize(dof, dim);
293 
294  double dJ = J.Det();
295  double sJ = dJ/g;
296  double a = mu*pow(dJ, -2.0/dim);
297  double bc = a*(J*J)/dim;
298  double b = bc - K*sJ*(sJ - 1.0);
299  double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0);
300 
302  Z *= (1.0/dJ); // Z = J^{-t}
303 
304  MultABt(DS, J, C); // C = DS J^t
305  MultABt(DS, Z, G); // G = DS J^{-1}
306 
307  a *= weight;
308  b *= weight;
309  c *= weight;
310 
311  // 1.
312  for (int i = 0; i < dof; i++)
313  for (int k = 0; k <= i; k++)
314  {
315  double s = 0.0;
316  for (int d = 0; d < dim; d++)
317  {
318  s += DS(i,d)*DS(k,d);
319  }
320  s *= a;
321 
322  for (int d = 0; d < dim; d++)
323  {
324  A(i+d*dof,k+d*dof) += s;
325  }
326 
327  if (k != i)
328  for (int d = 0; d < dim; d++)
329  {
330  A(k+d*dof,i+d*dof) += s;
331  }
332  }
333 
334  a *= (-2.0/dim);
335 
336  // 2.
337  for (int i = 0; i < dof; i++)
338  for (int j = 0; j < dim; j++)
339  for (int k = 0; k < dof; k++)
340  for (int l = 0; l < dim; l++)
341  {
342  A(i+j*dof,k+l*dof) +=
343  a*(C(i,j)*G(k,l) + G(i,j)*C(k,l)) +
344  b*G(i,l)*G(k,j) + c*G(i,j)*G(k,l);
345  }
346 }
347 
348 
351  const Vector &elfun)
352 {
353  int dof = el.GetDof(), dim = el.GetDim();
354  double energy;
355 
356  DSh.SetSize(dof, dim);
357  Jrt.SetSize(dim);
358  Jpr.SetSize(dim);
359  Jpt.SetSize(dim);
360  PMatI.UseExternalData(elfun.GetData(), dof, dim);
361 
362  const IntegrationRule *ir = IntRule;
363  if (!ir)
364  {
365  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
366  }
367 
368  energy = 0.0;
369  model->SetTransformation(Ttr);
370  for (int i = 0; i < ir->GetNPoints(); i++)
371  {
372  const IntegrationPoint &ip = ir->IntPoint(i);
373  Ttr.SetIntPoint(&ip);
374  CalcInverse(Ttr.Jacobian(), Jrt);
375 
376  el.CalcDShape(ip, DSh);
377  MultAtB(PMatI, DSh, Jpr);
378  Mult(Jpr, Jrt, Jpt);
379 
380  energy += ip.weight * Ttr.Weight() * model->EvalW(Jpt);
381  }
382 
383  return energy;
384 }
385 
387  const FiniteElement &el, ElementTransformation &Ttr,
388  const Vector &elfun, Vector &elvect)
389 {
390  int dof = el.GetDof(), dim = el.GetDim();
391 
392  DSh.SetSize(dof, dim);
393  DS.SetSize(dof, dim);
394  Jrt.SetSize(dim);
395  Jpt.SetSize(dim);
396  P.SetSize(dim);
397  PMatI.UseExternalData(elfun.GetData(), dof, dim);
398  elvect.SetSize(dof*dim);
399  PMatO.UseExternalData(elvect.GetData(), dof, dim);
400 
401  const IntegrationRule *ir = IntRule;
402  if (!ir)
403  {
404  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
405  }
406 
407  elvect = 0.0;
408  model->SetTransformation(Ttr);
409  for (int i = 0; i < ir->GetNPoints(); i++)
410  {
411  const IntegrationPoint &ip = ir->IntPoint(i);
412  Ttr.SetIntPoint(&ip);
413  CalcInverse(Ttr.Jacobian(), Jrt);
414 
415  el.CalcDShape(ip, DSh);
416  Mult(DSh, Jrt, DS);
417  MultAtB(PMatI, DS, Jpt);
418 
419  model->EvalP(Jpt, P);
420 
421  P *= ip.weight * Ttr.Weight();
422  AddMultABt(DS, P, PMatO);
423  }
424 }
425 
428  const Vector &elfun,
429  DenseMatrix &elmat)
430 {
431  int dof = el.GetDof(), dim = el.GetDim();
432 
433  DSh.SetSize(dof, dim);
434  DS.SetSize(dof, dim);
435  Jrt.SetSize(dim);
436  Jpt.SetSize(dim);
437  PMatI.UseExternalData(elfun.GetData(), dof, dim);
438  elmat.SetSize(dof*dim);
439 
440  const IntegrationRule *ir = IntRule;
441  if (!ir)
442  {
443  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
444  }
445 
446  elmat = 0.0;
447  model->SetTransformation(Ttr);
448  for (int i = 0; i < ir->GetNPoints(); i++)
449  {
450  const IntegrationPoint &ip = ir->IntPoint(i);
451  Ttr.SetIntPoint(&ip);
452  CalcInverse(Ttr.Jacobian(), Jrt);
453 
454  el.CalcDShape(ip, DSh);
455  Mult(DSh, Jrt, DS);
456  MultAtB(PMatI, DS, Jpt);
457 
458  model->AssembleH(Jpt, DS, ip.weight * Ttr.Weight(), elmat);
459  }
460 }
461 
465  const Array<const Vector *>&elfun)
466 {
467  if (el.Size() != 2)
468  {
469  mfem_error("IncompressibleNeoHookeanIntegrator::GetElementEnergy"
470  " has incorrect block finite element space size!");
471  }
472 
473  int dof_u = el[0]->GetDof();
474  int dim = el[0]->GetDim();
475 
476  DSh_u.SetSize(dof_u, dim);
477  J0i.SetSize(dim);
478  J1.SetSize(dim);
479  J.SetSize(dim);
480  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
481 
482  int intorder = 2*el[0]->GetOrder() + 3; // <---
483  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
484 
485  double energy = 0.0;
486  double mu = 0.0;
487 
488  for (int i = 0; i < ir.GetNPoints(); ++i)
489  {
490  const IntegrationPoint &ip = ir.IntPoint(i);
491  Tr.SetIntPoint(&ip);
492  CalcInverse(Tr.Jacobian(), J0i);
493 
494  el[0]->CalcDShape(ip, DSh_u);
495  MultAtB(PMatI_u, DSh_u, J1);
496  Mult(J1, J0i, J);
497 
498  mu = c_mu->Eval(Tr, ip);
499 
500  energy += ip.weight*Tr.Weight()*(mu/2.0)*(J*J - 3);
501  }
502 
503  return energy;
504 }
505 
509  const Array<const Vector *> &elfun,
510  const Array<Vector *> &elvec)
511 {
512  if (el.Size() != 2)
513  {
514  mfem_error("IncompressibleNeoHookeanIntegrator::AssembleElementVector"
515  " has finite element space of incorrect block number");
516  }
517 
518  int dof_u = el[0]->GetDof();
519  int dof_p = el[1]->GetDof();
520 
521  int dim = el[0]->GetDim();
522  int spaceDim = Tr.GetSpaceDim();
523 
524  if (dim != spaceDim)
525  {
526  mfem_error("IncompressibleNeoHookeanIntegrator::AssembleElementVector"
527  " is not defined on manifold meshes");
528  }
529 
530 
531  DSh_u.SetSize(dof_u, dim);
532  DS_u.SetSize(dof_u, dim);
533  J0i.SetSize(dim);
534  F.SetSize(dim);
535  FinvT.SetSize(dim);
536  P.SetSize(dim);
537  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
538  elvec[0]->SetSize(dof_u*dim);
539  PMatO_u.UseExternalData(elvec[0]->GetData(), dof_u, dim);
540 
541  Sh_p.SetSize(dof_p);
542  elvec[1]->SetSize(dof_p);
543 
544  int intorder = 2*el[0]->GetOrder() + 3; // <---
545  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
546 
547  *elvec[0] = 0.0;
548  *elvec[1] = 0.0;
549 
550  for (int i = 0; i < ir.GetNPoints(); ++i)
551  {
552  const IntegrationPoint &ip = ir.IntPoint(i);
553  Tr.SetIntPoint(&ip);
554  CalcInverse(Tr.Jacobian(), J0i);
555 
556  el[0]->CalcDShape(ip, DSh_u);
557  Mult(DSh_u, J0i, DS_u);
558  MultAtB(PMatI_u, DS_u, F);
559 
560  el[1]->CalcShape(ip, Sh_p);
561 
562  double pres = Sh_p * *elfun[1];
563  double mu = c_mu->Eval(Tr, ip);
564  double dJ = F.Det();
565 
566  CalcInverseTranspose(F, FinvT);
567 
568  P = 0.0;
569  P.Add(mu * dJ, F);
570  P.Add(-1.0 * pres * dJ, FinvT);
571  P *= ip.weight*Tr.Weight();
572 
573  AddMultABt(DS_u, P, PMatO_u);
574 
575  elvec[1]->Add(ip.weight * Tr.Weight() * (dJ - 1.0), Sh_p);
576  }
577 
578 }
579 
581  const Array<const FiniteElement*> &el,
583  const Array<const Vector *> &elfun,
584  const Array2D<DenseMatrix *> &elmats)
585 {
586  int dof_u = el[0]->GetDof();
587  int dof_p = el[1]->GetDof();
588 
589  int dim = el[0]->GetDim();
590 
591  elmats(0,0)->SetSize(dof_u*dim, dof_u*dim);
592  elmats(0,1)->SetSize(dof_u*dim, dof_p);
593  elmats(1,0)->SetSize(dof_p, dof_u*dim);
594  elmats(1,1)->SetSize(dof_p, dof_p);
595 
596  *elmats(0,0) = 0.0;
597  *elmats(0,1) = 0.0;
598  *elmats(1,0) = 0.0;
599  *elmats(1,1) = 0.0;
600 
601  DSh_u.SetSize(dof_u, dim);
602  DS_u.SetSize(dof_u, dim);
603  J0i.SetSize(dim);
604  F.SetSize(dim);
605  FinvT.SetSize(dim);
606  Finv.SetSize(dim);
607  P.SetSize(dim);
608  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
609  Sh_p.SetSize(dof_p);
610 
611  int intorder = 2*el[0]->GetOrder() + 3; // <---
612  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
613 
614  for (int i = 0; i < ir.GetNPoints(); ++i)
615  {
616  const IntegrationPoint &ip = ir.IntPoint(i);
617  Tr.SetIntPoint(&ip);
618  CalcInverse(Tr.Jacobian(), J0i);
619 
620  el[0]->CalcDShape(ip, DSh_u);
621  Mult(DSh_u, J0i, DS_u);
622  MultAtB(PMatI_u, DS_u, F);
623 
624  el[1]->CalcShape(ip, Sh_p);
625  double pres = Sh_p * *elfun[1];
626  double mu = c_mu->Eval(Tr, ip);
627  double dJ = F.Det();
628  double dJ_FinvT_DS;
629 
630  CalcInverseTranspose(F, FinvT);
631 
632  // u,u block
633  for (int i_u = 0; i_u < dof_u; ++i_u)
634  {
635  for (int i_dim = 0; i_dim < dim; ++i_dim)
636  {
637  for (int j_u = 0; j_u < dof_u; ++j_u)
638  {
639  for (int j_dim = 0; j_dim < dim; ++j_dim)
640  {
641 
642  // m = j_dim;
643  // k = i_dim;
644 
645  for (int n=0; n<dim; ++n)
646  {
647  for (int l=0; l<dim; ++l)
648  {
649  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
650  dJ * (mu * F(i_dim, l) - pres * FinvT(i_dim,l)) *
651  FinvT(j_dim,n) * DS_u(i_u,l) * DS_u(j_u, n) *
652  ip.weight * Tr.Weight();
653 
654  if (j_dim == i_dim && n==l)
655  {
656  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
657  dJ * mu * DS_u(i_u, l) * DS_u(j_u,n) *
658  ip.weight * Tr.Weight();
659  }
660 
661  // a = n;
662  // b = m;
663  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
664  dJ * pres * FinvT(i_dim, n) *
665  FinvT(j_dim,l) * DS_u(i_u,l) * DS_u(j_u,n) *
666  ip.weight * Tr.Weight();
667  }
668  }
669  }
670  }
671  }
672  }
673 
674  // u,p and p,u blocks
675  for (int i_p = 0; i_p < dof_p; ++i_p)
676  {
677  for (int j_u = 0; j_u < dof_u; ++j_u)
678  {
679  for (int dim_u = 0; dim_u < dim; ++dim_u)
680  {
681  for (int l=0; l<dim; ++l)
682  {
683  dJ_FinvT_DS = dJ * FinvT(dim_u,l) * DS_u(j_u, l) * Sh_p(i_p) *
684  ip.weight * Tr.Weight();
685  (*elmats(1,0))(i_p, j_u + dof_u * dim_u) += dJ_FinvT_DS;
686  (*elmats(0,1))(j_u + dof_u * dim_u, i_p) -= dJ_FinvT_DS;
687 
688  }
689  }
690  }
691  }
692  }
693 
694 }
695 
696 const IntegrationRule&
699 {
700  const int order = 2 * fe.GetOrder() + T.OrderGrad(&fe);
701  return IntRules.Get(fe.GetGeomType(), order);
702 }
703 
705  const FiniteElement &el,
707  const Vector &elfun,
708  Vector &elvect)
709 {
710  const int nd = el.GetDof();
711  const int dim = el.GetDim();
712 
713  shape.SetSize(nd);
714  dshape.SetSize(nd, dim);
715  elvect.SetSize(nd * dim);
716  gradEF.SetSize(dim);
717 
718  EF.UseExternalData(elfun.GetData(), nd, dim);
719  ELV.UseExternalData(elvect.GetData(), nd, dim);
720 
721  Vector vec1(dim), vec2(dim);
722  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, T);
723  ELV = 0.0;
724  for (int i = 0; i < ir->GetNPoints(); i++)
725  {
726  const IntegrationPoint &ip = ir->IntPoint(i);
727  T.SetIntPoint(&ip);
728  el.CalcShape(ip, shape);
729  el.CalcPhysDShape(T, dshape);
730  double w = ip.weight * T.Weight();
731  if (Q) { w *= Q->Eval(T, ip); }
732  MultAtB(EF, dshape, gradEF);
733  EF.MultTranspose(shape, vec1);
734  gradEF.Mult(vec1, vec2);
735  vec2 *= w;
736  AddMultVWt(shape, vec2, ELV);
737  }
738 }
739 
741  const FiniteElement &el,
743  const Vector &elfun,
744  DenseMatrix &elmat)
745 {
746  int nd = el.GetDof();
747  int dim = el.GetDim();
748 
749  shape.SetSize(nd);
750  dshape.SetSize(nd, dim);
751  dshapex.SetSize(nd, dim);
752  elmat.SetSize(nd * dim);
753  elmat_comp.SetSize(nd);
754  gradEF.SetSize(dim);
755 
756  EF.UseExternalData(elfun.GetData(), nd, dim);
757 
758  double w;
759  Vector vec1(dim), vec2(dim), vec3(nd);
760 
761  const IntegrationRule *ir = IntRule;
762  if (ir == nullptr)
763  {
764  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
765  ir = &IntRules.Get(el.GetGeomType(), order);
766  }
767 
768  elmat = 0.0;
769  for (int i = 0; i < ir->GetNPoints(); i++)
770  {
771  const IntegrationPoint &ip = ir->IntPoint(i);
772  trans.SetIntPoint(&ip);
773 
774  el.CalcShape(ip, shape);
775  el.CalcDShape(ip, dshape);
776 
777  Mult(dshape, trans.InverseJacobian(), dshapex);
778 
779  w = ip.weight;
780 
781  if (Q)
782  {
783  w *= Q->Eval(trans, ip);
784  }
785 
786  MultAtB(EF, dshapex, gradEF);
787  EF.MultTranspose(shape, vec1);
788 
789  trans.AdjugateJacobian().Mult(vec1, vec2);
790 
791  vec2 *= w;
792  dshape.Mult(vec2, vec3);
793  MultVWt(shape, vec3, elmat_comp);
794 
795  for (int i = 0; i < dim; i++)
796  {
797  elmat.AddMatrix(elmat_comp, i * nd, i * nd);
798  }
799 
800  MultVVt(shape, elmat_comp);
801  w = ip.weight * trans.Weight();
802  if (Q)
803  {
804  w *= Q->Eval(trans, ip);
805  }
806  for (int i = 0; i < dim; i++)
807  {
808  for (int j = 0; j < dim; j++)
809  {
810  elmat.AddMatrix(w * gradEF(i, j), elmat_comp, i * nd, j * nd);
811  }
812  }
813  }
814 }
815 
816 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Abstract class for all finite elements.
Definition: fe.hpp:235
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
Definition: nonlininteg.cpp:31
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:2393
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition: eltrans.hpp:127
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W = W(Jpt).
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:309
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:421
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2778
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual int OrderGrad(const FiniteElement *fe) const =0
Return the order of .
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:2759
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:89
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:915
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
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...
double Det() const
Definition: densemat.cpp:451
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:71
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:319
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:132
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:90
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
virtual void AssembleFaceGrad(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:312
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
void SetSize(int m, int n)
Definition: array.hpp:352
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:203
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2163
void EvalCoeffs() const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:65
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:544
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:153
double b
Definition: lissajous.cpp:42
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
Definition: nonlininteg.cpp:18
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:61
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe.cpp:201
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: nonlininteg.cpp:37
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun)
Computes the integral of W(Jacobian(Trt)) over a target zone.
Dynamic 2D array using row-major layout.
Definition: array.hpp:333
void SetTransformation(ElementTransformation &_Ttr)
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2199
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:53
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:2748
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:427
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2498
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:315
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:654
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2331
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0&lt;=i&lt;A.Height, 0&lt;=j&lt;A.Width.
Definition: densemat.cpp:1667
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:45
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2266
ElementTransformation * Ttr
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:29
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
double a
Definition: lissajous.cpp:41
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...
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(Jpt).
double mu
Definition: ex25.cpp:135
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...
int dim
Definition: ex24.cpp:53
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.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
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:2650
Vector data type.
Definition: vector.hpp:51
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:175
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
Definition: nonlininteg.cpp:70
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:65
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
static const IntegrationRule & GetRule(const FiniteElement &fe, ElementTransformation &T)
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...
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:79
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:378