MFEM  v4.5.1
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-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  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 }
Abstract class for all finite elements.
Definition: fe_base.hpp:235
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
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:587
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:1486
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:545
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:513
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:926
Array2D< SparseMatrix * > cGrads
virtual void Setup()
Setup the NonlinearForm: based on the current AssemblyLevel and the current mesh, optionally...
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
void SyncToBlocks() const
Synchronize the memory location flags (i.e. the memory validity flags) of the big/monolithic block-ve...
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.
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:1268
Array< const SparseMatrix * > cP
Array of results of dynamic-casting P to SparseMatrix pointer.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:118
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:598
void BooleanMultTranspose(int alpha, const int *x, int beta, int *y)
The &quot;Boolean&quot; analog of y = alpha * A^T * x + beta * y, where elements in the sparsity pattern of the...
Definition: hypre.hpp:724
virtual double GetEnergy(const Vector &x) const
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:289
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5376
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Boundary Face Integrator.
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
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 Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
virtual void Assemble()=0
Assemble at the AssemblyLevel of the subclass.
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
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.
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
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:479
const SparseMatrix * cP
The result of dynamic-casting P to SparseMatrix pointer.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
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 int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
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:50
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3213
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.
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:639
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2718
Dynamic 2D array using row-major layout.
Definition: array.hpp:351
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:551
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:647
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:270
virtual void Update()
Update the NonlinearForm to propagate updates of the associated FE space.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
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:1850
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
virtual ~BlockNonlinearForm()
Destructor.
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual ~NonlinearForm()
Destroy the NonlinearForm 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:630
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:465
double GetGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
virtual void Mult(const Vector &x, Vector &y) const
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:164
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:904
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:2781
virtual void Update()=0
Called by NonlinearForm::Update() to reflect changes in the FE space.
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.
void SyncFromBlocks() const
Synchronize the memory location flags (i.e. the memory validity flags) of the big/monolithic block-ve...
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
void ComputeGradientBlocked(const BlockVector &bx) const
Specialized version of GetGradient() for BlockVector.
SparseMatrix * Grad
void TransformDual(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition: doftrans.cpp:93
const BlockVector & Prolongate(const BlockVector &bx) const
RefCoord s[3]
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
Definition: fespace.hpp:922
A class to handle Block systems in a matrix-free implementation.
virtual void TransformDual(double *v) const =0
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).
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