MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nonlinearform.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 
14 namespace mfem
15 {
16 
18 {
19  if (ext)
20  {
21  MFEM_ABORT("the assembly level has already been set!");
22  }
23  assembly = assembly_level;
24  switch (assembly)
25  {
27  // This is the default behavior.
28  break;
30  ext = new PANonlinearFormExtension(this);
31  break;
32  default:
33  mfem_error("Unknown assembly level for this form.");
34  }
35 }
36 void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
37  Vector *rhs)
38 {
39  // virtual call, works in parallel too
40  fes->GetEssentialTrueDofs(bdr_attr_is_ess, ess_tdof_list);
41 
42  if (rhs)
43  {
44  for (int i = 0; i < ess_tdof_list.Size(); i++)
45  {
46  (*rhs)(ess_tdof_list[i]) = 0.0;
47  }
48  }
49 }
50 
51 void NonlinearForm::SetEssentialVDofs(const Array<int> &ess_vdofs_list)
52 {
53  if (!P)
54  {
55  ess_vdofs_list.Copy(ess_tdof_list); // ess_vdofs_list --> ess_tdof_list
56  }
57  else
58  {
59  Array<int> ess_vdof_marker, ess_tdof_marker;
61  ess_vdof_marker);
62  if (Serial())
63  {
64  fes->ConvertToConformingVDofs(ess_vdof_marker, ess_tdof_marker);
65  }
66  else
67  {
68 #ifdef MFEM_USE_MPI
69  ParFiniteElementSpace *pf = dynamic_cast<ParFiniteElementSpace*>(fes);
70  ess_tdof_marker.SetSize(pf->GetTrueVSize());
71  pf->Dof_TrueDof_Matrix()->BooleanMultTranspose(1, ess_vdof_marker,
72  0, ess_tdof_marker);
73 #else
74  MFEM_ABORT("internal MFEM error");
75 #endif
76  }
78  }
79 }
80 
82 {
83  Array<int> vdofs;
84  Vector el_x;
85  const FiniteElement *fe;
87  double energy = 0.0;
88 
89  if (dnfi.Size())
90  {
91  for (int i = 0; i < fes->GetNE(); i++)
92  {
93  fe = fes->GetFE(i);
94  fes->GetElementVDofs(i, vdofs);
96  x.GetSubVector(vdofs, el_x);
97  for (int k = 0; k < dnfi.Size(); k++)
98  {
99  energy += dnfi[k]->GetElementEnergy(*fe, *T, el_x);
100  }
101  }
102  }
103 
104  if (fnfi.Size())
105  {
106  MFEM_ABORT("TODO: add energy contribution from interior face terms");
107  }
108 
109  if (bfnfi.Size())
110  {
111  MFEM_ABORT("TODO: add energy contribution from boundary face terms");
112  }
113 
114  return energy;
115 }
116 
118 {
119  MFEM_VERIFY(x.Size() == Width(), "invalid input Vector size");
120  if (P)
121  {
122  aux1.SetSize(P->Height());
123  P->Mult(x, aux1);
124  return aux1;
125  }
126  return x;
127 }
128 
129 void NonlinearForm::Mult(const Vector &x, Vector &y) const
130 {
131  const Vector &px = Prolongate(x);
132  if (P) { aux2.SetSize(P->Height()); }
133 
134  // If we are in parallel, ParNonLinearForm::Mult uses the aux2 vector. In
135  // serial, place the result directly in y.
136  Vector &py = P ? aux2 : y;
137 
138  if (ext)
139  {
140  ext->Mult(px, py);
141  return;
142  }
143 
144  Array<int> vdofs;
145  Vector el_x, el_y;
146  const FiniteElement *fe;
148  Mesh *mesh = fes->GetMesh();
149 
150  py = 0.0;
151 
152  if (dnfi.Size())
153  {
154  for (int i = 0; i < fes->GetNE(); i++)
155  {
156  fe = fes->GetFE(i);
157  fes->GetElementVDofs(i, vdofs);
159  px.GetSubVector(vdofs, el_x);
160  for (int k = 0; k < dnfi.Size(); k++)
161  {
162  dnfi[k]->AssembleElementVector(*fe, *T, el_x, el_y);
163  py.AddElementVector(vdofs, el_y);
164  }
165  }
166  }
167 
168  if (fnfi.Size())
169  {
171  const FiniteElement *fe1, *fe2;
172  Array<int> vdofs2;
173 
174  for (int i = 0; i < mesh->GetNumFaces(); i++)
175  {
176  tr = mesh->GetInteriorFaceTransformations(i);
177  if (tr != NULL)
178  {
179  fes->GetElementVDofs(tr->Elem1No, vdofs);
180  fes->GetElementVDofs(tr->Elem2No, vdofs2);
181  vdofs.Append (vdofs2);
182 
183  px.GetSubVector(vdofs, el_x);
184 
185  fe1 = fes->GetFE(tr->Elem1No);
186  fe2 = fes->GetFE(tr->Elem2No);
187 
188  for (int k = 0; k < fnfi.Size(); k++)
189  {
190  fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
191  py.AddElementVector(vdofs, el_y);
192  }
193  }
194  }
195  }
196 
197  if (bfnfi.Size())
198  {
200  const FiniteElement *fe1, *fe2;
201 
202  // Which boundary attributes need to be processed?
203  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
204  mesh->bdr_attributes.Max() : 0);
205  bdr_attr_marker = 0;
206  for (int k = 0; k < bfnfi.Size(); k++)
207  {
208  if (bfnfi_marker[k] == NULL)
209  {
210  bdr_attr_marker = 1;
211  break;
212  }
213  Array<int> &bdr_marker = *bfnfi_marker[k];
214  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
215  "invalid boundary marker for boundary face integrator #"
216  << k << ", counting from zero");
217  for (int i = 0; i < bdr_attr_marker.Size(); i++)
218  {
219  bdr_attr_marker[i] |= bdr_marker[i];
220  }
221  }
222 
223  for (int i = 0; i < fes -> GetNBE(); i++)
224  {
225  const int bdr_attr = mesh->GetBdrAttribute(i);
226  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
227 
228  tr = mesh->GetBdrFaceTransformations (i);
229  if (tr != NULL)
230  {
231  fes->GetElementVDofs(tr->Elem1No, vdofs);
232  px.GetSubVector(vdofs, el_x);
233 
234  fe1 = fes->GetFE(tr->Elem1No);
235  // The fe2 object is really a dummy and not used on the boundaries,
236  // but we can't dereference a NULL pointer, and we don't want to
237  // actually make a fake element.
238  fe2 = fe1;
239  for (int k = 0; k < bfnfi.Size(); k++)
240  {
241  if (bfnfi_marker[k] &&
242  (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
243 
244  bfnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
245  py.AddElementVector(vdofs, el_y);
246  }
247  }
248  }
249  }
250 
251  if (Serial())
252  {
253  if (cP) { cP->MultTranspose(py, y); }
254 
255  for (int i = 0; i < ess_tdof_list.Size(); i++)
256  {
257  y(ess_tdof_list[i]) = 0.0;
258  }
259  // y(ess_tdof_list[i]) = x(ess_tdof_list[i]);
260  }
261 }
262 
264 {
265  if (ext)
266  {
267  MFEM_ABORT("Not yet implemented!");
268  }
269 
270  const int skip_zeros = 0;
271  Array<int> vdofs;
272  Vector el_x;
273  DenseMatrix elmat;
274  const FiniteElement *fe;
276  Mesh *mesh = fes->GetMesh();
277  const Vector &px = Prolongate(x);
278 
279  if (Grad == NULL)
280  {
281  Grad = new SparseMatrix(fes->GetVSize());
282  }
283  else
284  {
285  *Grad = 0.0;
286  }
287 
288  if (dnfi.Size())
289  {
290  for (int i = 0; i < fes->GetNE(); i++)
291  {
292  fe = fes->GetFE(i);
293  fes->GetElementVDofs(i, vdofs);
295  px.GetSubVector(vdofs, el_x);
296  for (int k = 0; k < dnfi.Size(); k++)
297  {
298  dnfi[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
299  Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
300  // Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
301  }
302  }
303  }
304 
305  if (fnfi.Size())
306  {
308  const FiniteElement *fe1, *fe2;
309  Array<int> vdofs2;
310 
311  for (int i = 0; i < mesh->GetNumFaces(); i++)
312  {
313  tr = mesh->GetInteriorFaceTransformations(i);
314  if (tr != NULL)
315  {
316  fes->GetElementVDofs(tr->Elem1No, vdofs);
317  fes->GetElementVDofs(tr->Elem2No, vdofs2);
318  vdofs.Append (vdofs2);
319 
320  px.GetSubVector(vdofs, el_x);
321 
322  fe1 = fes->GetFE(tr->Elem1No);
323  fe2 = fes->GetFE(tr->Elem2No);
324 
325  for (int k = 0; k < fnfi.Size(); k++)
326  {
327  fnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
328  Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
329  }
330  }
331  }
332  }
333 
334  if (bfnfi.Size())
335  {
337  const FiniteElement *fe1, *fe2;
338 
339  // Which boundary attributes need to be processed?
340  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
341  mesh->bdr_attributes.Max() : 0);
342  bdr_attr_marker = 0;
343  for (int k = 0; k < bfnfi.Size(); k++)
344  {
345  if (bfnfi_marker[k] == NULL)
346  {
347  bdr_attr_marker = 1;
348  break;
349  }
350  Array<int> &bdr_marker = *bfnfi_marker[k];
351  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
352  "invalid boundary marker for boundary face integrator #"
353  << k << ", counting from zero");
354  for (int i = 0; i < bdr_attr_marker.Size(); i++)
355  {
356  bdr_attr_marker[i] |= bdr_marker[i];
357  }
358  }
359 
360  for (int i = 0; i < fes -> GetNBE(); i++)
361  {
362  const int bdr_attr = mesh->GetBdrAttribute(i);
363  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
364 
365  tr = mesh->GetBdrFaceTransformations (i);
366  if (tr != NULL)
367  {
368  fes->GetElementVDofs(tr->Elem1No, vdofs);
369  px.GetSubVector(vdofs, el_x);
370 
371  fe1 = fes->GetFE(tr->Elem1No);
372  // The fe2 object is really a dummy and not used on the boundaries,
373  // but we can't dereference a NULL pointer, and we don't want to
374  // actually make a fake element.
375  fe2 = fe1;
376  for (int k = 0; k < bfnfi.Size(); k++)
377  {
378  if (bfnfi_marker[k] &&
379  (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
380 
381  bfnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
382  Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
383  }
384  }
385  }
386  }
387 
388  if (!Grad->Finalized())
389  {
390  Grad->Finalize(skip_zeros);
391  }
392 
393  SparseMatrix *mGrad = Grad;
394  if (Serial())
395  {
396  if (cP)
397  {
398  delete cGrad;
399  cGrad = RAP(*cP, *Grad, *cP);
400  mGrad = cGrad;
401  }
402  for (int i = 0; i < ess_tdof_list.Size(); i++)
403  {
404  mGrad->EliminateRowCol(ess_tdof_list[i]);
405  }
406  }
407 
408  return *mGrad;
409 }
410 
412 {
413  if (ext) { MFEM_ABORT("Not yet implemented!"); }
414 
415  if (sequence == fes->GetSequence()) { return; }
416 
417  height = width = fes->GetTrueVSize();
418  delete cGrad; cGrad = NULL;
419  delete Grad; Grad = NULL;
420  ess_tdof_list.SetSize(0); // essential b.c. will need to be set again
421  sequence = fes->GetSequence();
422  // Do not modify aux1 and aux2, their size will be set before use.
424  cP = dynamic_cast<const SparseMatrix*>(P);
425 }
426 
428 {
429  if (ext) { return ext->AssemblePA(); }
430 }
431 
433 {
434  delete cGrad;
435  delete Grad;
436  for (int i = 0; i < dnfi.Size(); i++) { delete dnfi[i]; }
437  for (int i = 0; i < fnfi.Size(); i++) { delete fnfi[i]; }
438  for (int i = 0; i < bfnfi.Size(); i++) { delete bfnfi[i]; }
439  delete ext;
440 }
441 
442 
444  fes(0), BlockGrad(NULL)
445 {
446  height = 0;
447  width = 0;
448 }
449 
451 {
452  delete BlockGrad;
453  BlockGrad = NULL;
454  for (int i=0; i<Grads.NumRows(); ++i)
455  {
456  for (int j=0; j<Grads.NumCols(); ++j)
457  {
458  delete Grads(i,j);
459  delete cGrads(i,j);
460  }
461  }
462  for (int i = 0; i < ess_tdofs.Size(); ++i)
463  {
464  delete ess_tdofs[i];
465  }
466 
467  height = 0;
468  width = 0;
469  f.Copy(fes);
470  block_offsets.SetSize(f.Size() + 1);
471  block_trueOffsets.SetSize(f.Size() + 1);
472  block_offsets[0] = 0;
473  block_trueOffsets[0] = 0;
474 
475  for (int i=0; i<fes.Size(); ++i)
476  {
477  block_offsets[i+1] = fes[i]->GetVSize();
478  block_trueOffsets[i+1] = fes[i]->GetTrueVSize();
479  }
480 
483 
486 
487  Grads.SetSize(fes.Size(), fes.Size());
488  Grads = NULL;
489 
490  cGrads.SetSize(fes.Size(), fes.Size());
491  cGrads = NULL;
492 
493  P.SetSize(fes.Size());
494  cP.SetSize(fes.Size());
495  ess_tdofs.SetSize(fes.Size());
496  for (int s = 0; s < fes.Size(); ++s)
497  {
498  // Retrieve prolongation matrix for each FE space
499  P[s] = fes[s]->GetProlongationMatrix();
500  cP[s] = dynamic_cast<const SparseMatrix *>(P[s]);
501 
502  // If the P Operator exists and its type is not SparseMatrix, this
503  // indicates the Operator is part of parallel run.
504  if (P[s] && !cP[s])
505  {
506  is_serial = false;
507  }
508 
509  // If the P Operator exists and its type is SparseMatrix, this indicates
510  // the Operator is serial but needs prolongation on assembly.
511  if (cP[s])
512  {
513  needs_prolongation = true;
514  }
515 
516  ess_tdofs[s] = new Array<int>;
517  }
518 }
519 
521  fes(0), BlockGrad(NULL)
522 {
523  SetSpaces(f);
524 }
525 
527  Array<int> &bdr_attr_marker)
528 {
529  bfnfi.Append(nfi);
530  bfnfi_marker.Append(&bdr_attr_marker);
531 }
532 
534  const Array<Array<int> *> &bdr_attr_is_ess, Array<Vector *> &rhs)
535 {
536  for (int s = 0; s < fes.Size(); ++s)
537  {
538  ess_tdofs[s]->SetSize(ess_tdofs.Size());
539 
540  fes[s]->GetEssentialTrueDofs(*bdr_attr_is_ess[s], *ess_tdofs[s]);
541 
542  if (rhs[s])
543  {
544  rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
545  }
546  }
547 }
548 
550 {
551  Array<Array<int> *> vdofs(fes.Size());
552  Array<Vector *> el_x(fes.Size());
553  Array<const Vector *> el_x_const(fes.Size());
556  double energy = 0.0;
557 
558  for (int i=0; i<fes.Size(); ++i)
559  {
560  el_x_const[i] = el_x[i] = new Vector();
561  vdofs[i] = new Array<int>;
562  }
563 
564  if (dnfi.Size())
565  for (int i = 0; i < fes[0]->GetNE(); ++i)
566  {
567  T = fes[0]->GetElementTransformation(i);
568  for (int s=0; s<fes.Size(); ++s)
569  {
570  fe[s] = fes[s]->GetFE(i);
571  fes[s]->GetElementVDofs(i, *vdofs[s]);
572  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
573  }
574 
575  for (int k = 0; k < dnfi.Size(); ++k)
576  {
577  energy += dnfi[k]->GetElementEnergy(fe, *T, el_x_const);
578  }
579  }
580 
581  // free the allocated memory
582  for (int i = 0; i < fes.Size(); ++i)
583  {
584  delete el_x[i];
585  delete vdofs[i];
586  }
587 
588  if (fnfi.Size())
589  {
590  MFEM_ABORT("TODO: add energy contribution from interior face terms");
591  }
592 
593  if (bfnfi.Size())
594  {
595  MFEM_ABORT("TODO: add energy contribution from boundary face terms");
596  }
597 
598  return energy;
599 }
600 
601 double BlockNonlinearForm::GetEnergy(const Vector &x) const
602 {
604  return GetEnergyBlocked(xs);
605 }
606 
608  BlockVector &by) const
609 {
610  Array<Array<int> *>vdofs(fes.Size());
611  Array<Array<int> *>vdofs2(fes.Size());
612  Array<Vector *> el_x(fes.Size());
613  Array<const Vector *> el_x_const(fes.Size());
614  Array<Vector *> el_y(fes.Size());
618 
619  by = 0.0;
620  for (int s=0; s<fes.Size(); ++s)
621  {
622  el_x_const[s] = el_x[s] = new Vector();
623  el_y[s] = new Vector();
624  vdofs[s] = new Array<int>;
625  vdofs2[s] = new Array<int>;
626  }
627 
628  if (dnfi.Size())
629  {
630  for (int i = 0; i < fes[0]->GetNE(); ++i)
631  {
632  T = fes[0]->GetElementTransformation(i);
633  for (int s = 0; s < fes.Size(); ++s)
634  {
635  fes[s]->GetElementVDofs(i, *(vdofs[s]));
636  fe[s] = fes[s]->GetFE(i);
637  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
638  }
639 
640  for (int k = 0; k < dnfi.Size(); ++k)
641  {
642  dnfi[k]->AssembleElementVector(fe, *T,
643  el_x_const, el_y);
644 
645  for (int s=0; s<fes.Size(); ++s)
646  {
647  if (el_y[s]->Size() == 0) { continue; }
648  by.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
649  }
650  }
651  }
652  }
653 
654  if (fnfi.Size())
655  {
656  Mesh *mesh = fes[0]->GetMesh();
658 
659  for (int i = 0; i < mesh->GetNumFaces(); ++i)
660  {
661  tr = mesh->GetInteriorFaceTransformations(i);
662  if (tr != NULL)
663  {
664  for (int s=0; s<fes.Size(); ++s)
665  {
666  fe[s] = fes[s]->GetFE(tr->Elem1No);
667  fe2[s] = fes[s]->GetFE(tr->Elem2No);
668 
669  fes[s]->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
670  fes[s]->GetElementVDofs(tr->Elem2No, *(vdofs2[s]));
671 
672  vdofs[s]->Append(*(vdofs2[s]));
673 
674  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
675  }
676 
677  for (int k = 0; k < fnfi.Size(); ++k)
678  {
679 
680  fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
681 
682  for (int s=0; s<fes.Size(); ++s)
683  {
684  if (el_y[s]->Size() == 0) { continue; }
685  by.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
686  }
687  }
688  }
689  }
690  }
691 
692  if (bfnfi.Size())
693  {
694  Mesh *mesh = fes[0]->GetMesh();
696  // Which boundary attributes need to be processed?
697  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
698  mesh->bdr_attributes.Max() : 0);
699  bdr_attr_marker = 0;
700  for (int k = 0; k < bfnfi.Size(); ++k)
701  {
702  if (bfnfi_marker[k] == NULL)
703  {
704  bdr_attr_marker = 1;
705  break;
706  }
707  Array<int> &bdr_marker = *bfnfi_marker[k];
708  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
709  "invalid boundary marker for boundary face integrator #"
710  << k << ", counting from zero");
711  for (int i = 0; i < bdr_attr_marker.Size(); ++i)
712  {
713  bdr_attr_marker[i] |= bdr_marker[i];
714  }
715  }
716 
717  for (int i = 0; i < mesh->GetNBE(); ++i)
718  {
719  const int bdr_attr = mesh->GetBdrAttribute(i);
720  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
721 
722  tr = mesh->GetBdrFaceTransformations(i);
723  if (tr != NULL)
724  {
725  for (int s=0; s<fes.Size(); ++s)
726  {
727  fe[s] = fes[s]->GetFE(tr->Elem1No);
728  fe2[s] = fes[s]->GetFE(tr->Elem1No);
729 
730  fes[s]->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
731  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
732  }
733 
734  for (int k = 0; k < bfnfi.Size(); ++k)
735  {
736  if (bfnfi_marker[k] &&
737  (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
738 
739  bfnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
740 
741  for (int s=0; s<fes.Size(); ++s)
742  {
743  if (el_y[s]->Size() == 0) { continue; }
744  by.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
745  }
746  }
747  }
748  }
749  }
750 
751  for (int s=0; s<fes.Size(); ++s)
752  {
753  delete vdofs2[s];
754  delete vdofs[s];
755  delete el_y[s];
756  delete el_x[s];
757  }
758 }
759 
761 {
762  MFEM_VERIFY(bx.Size() == Width(), "invalid input BlockVector size");
763 
764  if (needs_prolongation)
765  {
767  for (int s = 0; s < fes.Size(); s++)
768  {
769  P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
770  }
771  return aux1;
772  }
773  return bx;
774 }
775 
776 void BlockNonlinearForm::Mult(const Vector &x, Vector &y) const
777 {
780 
781  const BlockVector &pbx = Prolongate(bx);
782  if (needs_prolongation)
783  {
785  }
786  BlockVector &pby = needs_prolongation ? aux2 : by;
787 
788  xs.Update(pbx.GetData(), block_offsets);
789  ys.Update(pby.GetData(), block_offsets);
790  MultBlocked(xs, ys);
791 
792  for (int s = 0; s < fes.Size(); s++)
793  {
794  if (cP[s])
795  {
796  cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
797  }
798  by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
799  }
800 }
801 
803 {
804  const int skip_zeros = 0;
805  Array<Array<int> *> vdofs(fes.Size());
806  Array<Array<int> *> vdofs2(fes.Size());
807  Array<Vector *> el_x(fes.Size());
808  Array<const Vector *> el_x_const(fes.Size());
809  Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
813 
814  for (int i=0; i<fes.Size(); ++i)
815  {
816  el_x_const[i] = el_x[i] = new Vector();
817  vdofs[i] = new Array<int>;
818  vdofs2[i] = new Array<int>;
819  for (int j=0; j<fes.Size(); ++j)
820  {
821  elmats(i,j) = new DenseMatrix();
822  }
823  }
824 
825  for (int i=0; i<fes.Size(); ++i)
826  {
827  for (int j=0; j<fes.Size(); ++j)
828  {
829  if (Grads(i,j) != NULL)
830  {
831  *Grads(i,j) = 0.0;
832  }
833  else
834  {
835  Grads(i,j) = new SparseMatrix(fes[i]->GetVSize(),
836  fes[j]->GetVSize());
837  }
838  }
839  }
840 
841  if (dnfi.Size())
842  {
843  for (int i = 0; i < fes[0]->GetNE(); ++i)
844  {
845  T = fes[0]->GetElementTransformation(i);
846  for (int s = 0; s < fes.Size(); ++s)
847  {
848  fe[s] = fes[s]->GetFE(i);
849  fes[s]->GetElementVDofs(i, *vdofs[s]);
850  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
851  }
852 
853  for (int k = 0; k < dnfi.Size(); ++k)
854  {
855  dnfi[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
856 
857  for (int j=0; j<fes.Size(); ++j)
858  {
859  for (int l=0; l<fes.Size(); ++l)
860  {
861  if (elmats(j,l)->Height() == 0) { continue; }
862  Grads(j,l)->AddSubMatrix(*vdofs[j], *vdofs[l],
863  *elmats(j,l), skip_zeros);
864  }
865  }
866  }
867  }
868  }
869 
870  if (fnfi.Size())
871  {
873  Mesh *mesh = fes[0]->GetMesh();
874 
875  for (int i = 0; i < mesh->GetNumFaces(); ++i)
876  {
877  tr = mesh->GetInteriorFaceTransformations(i);
878 
879  for (int s=0; s < fes.Size(); ++s)
880  {
881  fe[s] = fes[s]->GetFE(tr->Elem1No);
882  fe2[s] = fes[s]->GetFE(tr->Elem2No);
883 
884  fes[s]->GetElementVDofs(tr->Elem1No, *vdofs[s]);
885  fes[s]->GetElementVDofs(tr->Elem2No, *vdofs2[s]);
886  vdofs[s]->Append(*(vdofs2[s]));
887 
888  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
889  }
890 
891  for (int k = 0; k < fnfi.Size(); ++k)
892  {
893  fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
894  for (int j=0; j<fes.Size(); ++j)
895  {
896  for (int l=0; l<fes.Size(); ++l)
897  {
898  if (elmats(j,l)->Height() == 0) { continue; }
899  Grads(j,l)->AddSubMatrix(*vdofs[j], *vdofs[l],
900  *elmats(j,l), skip_zeros);
901  }
902  }
903  }
904  }
905  }
906 
907  if (bfnfi.Size())
908  {
910  Mesh *mesh = fes[0]->GetMesh();
911 
912  // Which boundary attributes need to be processed?
913  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
914  mesh->bdr_attributes.Max() : 0);
915  bdr_attr_marker = 0;
916  for (int k = 0; k < bfnfi.Size(); ++k)
917  {
918  if (bfnfi_marker[k] == NULL)
919  {
920  bdr_attr_marker = 1;
921  break;
922  }
923  Array<int> &bdr_marker = *bfnfi_marker[k];
924  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
925  "invalid boundary marker for boundary face integrator #"
926  << k << ", counting from zero");
927  for (int i = 0; i < bdr_attr_marker.Size(); ++i)
928  {
929  bdr_attr_marker[i] |= bdr_marker[i];
930  }
931  }
932 
933  for (int i = 0; i < mesh->GetNBE(); ++i)
934  {
935  const int bdr_attr = mesh->GetBdrAttribute(i);
936  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
937 
938  tr = mesh->GetBdrFaceTransformations(i);
939  if (tr != NULL)
940  {
941  for (int s = 0; s < fes.Size(); ++s)
942  {
943  fe[s] = fes[s]->GetFE(tr->Elem1No);
944  fe2[s] = fe[s];
945 
946  fes[s]->GetElementVDofs(i, *vdofs[s]);
947  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
948  }
949 
950  for (int k = 0; k < bfnfi.Size(); ++k)
951  {
952  bfnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
953  for (int l=0; l<fes.Size(); ++l)
954  {
955  for (int j=0; j<fes.Size(); ++j)
956  {
957  if (elmats(j,l)->Height() == 0) { continue; }
958  Grads(j,l)->AddSubMatrix(*vdofs[j], *vdofs[l],
959  *elmats(j,l), skip_zeros);
960  }
961  }
962  }
963  }
964  }
965  }
966 
967  if (!Grads(0,0)->Finalized())
968  {
969  for (int i=0; i<fes.Size(); ++i)
970  {
971  for (int j=0; j<fes.Size(); ++j)
972  {
973  Grads(i,j)->Finalize(skip_zeros);
974  }
975  }
976  }
977 
978  for (int i=0; i<fes.Size(); ++i)
979  {
980  for (int j=0; j<fes.Size(); ++j)
981  {
982  delete elmats(i,j);
983  }
984  delete vdofs2[i];
985  delete vdofs[i];
986  delete el_x[i];
987  }
988 }
989 
991 {
993  const BlockVector &pbx = Prolongate(bx);
994 
996 
997  Array2D<SparseMatrix *> mGrads(fes.Size(), fes.Size());
998  mGrads = Grads;
999  if (needs_prolongation)
1000  {
1001  for (int s1 = 0; s1 < fes.Size(); ++s1)
1002  {
1003  for (int s2 = 0; s2 < fes.Size(); ++s2)
1004  {
1005  delete cGrads(s1, s2);
1006  cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
1007  mGrads(s1, s2) = cGrads(s1, s2);
1008  }
1009  }
1010  }
1011 
1012  for (int s = 0; s < fes.Size(); ++s)
1013  {
1014  for (int i = 0; i < ess_tdofs[s]->Size(); ++i)
1015  {
1016  for (int j = 0; j < fes.Size(); ++j)
1017  {
1018  if (s == j)
1019  {
1020  mGrads(s, s)->EliminateRowCol((*ess_tdofs[s])[i],
1022  }
1023  else
1024  {
1025  mGrads(s, j)->EliminateRow((*ess_tdofs[s])[i]);
1026  mGrads(j, s)->EliminateCol((*ess_tdofs[s])[i]);
1027  }
1028  }
1029  }
1030  }
1031 
1032  delete BlockGrad;
1034  for (int i = 0; i < fes.Size(); ++i)
1035  {
1036  for (int j = 0; j < fes.Size(); ++j)
1037  {
1038  BlockGrad->SetBlock(i, j, mGrads(i, j));
1039  }
1040  }
1041  return *BlockGrad;
1042 }
1043 
1045 {
1046  delete BlockGrad;
1047  for (int i=0; i<fes.Size(); ++i)
1048  {
1049  for (int j=0; j<fes.Size(); ++j)
1050  {
1051  delete Grads(i,j);
1052  delete cGrads(i,j);
1053  }
1054  delete ess_tdofs[i];
1055  }
1056 
1057  for (int i = 0; i < dnfi.Size(); ++i)
1058  {
1059  delete dnfi[i];
1060  }
1061 
1062  for (int i = 0; i < fnfi.Size(); ++i)
1063  {
1064  delete fnfi[i];
1065  }
1066 
1067  for (int i = 0; i < bfnfi.Size(); ++i)
1068  {
1069  delete bfnfi[i];
1070  }
1071 
1072 }
1073 
1074 }
Abstract class for all finite elements.
Definition: fe.hpp:235
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
Array< const Operator * > P
Array of pointers to the prolongation matrix of fes, may be NULL.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:400
BlockNonlinearForm()
Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
bool is_serial
Indicator if the Operator is part of a parallel run.
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1053
Array< NonlinearFormIntegrator * > dnfi
Set of Domain Integrators to be assembled (added).
Array< BlockNonlinearFormIntegrator * > dnfi
Set of Domain Integrators to be assembled (added).
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1640
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:460
SparseMatrix * cGrad
virtual void SetEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:740
Array2D< SparseMatrix * > cGrads
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:173
virtual void Setup()
Setup the NonlinearForm.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:71
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:831
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:495
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
void Update(double *data, const Array< int > &bOffsets)
Update method.
Definition: blockvector.cpp:77
Array< NonlinearFormIntegrator * > bfnfi
Set of boundary face Integrators to be assembled (added).
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:415
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Array< BlockNonlinearFormIntegrator * > bfnfi
Set of Boundary Face Integrators to be assembled (added).
Data and methods for partially-assembled nonlinear forms.
Array< int > ess_tdof_list
A list of all essential true dofs.
FiniteElementSpace * fes
FE space on which the form lives.
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
Definition: mesh.hpp:1020
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
Array< const SparseMatrix * > cP
Array of results of dynamic-casting P to SparseMatrix pointer.
Array2D< SparseMatrix * > Grads
void SetSpaces(Array< FiniteElementSpace * > &f)
(Re)initialize the BlockNonlinearForm.
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition: fespace.cpp:434
void BooleanMultTranspose(int alpha, const int *x, int beta, int *y)
Definition: hypre.hpp:466
virtual double GetEnergy(const Vector &x) const
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:255
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:4184
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:427
void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Boundary Face Integrator.
Array< FiniteElementSpace * > fes
FE spaces on which the form lives.
Data type sparse matrix.
Definition: sparsemat.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:65
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:726
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:314
Array< NonlinearFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
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
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:281
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:330
const SparseMatrix * cP
The result of dynamic-casting P to SparseMatrix pointer.
Vector aux1
Auxiliary Vectors.
Array< BlockNonlinearFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
AssemblyLevel assembly
The assembly level.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:403
const Operator * P
Pointer to the prolongation matrix of fes, may be NULL.
Array< Array< int > * > bfnfi_marker
Set the diagonal value to one.
Definition: operator.hpp:48
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1003
Array< Array< int > * > bfnfi_marker
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::NONE.
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition: vector.cpp:587
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2436
Dynamic 2D array using row-major layout.
Definition: array.hpp:333
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:466
bool Serial() const
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:460
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
virtual void Update()
Update the NonlinearForm to propagate updates of the associated FE space.
virtual void AssemblePA()=0
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:654
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1624
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
NonlinearFormExtension * ext
virtual ~BlockNonlinearForm()
Destructor.
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual ~NonlinearForm()
Destroy the NoninearForm including the owned NonlinearFormIntegrators and gradient Operator...
void ConvertToConformingVDofs(const Array< int > &dofs, Array< int > &cdofs)
For a partially conforming FE space, convert a marker array (nonzero entries are true) on the partial...
Definition: fespace.cpp:466
double GetGridFunctionEnergy(const Vector &x) const
Compute the enery corresponding to the state x.
virtual void Mult(const Vector &x, Vector &y) const
virtual Operator & GetGradient(const Vector &x) const
Compute the gradient Operator of the NonlinearForm corresponding to the state x.
double GetEnergyBlocked(const BlockVector &bx) const
Specialized version of GetEnergy() for BlockVectors.
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:721
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:1798
long sequence
Counter for updates propagated from the FiniteElementSpace.
bool needs_prolongation
Indicator if the Operator needs prolongation on assembly.
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
Definition: fespace.cpp:453
Vector data type.
Definition: vector.hpp:51
void ComputeGradientBlocked(const BlockVector &bx) const
Specialized version of GetGradient() for BlockVector.
SparseMatrix * Grad
const BlockVector & Prolongate(const BlockVector &bx) const
virtual Operator & GetGradient(const Vector &x) const
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
Array< Array< int > * > ess_tdofs
Abstract operator.
Definition: operator.hpp:24
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:707
A class to handle Block systems in a matrix-free implementation.
const Vector & Prolongate(const Vector &x) const
void MultBlocked(const BlockVector &bx, BlockVector &by) const
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:87
double f(const Vector &p)