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