MFEM  v4.5.1
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-2022, 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::GetLocalStateEnergyPA(...)\n"
21  " is not implemented for this class.");
22  return 0.0;
23 }
24 
26 {
27  mfem_error ("NonlinearFormIntegrator::AssemblePA(...)\n"
28  " is not implemented for this class.");
29 }
30 
32  const FiniteElementSpace &)
33 {
34  mfem_error ("NonlinearFormIntegrator::AssemblePA(...)\n"
35  " is not implemented for this class.");
36 }
37 
39  const FiniteElementSpace &fes)
40 {
41  mfem_error ("NonlinearFormIntegrator::AssembleGradPA(...)\n"
42  " is not implemented for this class.");
43 }
44 
46 {
47  mfem_error ("NonlinearFormIntegrator::AddMultPA(...)\n"
48  " is not implemented for this class.");
49 }
50 
52 {
53  mfem_error ("NonlinearFormIntegrator::AddMultGradPA(...)\n"
54  " is not implemented for this class.");
55 }
56 
58 {
59  mfem_error ("NonlinearFormIntegrator::AssembleGradDiagonalPA(...)\n"
60  " is not implemented for this class.");
61 }
62 
64 {
65  mfem_error ("NonlinearFormIntegrator::AssembleMF(...)\n"
66  " is not implemented for this class.");
67 }
68 
70 {
71  mfem_error ("NonlinearFormIntegrator::AddMultMF(...)\n"
72  " is not implemented for this class.");
73 }
74 
76  const FiniteElement &el, ElementTransformation &Tr,
77  const Vector &elfun, Vector &elvect)
78 {
79  mfem_error("NonlinearFormIntegrator::AssembleElementVector"
80  " is not overloaded!");
81 }
82 
84  const FiniteElement &el1, const FiniteElement &el2,
85  FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
86 {
87  mfem_error("NonlinearFormIntegrator::AssembleFaceVector"
88  " is not overloaded!");
89 }
90 
92  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
93  DenseMatrix &elmat)
94 {
95  mfem_error("NonlinearFormIntegrator::AssembleElementGrad"
96  " is not overloaded!");
97 }
98 
100  const FiniteElement &el1, const FiniteElement &el2,
101  FaceElementTransformations &Tr, const Vector &elfun,
102  DenseMatrix &elmat)
103 {
104  mfem_error("NonlinearFormIntegrator::AssembleFaceGrad"
105  " is not overloaded!");
106 }
107 
109  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
110 {
111  mfem_error("NonlinearFormIntegrator::GetElementEnergy"
112  " is not overloaded!");
113  return 0.0;
114 }
115 
116 
120  const Array<const Vector *> &elfun,
121  const Array<Vector *> &elvec)
122 {
123  mfem_error("BlockNonlinearFormIntegrator::AssembleElementVector"
124  " is not overloaded!");
125 }
126 
128  const Array<const FiniteElement *> &el1,
129  const Array<const FiniteElement *> &el2,
131  const Array<const Vector *> &elfun,
132  const Array<Vector *> &elvect)
133 {
134  mfem_error("BlockNonlinearFormIntegrator::AssembleFaceVector"
135  " is not overloaded!");
136 }
137 
139  const Array<const FiniteElement*> &el,
141  const Array<const Vector *> &elfun,
142  const Array2D<DenseMatrix *> &elmats)
143 {
144  mfem_error("BlockNonlinearFormIntegrator::AssembleElementGrad"
145  " is not overloaded!");
146 }
147 
152  const Array<const Vector *> &elfun,
153  const Array2D<DenseMatrix *> &elmats)
154 {
155  mfem_error("BlockNonlinearFormIntegrator::AssembleFaceGrad"
156  " is not overloaded!");
157 }
158 
162  const Array<const Vector *>&elfun)
163 {
164  mfem_error("BlockNonlinearFormIntegrator::GetElementEnergy"
165  " is not overloaded!");
166  return 0.0;
167 }
168 
169 
171 {
172  Z.SetSize(J.Width());
174  return 0.5*(Z*Z)/J.Det();
175 }
176 
178 {
179  int dim = J.Width();
180  double t;
181 
182  Z.SetSize(dim);
183  S.SetSize(dim);
185  MultAAt(Z, S);
186  t = 0.5*S.Trace();
187  for (int i = 0; i < dim; i++)
188  {
189  S(i,i) -= t;
190  }
191  t = J.Det();
192  S *= -1.0/(t*t);
193  Mult(S, Z, P);
194 }
195 
197  const DenseMatrix &J, const DenseMatrix &DS, const double weight,
198  DenseMatrix &A) const
199 {
200  int dof = DS.Height(), dim = DS.Width();
201  double t;
202 
203  Z.SetSize(dim);
204  S.SetSize(dim);
205  G.SetSize(dof, dim);
206  C.SetSize(dof, dim);
207 
209  MultAAt(Z, S);
210 
211  t = 1.0/J.Det();
212  Z *= t; // Z = J^{-t}
213  S *= t; // S = |J| (J.J^t)^{-1}
214  t = 0.5*S.Trace();
215 
216  MultABt(DS, Z, G); // G = DS.J^{-1}
217  Mult(G, S, C);
218 
219  // 1.
220  for (int i = 0; i < dof; i++)
221  for (int j = 0; j <= i; j++)
222  {
223  double a = 0.0;
224  for (int d = 0; d < dim; d++)
225  {
226  a += G(i,d)*G(j,d);
227  }
228  a *= weight;
229  for (int k = 0; k < dim; k++)
230  for (int l = 0; l <= k; l++)
231  {
232  double b = a*S(k,l);
233  A(i+k*dof,j+l*dof) += b;
234  if (i != j)
235  {
236  A(j+k*dof,i+l*dof) += b;
237  }
238  if (k != l)
239  {
240  A(i+l*dof,j+k*dof) += b;
241  if (i != j)
242  {
243  A(j+l*dof,i+k*dof) += b;
244  }
245  }
246  }
247  }
248 
249  // 2.
250  for (int i = 1; i < dof; i++)
251  for (int j = 0; j < i; j++)
252  {
253  for (int k = 1; k < dim; k++)
254  for (int l = 0; l < k; l++)
255  {
256  double a =
257  weight*(C(i,l)*G(j,k) - C(i,k)*G(j,l) +
258  C(j,k)*G(i,l) - C(j,l)*G(i,k) +
259  t*(G(i,k)*G(j,l) - G(i,l)*G(j,k)));
260 
261  A(i+k*dof,j+l*dof) += a;
262  A(j+l*dof,i+k*dof) += a;
263 
264  A(i+l*dof,j+k*dof) -= a;
265  A(j+k*dof,i+l*dof) -= a;
266  }
267  }
268 }
269 
270 
271 inline void NeoHookeanModel::EvalCoeffs() const
272 {
273  mu = c_mu->Eval(*Ttr, Ttr->GetIntPoint());
274  K = c_K->Eval(*Ttr, Ttr->GetIntPoint());
275  if (c_g)
276  {
277  g = c_g->Eval(*Ttr, Ttr->GetIntPoint());
278  }
279 }
280 
281 double NeoHookeanModel::EvalW(const DenseMatrix &J) const
282 {
283  int dim = J.Width();
284 
285  if (have_coeffs)
286  {
287  EvalCoeffs();
288  }
289 
290  double dJ = J.Det();
291  double sJ = dJ/g;
292  double bI1 = pow(dJ, -2.0/dim)*(J*J); // \bar{I}_1
293 
294  return 0.5*(mu*(bI1 - dim) + K*(sJ - 1.0)*(sJ - 1.0));
295 }
296 
298 {
299  int dim = J.Width();
300 
301  if (have_coeffs)
302  {
303  EvalCoeffs();
304  }
305 
306  Z.SetSize(dim);
308 
309  double dJ = J.Det();
310  double a = mu*pow(dJ, -2.0/dim);
311  double b = K*(dJ/g - 1.0)/g - a*(J*J)/(dim*dJ);
312 
313  P = 0.0;
314  P.Add(a, J);
315  P.Add(b, Z);
316 }
317 
319  const double weight, DenseMatrix &A) const
320 {
321  int dof = DS.Height(), dim = DS.Width();
322 
323  if (have_coeffs)
324  {
325  EvalCoeffs();
326  }
327 
328  Z.SetSize(dim);
329  G.SetSize(dof, dim);
330  C.SetSize(dof, dim);
331 
332  double dJ = J.Det();
333  double sJ = dJ/g;
334  double a = mu*pow(dJ, -2.0/dim);
335  double bc = a*(J*J)/dim;
336  double b = bc - K*sJ*(sJ - 1.0);
337  double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0);
338 
340  Z *= (1.0/dJ); // Z = J^{-t}
341 
342  MultABt(DS, J, C); // C = DS J^t
343  MultABt(DS, Z, G); // G = DS J^{-1}
344 
345  a *= weight;
346  b *= weight;
347  c *= weight;
348 
349  // 1.
350  for (int i = 0; i < dof; i++)
351  for (int k = 0; k <= i; k++)
352  {
353  double s = 0.0;
354  for (int d = 0; d < dim; d++)
355  {
356  s += DS(i,d)*DS(k,d);
357  }
358  s *= a;
359 
360  for (int d = 0; d < dim; d++)
361  {
362  A(i+d*dof,k+d*dof) += s;
363  }
364 
365  if (k != i)
366  for (int d = 0; d < dim; d++)
367  {
368  A(k+d*dof,i+d*dof) += s;
369  }
370  }
371 
372  a *= (-2.0/dim);
373 
374  // 2.
375  for (int i = 0; i < dof; i++)
376  for (int j = 0; j < dim; j++)
377  for (int k = 0; k < dof; k++)
378  for (int l = 0; l < dim; l++)
379  {
380  A(i+j*dof,k+l*dof) +=
381  a*(C(i,j)*G(k,l) + G(i,j)*C(k,l)) +
382  b*G(i,l)*G(k,j) + c*G(i,j)*G(k,l);
383  }
384 }
385 
386 
389  const Vector &elfun)
390 {
391  int dof = el.GetDof(), dim = el.GetDim();
392  double energy;
393 
394  DSh.SetSize(dof, dim);
395  Jrt.SetSize(dim);
396  Jpr.SetSize(dim);
397  Jpt.SetSize(dim);
398  PMatI.UseExternalData(elfun.GetData(), dof, dim);
399 
400  const IntegrationRule *ir = IntRule;
401  if (!ir)
402  {
403  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
404  }
405 
406  energy = 0.0;
407  model->SetTransformation(Ttr);
408  for (int i = 0; i < ir->GetNPoints(); i++)
409  {
410  const IntegrationPoint &ip = ir->IntPoint(i);
411  Ttr.SetIntPoint(&ip);
412  CalcInverse(Ttr.Jacobian(), Jrt);
413 
414  el.CalcDShape(ip, DSh);
415  MultAtB(PMatI, DSh, Jpr);
416  Mult(Jpr, Jrt, Jpt);
417 
418  energy += ip.weight * Ttr.Weight() * model->EvalW(Jpt);
419  }
420 
421  return energy;
422 }
423 
425  const FiniteElement &el, ElementTransformation &Ttr,
426  const Vector &elfun, Vector &elvect)
427 {
428  int dof = el.GetDof(), dim = el.GetDim();
429 
430  DSh.SetSize(dof, dim);
431  DS.SetSize(dof, dim);
432  Jrt.SetSize(dim);
433  Jpt.SetSize(dim);
434  P.SetSize(dim);
435  PMatI.UseExternalData(elfun.GetData(), dof, dim);
436  elvect.SetSize(dof*dim);
437  PMatO.UseExternalData(elvect.GetData(), dof, dim);
438 
439  const IntegrationRule *ir = IntRule;
440  if (!ir)
441  {
442  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
443  }
444 
445  elvect = 0.0;
446  model->SetTransformation(Ttr);
447  for (int i = 0; i < ir->GetNPoints(); i++)
448  {
449  const IntegrationPoint &ip = ir->IntPoint(i);
450  Ttr.SetIntPoint(&ip);
451  CalcInverse(Ttr.Jacobian(), Jrt);
452 
453  el.CalcDShape(ip, DSh);
454  Mult(DSh, Jrt, DS);
455  MultAtB(PMatI, DS, Jpt);
456 
457  model->EvalP(Jpt, P);
458 
459  P *= ip.weight * Ttr.Weight();
460  AddMultABt(DS, P, PMatO);
461  }
462 }
463 
466  const Vector &elfun,
467  DenseMatrix &elmat)
468 {
469  int dof = el.GetDof(), dim = el.GetDim();
470 
471  DSh.SetSize(dof, dim);
472  DS.SetSize(dof, dim);
473  Jrt.SetSize(dim);
474  Jpt.SetSize(dim);
475  PMatI.UseExternalData(elfun.GetData(), dof, dim);
476  elmat.SetSize(dof*dim);
477 
478  const IntegrationRule *ir = IntRule;
479  if (!ir)
480  {
481  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
482  }
483 
484  elmat = 0.0;
485  model->SetTransformation(Ttr);
486  for (int i = 0; i < ir->GetNPoints(); i++)
487  {
488  const IntegrationPoint &ip = ir->IntPoint(i);
489  Ttr.SetIntPoint(&ip);
490  CalcInverse(Ttr.Jacobian(), Jrt);
491 
492  el.CalcDShape(ip, DSh);
493  Mult(DSh, Jrt, DS);
494  MultAtB(PMatI, DS, Jpt);
495 
496  model->AssembleH(Jpt, DS, ip.weight * Ttr.Weight(), elmat);
497  }
498 }
499 
503  const Array<const Vector *>&elfun)
504 {
505  if (el.Size() != 2)
506  {
507  mfem_error("IncompressibleNeoHookeanIntegrator::GetElementEnergy"
508  " has incorrect block finite element space size!");
509  }
510 
511  int dof_u = el[0]->GetDof();
512  int dim = el[0]->GetDim();
513 
514  DSh_u.SetSize(dof_u, dim);
515  J0i.SetSize(dim);
516  J1.SetSize(dim);
517  J.SetSize(dim);
518  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
519 
520  int intorder = 2*el[0]->GetOrder() + 3; // <---
521  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
522 
523  double energy = 0.0;
524  double mu = 0.0;
525 
526  for (int i = 0; i < ir.GetNPoints(); ++i)
527  {
528  const IntegrationPoint &ip = ir.IntPoint(i);
529  Tr.SetIntPoint(&ip);
530  CalcInverse(Tr.Jacobian(), J0i);
531 
532  el[0]->CalcDShape(ip, DSh_u);
533  MultAtB(PMatI_u, DSh_u, J1);
534  Mult(J1, J0i, J);
535 
536  mu = c_mu->Eval(Tr, ip);
537 
538  energy += ip.weight*Tr.Weight()*(mu/2.0)*(J*J - 3);
539  }
540 
541  return energy;
542 }
543 
547  const Array<const Vector *> &elfun,
548  const Array<Vector *> &elvec)
549 {
550  if (el.Size() != 2)
551  {
552  mfem_error("IncompressibleNeoHookeanIntegrator::AssembleElementVector"
553  " has finite element space of incorrect block number");
554  }
555 
556  int dof_u = el[0]->GetDof();
557  int dof_p = el[1]->GetDof();
558 
559  int dim = el[0]->GetDim();
560  int spaceDim = Tr.GetSpaceDim();
561 
562  if (dim != spaceDim)
563  {
564  mfem_error("IncompressibleNeoHookeanIntegrator::AssembleElementVector"
565  " is not defined on manifold meshes");
566  }
567 
568 
569  DSh_u.SetSize(dof_u, dim);
570  DS_u.SetSize(dof_u, dim);
571  J0i.SetSize(dim);
572  F.SetSize(dim);
573  FinvT.SetSize(dim);
574  P.SetSize(dim);
575  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
576  elvec[0]->SetSize(dof_u*dim);
577  PMatO_u.UseExternalData(elvec[0]->GetData(), dof_u, dim);
578 
579  Sh_p.SetSize(dof_p);
580  elvec[1]->SetSize(dof_p);
581 
582  int intorder = 2*el[0]->GetOrder() + 3; // <---
583  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
584 
585  *elvec[0] = 0.0;
586  *elvec[1] = 0.0;
587 
588  for (int i = 0; i < ir.GetNPoints(); ++i)
589  {
590  const IntegrationPoint &ip = ir.IntPoint(i);
591  Tr.SetIntPoint(&ip);
592  CalcInverse(Tr.Jacobian(), J0i);
593 
594  el[0]->CalcDShape(ip, DSh_u);
595  Mult(DSh_u, J0i, DS_u);
596  MultAtB(PMatI_u, DS_u, F);
597 
598  el[1]->CalcShape(ip, Sh_p);
599 
600  double pres = Sh_p * *elfun[1];
601  double mu = c_mu->Eval(Tr, ip);
602  double dJ = F.Det();
603 
604  CalcInverseTranspose(F, FinvT);
605 
606  P = 0.0;
607  P.Add(mu * dJ, F);
608  P.Add(-1.0 * pres * dJ, FinvT);
609  P *= ip.weight*Tr.Weight();
610 
611  AddMultABt(DS_u, P, PMatO_u);
612 
613  elvec[1]->Add(ip.weight * Tr.Weight() * (dJ - 1.0), Sh_p);
614  }
615 
616 }
617 
619  const Array<const FiniteElement*> &el,
621  const Array<const Vector *> &elfun,
622  const Array2D<DenseMatrix *> &elmats)
623 {
624  int dof_u = el[0]->GetDof();
625  int dof_p = el[1]->GetDof();
626 
627  int dim = el[0]->GetDim();
628 
629  elmats(0,0)->SetSize(dof_u*dim, dof_u*dim);
630  elmats(0,1)->SetSize(dof_u*dim, dof_p);
631  elmats(1,0)->SetSize(dof_p, dof_u*dim);
632  elmats(1,1)->SetSize(dof_p, dof_p);
633 
634  *elmats(0,0) = 0.0;
635  *elmats(0,1) = 0.0;
636  *elmats(1,0) = 0.0;
637  *elmats(1,1) = 0.0;
638 
639  DSh_u.SetSize(dof_u, dim);
640  DS_u.SetSize(dof_u, dim);
641  J0i.SetSize(dim);
642  F.SetSize(dim);
643  FinvT.SetSize(dim);
644  Finv.SetSize(dim);
645  P.SetSize(dim);
646  PMatI_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
647  Sh_p.SetSize(dof_p);
648 
649  int intorder = 2*el[0]->GetOrder() + 3; // <---
650  const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
651 
652  for (int i = 0; i < ir.GetNPoints(); ++i)
653  {
654  const IntegrationPoint &ip = ir.IntPoint(i);
655  Tr.SetIntPoint(&ip);
656  CalcInverse(Tr.Jacobian(), J0i);
657 
658  el[0]->CalcDShape(ip, DSh_u);
659  Mult(DSh_u, J0i, DS_u);
660  MultAtB(PMatI_u, DS_u, F);
661 
662  el[1]->CalcShape(ip, Sh_p);
663  double pres = Sh_p * *elfun[1];
664  double mu = c_mu->Eval(Tr, ip);
665  double dJ = F.Det();
666  double dJ_FinvT_DS;
667 
668  CalcInverseTranspose(F, FinvT);
669 
670  // u,u block
671  for (int i_u = 0; i_u < dof_u; ++i_u)
672  {
673  for (int i_dim = 0; i_dim < dim; ++i_dim)
674  {
675  for (int j_u = 0; j_u < dof_u; ++j_u)
676  {
677  for (int j_dim = 0; j_dim < dim; ++j_dim)
678  {
679 
680  // m = j_dim;
681  // k = i_dim;
682 
683  for (int n=0; n<dim; ++n)
684  {
685  for (int l=0; l<dim; ++l)
686  {
687  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
688  dJ * (mu * F(i_dim, l) - pres * FinvT(i_dim,l)) *
689  FinvT(j_dim,n) * DS_u(i_u,l) * DS_u(j_u, n) *
690  ip.weight * Tr.Weight();
691 
692  if (j_dim == i_dim && n==l)
693  {
694  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
695  dJ * mu * DS_u(i_u, l) * DS_u(j_u,n) *
696  ip.weight * Tr.Weight();
697  }
698 
699  // a = n;
700  // b = m;
701  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
702  dJ * pres * FinvT(i_dim, n) *
703  FinvT(j_dim,l) * DS_u(i_u,l) * DS_u(j_u,n) *
704  ip.weight * Tr.Weight();
705  }
706  }
707  }
708  }
709  }
710  }
711 
712  // u,p and p,u blocks
713  for (int i_p = 0; i_p < dof_p; ++i_p)
714  {
715  for (int j_u = 0; j_u < dof_u; ++j_u)
716  {
717  for (int dim_u = 0; dim_u < dim; ++dim_u)
718  {
719  for (int l=0; l<dim; ++l)
720  {
721  dJ_FinvT_DS = dJ * FinvT(dim_u,l) * DS_u(j_u, l) * Sh_p(i_p) *
722  ip.weight * Tr.Weight();
723  (*elmats(1,0))(i_p, j_u + dof_u * dim_u) += dJ_FinvT_DS;
724  (*elmats(0,1))(j_u + dof_u * dim_u, i_p) -= dJ_FinvT_DS;
725 
726  }
727  }
728  }
729  }
730  }
731 
732 }
733 
734 const IntegrationRule&
737 {
738  const int order = 2 * fe.GetOrder() + T.OrderGrad(&fe);
739  return IntRules.Get(fe.GetGeomType(), order);
740 }
741 
743  const FiniteElement &el,
745  const Vector &elfun,
746  Vector &elvect)
747 {
748  const int nd = el.GetDof();
749  dim = el.GetDim();
750 
751  shape.SetSize(nd);
752  dshape.SetSize(nd, dim);
753  elvect.SetSize(nd * dim);
754  gradEF.SetSize(dim);
755 
756  EF.UseExternalData(elfun.GetData(), nd, dim);
757  ELV.UseExternalData(elvect.GetData(), nd, dim);
758 
759  Vector vec1(dim), vec2(dim);
760  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, T);
761  ELV = 0.0;
762  for (int i = 0; i < ir->GetNPoints(); i++)
763  {
764  const IntegrationPoint &ip = ir->IntPoint(i);
765  T.SetIntPoint(&ip);
766  el.CalcShape(ip, shape);
767  el.CalcPhysDShape(T, dshape);
768  double w = ip.weight * T.Weight();
769  if (Q) { w *= Q->Eval(T, ip); }
770 
771  MultAtB(EF, dshape, gradEF);
772  EF.MultTranspose(shape, vec1);
773  gradEF.Mult(vec1, vec2);
774  vec2 *= w;
775  AddMultVWt(shape, vec2, ELV);
776  }
777 }
778 
780  const FiniteElement &el,
782  const Vector &elfun,
783  DenseMatrix &elmat)
784 {
785  const int nd = el.GetDof();
786  dim = el.GetDim();
787 
788  shape.SetSize(nd);
789  dshape.SetSize(nd, dim);
790  dshapex.SetSize(nd, dim);
791  elmat.SetSize(nd * dim);
792  elmat_comp.SetSize(nd);
793  gradEF.SetSize(dim);
794 
795  EF.UseExternalData(elfun.GetData(), nd, dim);
796 
797  double w;
798  Vector vec1(dim), vec2(dim), vec3(nd);
799 
800  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, trans);
801 
802  elmat = 0.0;
803  for (int i = 0; i < ir->GetNPoints(); i++)
804  {
805  const IntegrationPoint &ip = ir->IntPoint(i);
806  trans.SetIntPoint(&ip);
807 
808  el.CalcShape(ip, shape);
809  el.CalcDShape(ip, dshape);
810 
811  Mult(dshape, trans.InverseJacobian(), dshapex);
812 
813  w = ip.weight;
814 
815  if (Q)
816  {
817  w *= Q->Eval(trans, ip);
818  }
819 
820  MultAtB(EF, dshapex, gradEF);
821  EF.MultTranspose(shape, vec1);
822 
823  trans.AdjugateJacobian().Mult(vec1, vec2);
824 
825  vec2 *= w;
826  dshape.Mult(vec2, vec3);
827  MultVWt(shape, vec3, elmat_comp);
828 
829  for (int ii = 0; ii < dim; ii++)
830  {
831  elmat.AddMatrix(elmat_comp, ii * nd, ii * nd);
832  }
833 
834  MultVVt(shape, elmat_comp);
835  w = ip.weight * trans.Weight();
836  if (Q)
837  {
838  w *= Q->Eval(trans, ip);
839  }
840  for (int ii = 0; ii < dim; ii++)
841  {
842  for (int jj = 0; jj < dim; jj++)
843  {
844  elmat.AddMatrix(w * gradEF(ii, jj), elmat_comp, ii * nd, jj * nd);
845  }
846  }
847  }
848 }
849 
850 
852  const FiniteElement &el,
854  const Vector &elfun,
855  DenseMatrix &elmat)
856 {
857  const int nd = el.GetDof();
858  const int dim = el.GetDim();
859 
860  shape.SetSize(nd);
861  dshape.SetSize(nd, dim);
862  dshapex.SetSize(nd, dim);
863  elmat.SetSize(nd * dim);
864  elmat_comp.SetSize(nd);
865  gradEF.SetSize(dim);
866 
867  EF.UseExternalData(elfun.GetData(), nd, dim);
868 
869  Vector vec1(dim), vec2(dim), vec3(nd);
870 
871  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, trans);
872 
873  elmat = 0.0;
874  for (int i = 0; i < ir->GetNPoints(); i++)
875  {
876  const IntegrationPoint &ip = ir->IntPoint(i);
877  trans.SetIntPoint(&ip);
878 
879  el.CalcShape(ip, shape);
880  el.CalcDShape(ip, dshape);
881 
882  const double w = Q ? Q->Eval(trans, ip) * ip.weight : ip.weight;
883 
884  EF.MultTranspose(shape, vec1); // u^n
885 
886  trans.AdjugateJacobian().Mult(vec1, vec2);
887 
888  vec2 *= w;
889  dshape.Mult(vec2, vec3); // (u^n \cdot grad u^{n+1})
890  MultVWt(shape, vec3, elmat_comp); // (u^n \cdot grad u^{n+1},v)
891 
892  for (int ii = 0; ii < dim; ii++)
893  {
894  elmat.AddMatrix(elmat_comp, ii * nd, ii * nd);
895  }
896  }
897 }
898 
899 
901  const FiniteElement &el,
903  const Vector &elfun,
904  DenseMatrix &elmat)
905 {
906  const int nd = el.GetDof();
907  const int dim = el.GetDim();
908 
909  shape.SetSize(nd);
910  dshape.SetSize(nd, dim);
911  dshapex.SetSize(nd, dim);
912  elmat.SetSize(nd * dim);
913  elmat_comp.SetSize(nd);
914  gradEF.SetSize(dim);
915 
916  DenseMatrix elmat_comp_T(nd);
917 
918  EF.UseExternalData(elfun.GetData(), nd, dim);
919 
920  Vector vec1(dim), vec2(dim), vec3(nd), vec4(dim), vec5(nd);
921 
922  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, trans);
923 
924  elmat = 0.0;
925  elmat_comp_T = 0.0;
926  for (int i = 0; i < ir->GetNPoints(); i++)
927  {
928  const IntegrationPoint &ip = ir->IntPoint(i);
929  trans.SetIntPoint(&ip);
930 
931  el.CalcShape(ip, shape);
932  el.CalcDShape(ip, dshape);
933 
934  Mult(dshape, trans.InverseJacobian(), dshapex);
935 
936  const double w = Q ? Q->Eval(trans, ip) * ip.weight : ip.weight;
937 
938  EF.MultTranspose(shape, vec1); // u^n
939 
940  trans.AdjugateJacobian().Mult(vec1, vec2);
941 
942  vec2 *= w;
943  dshape.Mult(vec2, vec3); // (u^n \cdot grad u^{n+1})
944  MultVWt(shape, vec3, elmat_comp); // (u^n \cdot grad u^{n+1},v)
945  elmat_comp_T.Transpose(elmat_comp);
946 
947  for (int ii = 0; ii < dim; ii++)
948  {
949  elmat.AddMatrix(.5, elmat_comp, ii * nd, ii * nd);
950  elmat.AddMatrix(-.5, elmat_comp_T, ii * nd, ii * nd);
951  }
952  }
953 }
954 
955 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
Definition: nonlininteg.cpp:45
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:2742
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition: eltrans.hpp:135
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_base.hpp:311
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
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:3127
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:3108
virtual void AssembleFaceVector(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvect)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
virtual void AssembleGradDiagonalPA(Vector &diag) const
Method for computing the diagonal of the gradient with partial assembly.
Definition: nonlininteg.cpp:57
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
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:449
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
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:145
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:98
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
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.
virtual void AssembleMF(const FiniteElementSpace &fes)
Method defining fully unassembled operator.
Definition: nonlininteg.cpp:63
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
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:370
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:201
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2512
void EvalCoeffs() const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:542
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:154
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:119
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:25
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:99
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
virtual void AddMultGradPA(const Vector &x, Vector &y) const
Method for partially assembled gradient action.
Definition: nonlininteg.cpp:51
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_base.cpp:191
virtual double GetLocalStateEnergyPA(const Vector &x) const
Compute the local (to the MPI rank) energy with partial assembly.
Definition: nonlininteg.cpp:18
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:75
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:351
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2548
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:91
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1381
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:3097
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:425
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2847
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void SetTransformation(ElementTransformation &Ttr_)
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2680
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:1684
virtual void AddMultMF(const Vector &x, Vector &y) const
Definition: nonlininteg.cpp:69
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:83
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2615
ElementTransformation * Ttr
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:30
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:139
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.
RefCoord t[3]
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.
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
RefCoord s[3]
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
virtual void AssembleGradPA(const Vector &x, const FiniteElementSpace &fes)
Prepare the integrator for partial assembly (PA) gradient evaluations on the given FE space fes at th...
Definition: nonlininteg.cpp:38