MFEM  v4.2.0 Finite element discretization library
bilinearform.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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 // Implementation of class BilinearForm
13
14 #include "fem.hpp"
15 #include "../general/device.hpp"
16 #include <cmath>
17
18 namespace mfem
19 {
20
22 {
23  if (static_cond) { return; }
24
25  if (precompute_sparsity == 0 || fes->GetVDim() > 1)
26  {
27  mat = new SparseMatrix(height);
28  return;
29  }
30
31  const Table &elem_dof = fes->GetElementToDofTable();
32  Table dof_dof;
33
34  if (fbfi.Size() > 0)
35  {
36  // the sparsity pattern is defined from the map: face->element->dof
37  Table face_dof, dof_face;
38  {
39  Table *face_elem = fes->GetMesh()->GetFaceToElementTable();
40  mfem::Mult(*face_elem, elem_dof, face_dof);
41  delete face_elem;
42  }
43  Transpose(face_dof, dof_face, height);
44  mfem::Mult(dof_face, face_dof, dof_dof);
45  }
46  else
47  {
48  // the sparsity pattern is defined from the map: element->dof
49  Table dof_elem;
50  Transpose(elem_dof, dof_elem, height);
51  mfem::Mult(dof_elem, elem_dof, dof_dof);
52  }
53
54  dof_dof.SortRows();
55
56  int *I = dof_dof.GetI();
57  int *J = dof_dof.GetJ();
58  double *data = Memory<double>(I[height]);
59
60  mat = new SparseMatrix(I, J, data, height, height, true, true, true);
61  *mat = 0.0;
62
63  dof_dof.LoseData();
64 }
65
67  : Matrix (f->GetVSize())
68 {
69  fes = f;
70  sequence = f->GetSequence();
71  mat = mat_e = NULL;
72  extern_bfs = 0;
73  element_matrices = NULL;
74  static_cond = NULL;
75  hybridization = NULL;
78
80  batch = 1;
81  ext = NULL;
82 }
83
85  : Matrix (f->GetVSize())
86 {
87  fes = f;
88  sequence = f->GetSequence();
89  mat_e = NULL;
90  extern_bfs = 1;
91  element_matrices = NULL;
92  static_cond = NULL;
93  hybridization = NULL;
96
98  batch = 1;
99  ext = NULL;
100
101  // Copy the pointers to the integrators
102  dbfi = bf->dbfi;
103
104  bbfi = bf->bbfi;
105  bbfi_marker = bf->bbfi_marker;
106
107  fbfi = bf->fbfi;
108
109  bfbfi = bf->bfbfi;
111
112  AllocMat();
113 }
114
116 {
117  if (ext)
118  {
119  MFEM_ABORT("the assembly level has already been set!");
120  }
121  assembly = assembly_level;
122  switch (assembly)
123  {
125  break;
126  case AssemblyLevel::FULL:
127  ext = new FABilinearFormExtension(this);
128  break;
130  ext = new EABilinearFormExtension(this);
131  break;
133  ext = new PABilinearFormExtension(this);
134  break;
135  case AssemblyLevel::NONE:
136  ext = new MFBilinearFormExtension(this);
137  break;
138  default:
139  mfem_error("Unknown assembly level");
140  }
141 }
142
144 {
145  delete static_cond;
147  {
148  static_cond = NULL;
149  MFEM_WARNING("Static condensation not supported for this assembly level");
150  return;
151  }
154  {
155  bool symmetric = false; // TODO
156  bool block_diagonal = false; // TODO
157  static_cond->Init(symmetric, block_diagonal);
158  }
159  else
160  {
161  delete static_cond;
162  static_cond = NULL;
163  }
164 }
165
167  BilinearFormIntegrator *constr_integ,
168  const Array<int> &ess_tdof_list)
169 {
170  delete hybridization;
172  {
173  delete constr_integ;
174  hybridization = NULL;
175  MFEM_WARNING("Hybridization not supported for this assembly level");
176  return;
177  }
178  hybridization = new Hybridization(fes, constr_space);
179  hybridization->SetConstraintIntegrator(constr_integ);
180  hybridization->Init(ess_tdof_list);
181 }
182
183 void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
184 {
185  if (static_cond) { return; }
186
187  if (mat)
188  {
189  if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
190  {
191  return; // mat is already using the given sparsity
192  }
193  delete mat;
194  }
195  height = width = fes->GetVSize();
196  mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
197 }
198
200 {
201  MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
202  "invalid matrix A dimensions: "
203  << A.Height() << " x " << A.Width());
204  MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
205
206  UseSparsity(A.GetI(), A.GetJ(), A.ColumnsAreSorted());
207 }
208
209 double& BilinearForm::Elem (int i, int j)
210 {
211  return mat -> Elem(i,j);
212 }
213
214 const double& BilinearForm::Elem (int i, int j) const
215 {
216  return mat -> Elem(i,j);
217 }
218
220 {
221  return mat -> Inverse();
222 }
223
224 void BilinearForm::Finalize (int skip_zeros)
225 {
227  {
228  if (!static_cond) { mat->Finalize(skip_zeros); }
229  if (mat_e) { mat_e->Finalize(skip_zeros); }
230  if (static_cond) { static_cond->Finalize(); }
232  }
233 }
234
236 {
237  dbfi.Append(bfi);
238 }
239
241 {
242  bbfi.Append (bfi);
243  bbfi_marker.Append(NULL); // NULL marker means apply everywhere
244 }
245
247  Array<int> &bdr_marker)
248 {
249  bbfi.Append (bfi);
250  bbfi_marker.Append(&bdr_marker);
251 }
252
254 {
255  fbfi.Append (bfi);
256 }
257
259 {
260  bfbfi.Append(bfi);
261  bfbfi_marker.Append(NULL); // NULL marker means apply everywhere
262 }
263
265  Array<int> &bdr_marker)
266 {
267  bfbfi.Append(bfi);
268  bfbfi_marker.Append(&bdr_marker);
269 }
270
272 {
273  if (element_matrices)
274  {
276  elmat = element_matrices->GetData(i);
277  return;
278  }
279
280  if (dbfi.Size())
281  {
282  const FiniteElement &fe = *fes->GetFE(i);
284  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
285  for (int k = 1; k < dbfi.Size(); k++)
286  {
287  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
288  elmat += elemmat;
289  }
290  }
291  else
292  {
294  elmat.SetSize(vdofs.Size());
295  elmat = 0.0;
296  }
297 }
298
300 {
301  if (bbfi.Size())
302  {
303  const FiniteElement &be = *fes->GetBE(i);
305  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
306  for (int k = 1; k < bbfi.Size(); k++)
307  {
308  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
309  elmat += elemmat;
310  }
311  }
312  else
313  {
315  elmat.SetSize(vdofs.Size());
316  elmat = 0.0;
317  }
318 }
319
321  int i, const DenseMatrix &elmat, int skip_zeros)
322 {
323  AssembleElementMatrix(i, elmat, vdofs, skip_zeros);
324 }
325
327  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
328 {
329  fes->GetElementVDofs(i, vdofs);
330  if (static_cond)
331  {
332  static_cond->AssembleMatrix(i, elmat);
333  }
334  else
335  {
336  if (mat == NULL)
337  {
338  AllocMat();
339  }
341  if (hybridization)
342  {
343  hybridization->AssembleMatrix(i, elmat);
344  }
345  }
346 }
347
349  int i, const DenseMatrix &elmat, int skip_zeros)
350 {
351  AssembleBdrElementMatrix(i, elmat, vdofs, skip_zeros);
352 }
353
355  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
356 {
357  fes->GetBdrElementVDofs(i, vdofs);
358  if (static_cond)
359  {
360  static_cond->AssembleBdrMatrix(i, elmat);
361  }
362  else
363  {
364  if (mat == NULL)
365  {
366  AllocMat();
367  }
369  if (hybridization)
370  {
371  hybridization->AssembleBdrMatrix(i, elmat);
372  }
373  }
374 }
375
376 void BilinearForm::Assemble(int skip_zeros)
377 {
378  if (ext)
379  {
380  ext->Assemble();
381  return;
382  }
383
384  ElementTransformation *eltrans;
385  Mesh *mesh = fes -> GetMesh();
386  DenseMatrix elmat, *elmat_p;
387
388  if (mat == NULL)
389  {
390  AllocMat();
391  }
392
393 #ifdef MFEM_USE_LEGACY_OPENMP
394  int free_element_matrices = 0;
395  if (!element_matrices)
396  {
398  free_element_matrices = 1;
399  }
400 #endif
401
402  if (dbfi.Size())
403  {
404  for (int i = 0; i < fes -> GetNE(); i++)
405  {
407  if (element_matrices)
408  {
409  elmat_p = &(*element_matrices)(i);
410  }
411  else
412  {
413  const FiniteElement &fe = *fes->GetFE(i);
414  eltrans = fes->GetElementTransformation(i);
415  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
416  for (int k = 1; k < dbfi.Size(); k++)
417  {
418  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
419  elmat += elemmat;
420  }
421  elmat_p = &elmat;
422  }
423  if (static_cond)
424  {
425  static_cond->AssembleMatrix(i, *elmat_p);
426  }
427  else
428  {
430  if (hybridization)
431  {
432  hybridization->AssembleMatrix(i, *elmat_p);
433  }
434  }
435  }
436  }
437
438  if (bbfi.Size())
439  {
440  // Which boundary attributes need to be processed?
441  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
442  mesh->bdr_attributes.Max() : 0);
443  bdr_attr_marker = 0;
444  for (int k = 0; k < bbfi.Size(); k++)
445  {
446  if (bbfi_marker[k] == NULL)
447  {
448  bdr_attr_marker = 1;
449  break;
450  }
451  Array<int> &bdr_marker = *bbfi_marker[k];
452  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
453  "invalid boundary marker for boundary integrator #"
454  << k << ", counting from zero");
455  for (int i = 0; i < bdr_attr_marker.Size(); i++)
456  {
457  bdr_attr_marker[i] |= bdr_marker[i];
458  }
459  }
460
461  for (int i = 0; i < fes -> GetNBE(); i++)
462  {
463  const int bdr_attr = mesh->GetBdrAttribute(i);
464  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
465
466  const FiniteElement &be = *fes->GetBE(i);
467  fes -> GetBdrElementVDofs (i, vdofs);
468  eltrans = fes -> GetBdrElementTransformation (i);
469  int k = 0;
470  for (; k < bbfi.Size(); k++)
471  {
472  if (bbfi_marker[k] &&
473  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
474
475  bbfi[k]->AssembleElementMatrix(be, *eltrans, elmat);
476  k++;
477  break;
478  }
479  for (; k < bbfi.Size(); k++)
480  {
481  if (bbfi_marker[k] &&
482  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
483
484  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
485  elmat += elemmat;
486  }
487  if (!static_cond)
488  {
490  if (hybridization)
491  {
492  hybridization->AssembleBdrMatrix(i, elmat);
493  }
494  }
495  else
496  {
497  static_cond->AssembleBdrMatrix(i, elmat);
498  }
499  }
500  }
501
502  if (fbfi.Size())
503  {
505  Array<int> vdofs2;
506
507  int nfaces = mesh->GetNumFaces();
508  for (int i = 0; i < nfaces; i++)
509  {
510  tr = mesh -> GetInteriorFaceTransformations (i);
511  if (tr != NULL)
512  {
513  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
514  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
515  vdofs.Append (vdofs2);
516  for (int k = 0; k < fbfi.Size(); k++)
517  {
518  fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
519  *fes -> GetFE (tr -> Elem2No),
520  *tr, elemmat);
521  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
522  }
523  }
524  }
525  }
526
527  if (bfbfi.Size())
528  {
530  const FiniteElement *fe1, *fe2;
531
532  // Which boundary attributes need to be processed?
533  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
534  mesh->bdr_attributes.Max() : 0);
535  bdr_attr_marker = 0;
536  for (int k = 0; k < bfbfi.Size(); k++)
537  {
538  if (bfbfi_marker[k] == NULL)
539  {
540  bdr_attr_marker = 1;
541  break;
542  }
543  Array<int> &bdr_marker = *bfbfi_marker[k];
544  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
545  "invalid boundary marker for boundary face integrator #"
546  << k << ", counting from zero");
547  for (int i = 0; i < bdr_attr_marker.Size(); i++)
548  {
549  bdr_attr_marker[i] |= bdr_marker[i];
550  }
551  }
552
553  for (int i = 0; i < fes -> GetNBE(); i++)
554  {
555  const int bdr_attr = mesh->GetBdrAttribute(i);
556  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
557
558  tr = mesh -> GetBdrFaceTransformations (i);
559  if (tr != NULL)
560  {
561  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
562  fe1 = fes -> GetFE (tr -> Elem1No);
563  // The fe2 object is really a dummy and not used on the boundaries,
564  // but we can't dereference a NULL pointer, and we don't want to
565  // actually make a fake element.
566  fe2 = fe1;
567  for (int k = 0; k < bfbfi.Size(); k++)
568  {
569  if (bfbfi_marker[k] &&
570  (*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
571
572  bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
573  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
574  }
575  }
576  }
577  }
578
579 #ifdef MFEM_USE_LEGACY_OPENMP
580  if (free_element_matrices)
581  {
583  }
584 #endif
585 }
586
588 {
589  // Do not remove zero entries to preserve the symmetric structure of the
590  // matrix which in turn will give rise to symmetric structure in the new
591  // matrix. This ensures that subsequent calls to EliminateRowCol will work
592  // correctly.
593  Finalize(0);
594  MFEM_ASSERT(mat, "the BilinearForm is not assembled");
595
597  if (!P) { return; } // conforming mesh
598
599  SparseMatrix *R = Transpose(*P);
600  SparseMatrix *RA = mfem::Mult(*R, *mat);
601  delete mat;
602  if (mat_e)
603  {
604  SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
605  delete mat_e;
606  mat_e = RAe;
607  }
608  delete R;
609  mat = mfem::Mult(*RA, *P);
610  delete RA;
611  if (mat_e)
612  {
613  SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
614  delete mat_e;
615  mat_e = RAeP;
616  }
617
618  height = mat->Height();
619  width = mat->Width();
620 }
621
623 {
624  if (ext)
625  {
626  MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
627  "Vector for holding diagonal has wrong size!");
628  const Operator *P = fes->GetProlongationMatrix();
629  // For an AMR mesh, a convergent diagonal is assembled with |P^T| d_e,
630  // where |P^T| has the entry-wise absolute values of the conforming
631  // prolongation transpose operator.
632  if (P && !fes->Conforming())
633  {
634  Vector local_diag(P->Height());
635  ext->AssembleDiagonal(local_diag);
636  const SparseMatrix *SP = dynamic_cast<const SparseMatrix*>(P);
637 #ifdef MFEM_USE_MPI
638  const HypreParMatrix *HP = dynamic_cast<const HypreParMatrix*>(P);
639 #endif
640  if (SP)
641  {
642  SP->AbsMultTranspose(local_diag, diag);
643  }
644 #ifdef MFEM_USE_MPI
645  else if (HP)
646  {
647  HP->AbsMultTranspose(1.0, local_diag, 0.0, diag);
648  }
649 #endif
650  else
651  {
652  MFEM_ABORT("Prolongation matrix has unexpected type.");
653  }
654  return;
655  }
656  if (!IsIdentityProlongation(P))
657  {
658  Vector local_diag(P->Height());
659  ext->AssembleDiagonal(local_diag);
660  P->MultTranspose(local_diag, diag);
661  }
662  else
663  {
664  ext->AssembleDiagonal(diag);
665  }
666  }
667  else
668  {
669  mat->GetDiag(diag);
670  }
671 }
672
673 void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
674  Vector &b, OperatorHandle &A, Vector &X,
675  Vector &B, int copy_interior)
676 {
677  if (ext)
678  {
679  ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
680  return;
681  }
683  FormSystemMatrix(ess_tdof_list, A);
684
685  // Transform the system and perform the elimination in B, based on the
686  // essential BC values from x. Restrict the BC part of x in X, and set the
687  // non-BC part to zero. Since there is no good initial guess for the Lagrange
688  // multipliers, set X = 0.0 for hybridization.
689  if (static_cond)
690  {
691  // Schur complement reduction to the exposed dofs
692  static_cond->ReduceSystem(x, b, X, B, copy_interior);
693  }
694  else if (!P) // conforming space
695  {
696  if (hybridization)
697  {
698  // Reduction to the Lagrange multipliers system
699  EliminateVDofsInRHS(ess_tdof_list, x, b);
700  hybridization->ReduceRHS(b, B);
701  X.SetSize(B.Size());
702  X = 0.0;
703  }
704  else
705  {
706  // A, X and B point to the same data as mat, x and b
707  EliminateVDofsInRHS(ess_tdof_list, x, b);
708  X.NewMemoryAndSize(x.GetMemory(), x.Size(), false);
709  B.NewMemoryAndSize(b.GetMemory(), b.Size(), false);
710  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
711  }
712  }
713  else // non-conforming space
714  {
715  if (hybridization)
716  {
717  // Reduction to the Lagrange multipliers system
719  Vector conf_b(P->Width()), conf_x(P->Width());
720  P->MultTranspose(b, conf_b);
721  R->Mult(x, conf_x);
722  EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
723  R->MultTranspose(conf_b, b); // store eliminated rhs in b
724  hybridization->ReduceRHS(conf_b, B);
725  X.SetSize(B.Size());
726  X = 0.0;
727  }
728  else
729  {
730  // Variational restriction with P
732  B.SetSize(P->Width());
733  P->MultTranspose(b, B);
734  X.SetSize(R->Height());
735  R->Mult(x, X);
736  EliminateVDofsInRHS(ess_tdof_list, X, B);
737  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
738  }
739  }
740 }
741
742 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
743  OperatorHandle &A)
744 {
745  if (ext)
746  {
747  ext->FormSystemMatrix(ess_tdof_list, A);
748  return;
749  }
750
751  // Finish the matrix assembly and perform BC elimination, storing the
752  // eliminated part of the matrix.
753  if (static_cond)
754  {
756  {
757  static_cond->SetEssentialTrueDofs(ess_tdof_list);
758  static_cond->Finalize(); // finalize Schur complement (to true dofs)
760  static_cond->Finalize(); // finalize eliminated part
761  }
762  A.Reset(&static_cond->GetMatrix(), false);
763  }
764  else
765  {
766  if (!mat_e)
767  {
769  if (P) { ConformingAssemble(); }
770  EliminateVDofs(ess_tdof_list, diag_policy);
771  const int remove_zeros = 0;
772  Finalize(remove_zeros);
773  }
774  if (hybridization)
775  {
776  A.Reset(&hybridization->GetMatrix(), false);
777  }
778  else
779  {
780  A.Reset(mat, false);
781  }
782  }
783 }
784
786  const Vector &b, Vector &x)
787 {
788  if (ext)
789  {
790  ext->RecoverFEMSolution(X, b, x);
791  return;
792  }
793
795  if (!P) // conforming space
796  {
797  if (static_cond)
798  {
799  // Private dofs back solve
800  static_cond->ComputeSolution(b, X, x);
801  }
802  else if (hybridization)
803  {
804  // Primal unknowns recovery
805  hybridization->ComputeSolution(b, X, x);
806  }
807  else
808  {
809  // X and x point to the same data
810
811  // If the validity flags of X's Memory were changed (e.g. if it was
812  // moved to device memory) then we need to tell x about that.
813  x.SyncMemory(X);
814  }
815  }
816  else // non-conforming space
817  {
818  if (static_cond)
819  {
820  // Private dofs back solve
821  static_cond->ComputeSolution(b, X, x);
822  }
823  else if (hybridization)
824  {
825  // Primal unknowns recovery
826  Vector conf_b(P->Width()), conf_x(P->Width());
827  P->MultTranspose(b, conf_b);
829  R->Mult(x, conf_x); // get essential b.c. from x
830  hybridization->ComputeSolution(conf_b, X, conf_x);
831  x.SetSize(P->Height());
832  P->Mult(conf_x, x);
833  }
834  else
835  {
836  // Apply conforming prolongation
837  x.SetSize(P->Height());
838  P->Mult(X, x);
839  }
840  }
841 }
842
844 {
845  if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
846  {
847  return;
848  }
849
850  int num_elements = fes->GetNE();
851  int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
852
853  element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
854  num_elements);
855
856  DenseMatrix tmp;
858
859 #ifdef MFEM_USE_LEGACY_OPENMP
860  #pragma omp parallel for private(tmp,eltrans)
861 #endif
862  for (int i = 0; i < num_elements; i++)
863  {
865  num_dofs_per_el, num_dofs_per_el);
866  const FiniteElement &fe = *fes->GetFE(i);
867 #ifdef MFEM_DEBUG
868  if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
869  mfem_error("BilinearForm::ComputeElementMatrices:"
870  " all elements must have same number of dofs");
871 #endif
872  fes->GetElementTransformation(i, &eltrans);
873
874  dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
875  for (int k = 1; k < dbfi.Size(); k++)
876  {
877  // note: some integrators may not be thread-safe
878  dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
879  elmat += tmp;
880  }
881  elmat.ClearExternalData();
882  }
883 }
884
885 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
886  const Vector &sol, Vector &rhs,
887  DiagonalPolicy dpolicy)
888 {
889  Array<int> ess_dofs, conf_ess_dofs;
890  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
891
892  if (fes->GetVSize() == height)
893  {
894  EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
895  }
896  else
897  {
898  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
899  EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
900  }
901 }
902
903 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
904  DiagonalPolicy dpolicy)
905 {
906  Array<int> ess_dofs, conf_ess_dofs;
907  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
908
909  if (fes->GetVSize() == height)
910  {
911  EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
912  }
913  else
914  {
915  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
916  EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
917  }
918 }
919
921  double value)
922 {
923  Array<int> ess_dofs, conf_ess_dofs;
924  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
925
926  if (fes->GetVSize() == height)
927  {
928  EliminateEssentialBCFromDofsDiag(ess_dofs, value);
929  }
930  else
931  {
932  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
933  EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
934  }
935 }
936
938  const Vector &sol, Vector &rhs,
939  DiagonalPolicy dpolicy)
940 {
941  for (int i = 0; i < vdofs.Size(); i++)
942  {
943  int vdof = vdofs[i];
944  if ( vdof >= 0 )
945  {
946  mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
947  }
948  else
949  {
950  mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
951  }
952  }
953 }
954
956  DiagonalPolicy dpolicy)
957 {
958  if (mat_e == NULL)
959  {
960  mat_e = new SparseMatrix(height);
961  }
962
963  for (int i = 0; i < vdofs.Size(); i++)
964  {
965  int vdof = vdofs[i];
966  if ( vdof >= 0 )
967  {
968  mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
969  }
970  else
971  {
972  mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
973  }
974  }
975 }
976
978  const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
979  DiagonalPolicy dpolicy)
980 {
981  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
982  MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
983  MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
984
985  for (int i = 0; i < ess_dofs.Size(); i++)
986  if (ess_dofs[i] < 0)
987  {
988  mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
989  }
990 }
991
993  DiagonalPolicy dpolicy)
994 {
995  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
996
997  for (int i = 0; i < ess_dofs.Size(); i++)
998  if (ess_dofs[i] < 0)
999  {
1000  mat -> EliminateRowCol (i, dpolicy);
1001  }
1002 }
1003
1005  double value)
1006 {
1007  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
1008
1009  for (int i = 0; i < ess_dofs.Size(); i++)
1010  if (ess_dofs[i] < 0)
1011  {
1012  mat -> EliminateRowColDiag (i, value);
1013  }
1014 }
1015
1017  const Array<int> &vdofs, const Vector &x, Vector &b)
1018 {
1020  mat->PartMult(vdofs, x, b);
1021 }
1022
1023 void BilinearForm::Mult(const Vector &x, Vector &y) const
1024 {
1025  if (ext)
1026  {
1027  ext->Mult(x, y);
1028  }
1029  else
1030  {
1031  mat->Mult(x, y);
1032  }
1033 }
1034
1036 {
1037  bool full_update;
1038
1039  if (nfes && nfes != fes)
1040  {
1041  full_update = true;
1042  fes = nfes;
1043  }
1044  else
1045  {
1046  // Check for different size (e.g. assembled form on non-conforming space)
1047  // or different sequence number.
1048  full_update = (fes->GetVSize() != Height() ||
1049  sequence < fes->GetSequence());
1050  }
1051
1052  delete mat_e;
1053  mat_e = NULL;
1055  delete static_cond;
1056  static_cond = NULL;
1057
1058  if (full_update)
1059  {
1060  delete mat;
1061  mat = NULL;
1062  delete hybridization;
1063  hybridization = NULL;
1064  sequence = fes->GetSequence();
1065  }
1066  else
1067  {
1068  if (mat) { *mat = 0.0; }
1069  if (hybridization) { hybridization->Reset(); }
1070  }
1071
1072  height = width = fes->GetVSize();
1073
1074  if (ext) { ext->Update(); }
1075 }
1076
1078 {
1079  diag_policy = policy;
1080 }
1081
1083 {
1084  delete mat_e;
1085  delete mat;
1086  delete element_matrices;
1087  delete static_cond;
1088  delete hybridization;
1089
1090  if (!extern_bfs)
1091  {
1092  int k;
1093  for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
1094  for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
1095  for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
1096  for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
1097  }
1098
1099  delete ext;
1100 }
1101
1102
1103 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1104  FiniteElementSpace *te_fes)
1105  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1106 {
1107  trial_fes = tr_fes;
1108  test_fes = te_fes;
1109  mat = NULL;
1110  mat_e = NULL;
1111  extern_bfs = 0;
1113  ext = NULL;
1114 }
1115
1116 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1117  FiniteElementSpace *te_fes,
1118  MixedBilinearForm * mbf)
1119  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1120 {
1121  trial_fes = tr_fes;
1122  test_fes = te_fes;
1123  mat = NULL;
1124  mat_e = NULL;
1125  extern_bfs = 1;
1126  ext = NULL;
1127
1128  // Copy the pointers to the integrators
1129  dbfi = mbf->dbfi;
1130  bbfi = mbf->bbfi;
1131  tfbfi = mbf->tfbfi;
1132  btfbfi = mbf->btfbfi;
1133
1134  bbfi_marker = mbf->bbfi_marker;
1136
1138  ext = NULL;
1139 }
1140
1142 {
1143  if (ext)
1144  {
1145  MFEM_ABORT("the assembly level has already been set!");
1146  }
1147  assembly = assembly_level;
1148  switch (assembly)
1149  {
1151  break;
1152  case AssemblyLevel::FULL:
1153  // ext = new FAMixedBilinearFormExtension(this);
1154  // Use the original BilinearForm implementation for now
1155  break;
1157  mfem_error("Element assembly not supported yet... stay tuned!");
1158  // ext = new EAMixedBilinearFormExtension(this);
1159  break;
1161  ext = new PAMixedBilinearFormExtension(this);
1162  break;
1163  case AssemblyLevel::NONE:
1164  mfem_error("Matrix-free action not supported yet... stay tuned!");
1165  // ext = new MFMixedBilinearFormExtension(this);
1166  break;
1167  default:
1168  mfem_error("Unknown assembly level");
1169  }
1170 }
1171
1172 double & MixedBilinearForm::Elem (int i, int j)
1173 {
1174  return (*mat)(i, j);
1175 }
1176
1177 const double & MixedBilinearForm::Elem (int i, int j) const
1178 {
1179  return (*mat)(i, j);
1180 }
1181
1182 void MixedBilinearForm::Mult(const Vector & x, Vector & y) const
1183 {
1184  y = 0.0;
1186 }
1187
1189  const double a) const
1190 {
1191  if (ext)
1192  {
1194  }
1195  else
1196  {
1198  }
1199 }
1200
1202 {
1203  y = 0.0;
1205 }
1206
1208  const double a) const
1209 {
1210  if (ext)
1211  {
1213  }
1214  else
1215  {
1217  }
1218 }
1219
1221 {
1223  {
1224  MFEM_WARNING("MixedBilinearForm::Inverse not possible with this assembly level!");
1225  return NULL;
1226  }
1227  else
1228  {
1229  return mat -> Inverse ();
1230  }
1231 }
1232
1233 void MixedBilinearForm::Finalize (int skip_zeros)
1234 {
1236  {
1237  mat -> Finalize (skip_zeros);
1238  }
1239 }
1240
1242 {
1243  MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
1245  "MixedBilinearForm::GetBlocks: both trial and test spaces "
1246  "must use Ordering::byNODES!");
1247
1248  blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
1249
1250  mat->GetBlocks(blocks);
1251 }
1252
1254 {
1255  dbfi.Append (bfi);
1256 }
1257
1259 {
1260  bbfi.Append (bfi);
1261  bbfi_marker.Append(NULL); // NULL marker means apply everywhere
1262 }
1263
1265  Array<int> &bdr_marker)
1266 {
1267  bbfi.Append (bfi);
1268  bbfi_marker.Append(&bdr_marker);
1269 }
1270
1272 {
1273  tfbfi.Append (bfi);
1274 }
1275
1277 {
1278  btfbfi.Append(bfi);
1279  btfbfi_marker.Append(NULL); // NULL marker means apply everywhere
1280 }
1281
1283  Array<int> &bdr_marker)
1284 {
1285  btfbfi.Append(bfi);
1286  btfbfi_marker.Append(&bdr_marker);
1287 }
1288
1289 void MixedBilinearForm::Assemble (int skip_zeros)
1290 {
1291  if (ext)
1292  {
1293  ext->Assemble();
1294  return;
1295  }
1296
1297  Array<int> tr_vdofs, te_vdofs;
1298  ElementTransformation *eltrans;
1300
1301  Mesh *mesh = test_fes -> GetMesh();
1302
1303  if (mat == NULL)
1304  {
1305  mat = new SparseMatrix(height, width);
1306  }
1307
1308  if (dbfi.Size())
1309  {
1310  for (int i = 0; i < test_fes -> GetNE(); i++)
1311  {
1312  trial_fes -> GetElementVDofs (i, tr_vdofs);
1313  test_fes -> GetElementVDofs (i, te_vdofs);
1314  eltrans = test_fes -> GetElementTransformation (i);
1315  for (int k = 0; k < dbfi.Size(); k++)
1316  {
1317  dbfi[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1318  *test_fes -> GetFE(i),
1319  *eltrans, elemmat);
1320  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1321  }
1322  }
1323  }
1324
1325  if (bbfi.Size())
1326  {
1327  // Which boundary attributes need to be processed?
1328  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1329  mesh->bdr_attributes.Max() : 0);
1330  bdr_attr_marker = 0;
1331  for (int k = 0; k < bbfi.Size(); k++)
1332  {
1333  if (bbfi_marker[k] == NULL)
1334  {
1335  bdr_attr_marker = 1;
1336  break;
1337  }
1338  Array<int> &bdr_marker = *bbfi_marker[k];
1339  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1340  "invalid boundary marker for boundary integrator #"
1341  << k << ", counting from zero");
1342  for (int i = 0; i < bdr_attr_marker.Size(); i++)
1343  {
1344  bdr_attr_marker[i] |= bdr_marker[i];
1345  }
1346  }
1347
1348  for (int i = 0; i < test_fes -> GetNBE(); i++)
1349  {
1350  const int bdr_attr = mesh->GetBdrAttribute(i);
1351  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1352
1353  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1354  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1355  eltrans = test_fes -> GetBdrElementTransformation (i);
1356  for (int k = 0; k < bbfi.Size(); k++)
1357  {
1358  if (bbfi_marker[k] &&
1359  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
1360
1361  bbfi[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1362  *test_fes -> GetBE(i),
1363  *eltrans, elemmat);
1364  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1365  }
1366  }
1367  }
1368
1369  if (tfbfi.Size())
1370  {
1372  Array<int> te_vdofs2;
1373  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1374
1375  int nfaces = mesh->GetNumFaces();
1376  for (int i = 0; i < nfaces; i++)
1377  {
1378  ftr = mesh->GetFaceElementTransformations(i);
1379  trial_fes->GetFaceVDofs(i, tr_vdofs);
1380  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1381  trial_face_fe = trial_fes->GetFaceElement(i);
1382  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1383  if (ftr->Elem2No >= 0)
1384  {
1385  test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
1386  te_vdofs.Append(te_vdofs2);
1387  test_fe2 = test_fes->GetFE(ftr->Elem2No);
1388  }
1389  else
1390  {
1391  // The test_fe2 object is really a dummy and not used on the
1392  // boundaries, but we can't dereference a NULL pointer, and we don't
1393  // want to actually make a fake element.
1394  test_fe2 = test_fe1;
1395  }
1396  for (int k = 0; k < tfbfi.Size(); k++)
1397  {
1398  tfbfi[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1399  *ftr, elemmat);
1401  }
1402  }
1403  }
1404
1405  if (btfbfi.Size())
1406  {
1408  Array<int> te_vdofs2;
1409  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1410
1411  // Which boundary attributes need to be processed?
1412  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1413  mesh->bdr_attributes.Max() : 0);
1414  bdr_attr_marker = 0;
1415  for (int k = 0; k < btfbfi.Size(); k++)
1416  {
1417  if (btfbfi_marker[k] == NULL)
1418  {
1419  bdr_attr_marker = 1;
1420  break;
1421  }
1422  Array<int> &bdr_marker = *btfbfi_marker[k];
1423  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1424  "invalid boundary marker for boundary trace face integrator #"
1425  << k << ", counting from zero");
1426  for (int i = 0; i < bdr_attr_marker.Size(); i++)
1427  {
1428  bdr_attr_marker[i] |= bdr_marker[i];
1429  }
1430  }
1431
1432  for (int i = 0; i < trial_fes -> GetNBE(); i++)
1433  {
1434  const int bdr_attr = mesh->GetBdrAttribute(i);
1435  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1436
1437  ftr = mesh->GetBdrFaceTransformations(i);
1438  if (ftr)
1439  {
1440  trial_fes->GetFaceVDofs(i, tr_vdofs);
1441  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1442  trial_face_fe = trial_fes->GetFaceElement(i);
1443  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1444  // The test_fe2 object is really a dummy and not used on the
1445  // boundaries, but we can't dereference a NULL pointer, and we don't
1446  // want to actually make a fake element.
1447  test_fe2 = test_fe1;
1448  for (int k = 0; k < btfbfi.Size(); k++)
1449  {
1450  if (btfbfi_marker[k] &&
1451  (*btfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
1452
1453  btfbfi[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1454  *ftr, elemmat);
1456  }
1457  }
1458  }
1459  }
1460 }
1461
1463  Vector &diag) const
1464 {
1465  if (ext)
1466  {
1467  MFEM_ASSERT(diag.Size() == test_fes->GetTrueVSize(),
1468  "Vector for holding diagonal has wrong size!");
1469  MFEM_ASSERT(D.Size() == trial_fes->GetTrueVSize(),
1470  "Vector for holding diagonal has wrong size!");
1471  const Operator *P_trial = trial_fes->GetProlongationMatrix();
1472  const Operator *P_test = test_fes->GetProlongationMatrix();
1473  if (!IsIdentityProlongation(P_trial))
1474  {
1475  Vector local_D(P_trial->Height());
1476  P_trial->Mult(D, local_D);
1477
1478  if (!IsIdentityProlongation(P_test))
1479  {
1480  Vector local_diag(P_test->Height());
1482  P_test->MultTranspose(local_diag, diag);
1483  }
1484  else
1485  {
1487  }
1488  }
1489  else
1490  {
1491  if (!IsIdentityProlongation(P_test))
1492  {
1493  Vector local_diag(P_test->Height());
1495  P_test->MultTranspose(local_diag, diag);
1496  }
1497  else
1498  {
1500  }
1501  }
1502  }
1503  else
1504  {
1505  MFEM_ABORT("Not implemented. Maybe assemble your bilinear form into a "
1506  "matrix and use SparseMatrix functions?");
1507  }
1508 }
1509
1511 {
1513  {
1514  MFEM_WARNING("Conforming assemble not supported for this assembly level!");
1515  return;
1516  }
1517
1518  Finalize();
1519
1521  if (P2)
1522  {
1523  SparseMatrix *R = Transpose(*P2);
1524  SparseMatrix *RA = mfem::Mult(*R, *mat);
1525  delete R;
1526  delete mat;
1527  mat = RA;
1528  }
1529
1531  if (P1)
1532  {
1533  SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1534  delete mat;
1535  mat = RAP;
1536  }
1537
1538  height = mat->Height();
1539  width = mat->Width();
1540 }
1541
1542
1544 {
1545  if (dbfi.Size())
1546  {
1547  const FiniteElement &trial_fe = *trial_fes->GetFE(i);
1548  const FiniteElement &test_fe = *test_fes->GetFE(i);
1550  dbfi[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans, elmat);
1551  for (int k = 1; k < dbfi.Size(); k++)
1552  {
1553  dbfi[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans, elemmat);
1554  elmat += elemmat;
1555  }
1556  }
1557  else
1558  {
1561  elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
1562  elmat = 0.0;
1563  }
1564 }
1565
1567 {
1568  if (bbfi.Size())
1569  {
1570  const FiniteElement &trial_be = *trial_fes->GetBE(i);
1571  const FiniteElement &test_be = *test_fes->GetBE(i);
1573  bbfi[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans, elmat);
1574  for (int k = 1; k < bbfi.Size(); k++)
1575  {
1576  bbfi[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans, elemmat);
1577  elmat += elemmat;
1578  }
1579  }
1580  else
1581  {
1584  elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
1585  elmat = 0.0;
1586  }
1587 }
1588
1590  int i, const DenseMatrix &elmat, int skip_zeros)
1591 {
1592  AssembleElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1593 }
1594
1596  int i, const DenseMatrix &elmat, Array<int> &trial_vdofs,
1597  Array<int> &test_vdofs, int skip_zeros)
1598 {
1599  trial_fes->GetElementVDofs(i, trial_vdofs);
1600  test_fes->GetElementVDofs(i, test_vdofs);
1601  if (mat == NULL)
1602  {
1603  mat = new SparseMatrix(height, width);
1604  }
1606 }
1607
1609  int i, const DenseMatrix &elmat, int skip_zeros)
1610 {
1611  AssembleBdrElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1612 }
1613
1615  int i, const DenseMatrix &elmat, Array<int> &trial_vdofs,
1616  Array<int> &test_vdofs, int skip_zeros)
1617 {
1618  trial_fes->GetBdrElementVDofs(i, trial_vdofs);
1619  test_fes->GetBdrElementVDofs(i, test_vdofs);
1620  if (mat == NULL)
1621  {
1622  mat = new SparseMatrix(height, width);
1623  }
1625 }
1626
1628  const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
1629 {
1630  int i, j, k;
1631  Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1632
1633  cols_marker = 0;
1634  for (i = 0; i < trial_fes -> GetNBE(); i++)
1635  if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1636  {
1637  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1638  for (j = 0; j < tr_vdofs.Size(); j++)
1639  {
1640  if ( (k = tr_vdofs[j]) < 0 )
1641  {
1642  k = -1-k;
1643  }
1644  cols_marker[k] = 1;
1645  }
1646  }
1647  mat -> EliminateCols (cols_marker, &sol, &rhs);
1648 }
1649
1651  const Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
1652 {
1653  mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1654 }
1655
1657 {
1658  int i, j, k;
1659  Array<int> te_vdofs;
1660
1661  for (i = 0; i < test_fes -> GetNBE(); i++)
1662  if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1663  {
1664  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1665  for (j = 0; j < te_vdofs.Size(); j++)
1666  {
1667  if ( (k = te_vdofs[j]) < 0 )
1668  {
1669  k = -1-k;
1670  }
1671  mat -> EliminateRow (k);
1672  }
1673  }
1674 }
1675
1677  &trial_tdof_list,
1678  const Array<int> &test_tdof_list,
1679  OperatorHandle &A)
1680
1681 {
1682  if (ext)
1683  {
1684  ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
1685  return;
1686  }
1687
1688  const SparseMatrix *test_P = test_fes->GetConformingProlongation();
1689  const SparseMatrix *trial_P = trial_fes->GetConformingProlongation();
1690
1691  mat->Finalize();
1692
1693  if (test_P) // TODO: Must actually check for trial_P too
1694  {
1695  SparseMatrix *m = RAP(*test_P, *mat, *trial_P);
1696  delete mat;
1697  mat = m;
1698  }
1699
1700  Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
1702  ess_trial_tdof_marker);
1704  ess_test_tdof_marker);
1705
1706  mat_e = new SparseMatrix(mat->Height(), mat->Width());
1707  mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
1708
1709  for (int i=0; i<test_tdof_list.Size(); ++i)
1710  {
1711  mat->EliminateRow(test_tdof_list[i]);
1712  }
1713  mat_e->Finalize();
1714  A.Reset(mat, false);
1715 }
1716
1718  &trial_tdof_list,
1719  const Array<int> &test_tdof_list,
1720  Vector &x, Vector &b,
1721  OperatorHandle &A,
1722  Vector &X, Vector &B)
1723 {
1724  if (ext)
1725  {
1726  ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list, x, b, A, X,
1727  B);
1728  return;
1729  }
1730
1731  const Operator *Pi = this->GetProlongation();
1732  const Operator *Po = this->GetOutputProlongation();
1733  const Operator *Ri = this->GetRestriction();
1734  InitTVectors(Po, Ri, Pi, x, b, X, B);
1735
1736  if (!mat_e)
1737  {
1738  FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list,
1739  A); // Set A = mat_e
1740  }
1741  // Eliminate essential BCs with B -= Ab xb
1743
1744  B.SetSubVector(test_tdof_list, 0.0);
1745 }
1746
1748 {
1749  delete mat;
1750  mat = NULL;
1751  delete mat_e;
1752  mat_e = NULL;
1753  height = test_fes->GetVSize();
1754  width = trial_fes->GetVSize();
1755  if (ext) { ext->Update(); }
1756 }
1757
1759 {
1760  if (mat) { delete mat; }
1761  if (mat_e) { delete mat_e; }
1762  if (!extern_bfs)
1763  {
1764  int i;
1765  for (i = 0; i < dbfi.Size(); i++) { delete dbfi[i]; }
1766  for (i = 0; i < bbfi.Size(); i++) { delete bbfi[i]; }
1767  for (i = 0; i < tfbfi.Size(); i++) { delete tfbfi[i]; }
1768  for (i = 0; i < btfbfi.Size(); i++) { delete btfbfi[i]; }
1769  }
1770  delete ext;
1771 }
1772
1773
1775 {
1776  Array<int> dom_vdofs, ran_vdofs;
1778  const FiniteElement *dom_fe, *ran_fe;
1779  DenseMatrix totelmat, elmat;
1780
1781  if (mat == NULL)
1782  {
1783  mat = new SparseMatrix(height, width);
1784  }
1785
1786  if (dbfi.Size() > 0)
1787  {
1788  for (int i = 0; i < test_fes->GetNE(); i++)
1789  {
1790  trial_fes->GetElementVDofs(i, dom_vdofs);
1791  test_fes->GetElementVDofs(i, ran_vdofs);
1793  dom_fe = trial_fes->GetFE(i);
1794  ran_fe = test_fes->GetFE(i);
1795
1796  dbfi[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1797  for (int j = 1; j < dbfi.Size(); j++)
1798  {
1799  dbfi[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1800  totelmat += elmat;
1801  }
1802  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1803  }
1804  }
1805
1806  if (tfbfi.Size())
1807  {
1808  const int nfaces = test_fes->GetMesh()->GetNumFaces();
1809  for (int i = 0; i < nfaces; i++)
1810  {
1811  trial_fes->GetFaceVDofs(i, dom_vdofs);
1812  test_fes->GetFaceVDofs(i, ran_vdofs);
1814  dom_fe = trial_fes->GetFaceElement(i);
1815  ran_fe = test_fes->GetFaceElement(i);
1816
1817  tfbfi[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1818  for (int j = 1; j < tfbfi.Size(); j++)
1819  {
1820  tfbfi[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1821  totelmat += elmat;
1822  }
1823  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1824  }
1825  }
1826 }
1827
1828 }
Abstract class for all finite elements.
Definition: fe.hpp:235
void EliminateEssentialBCFromDofs(const Array< int > &ess_dofs, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateVDofs(const Array&lt;int&gt; &amp;, const Vector &amp;, Vector &amp;, DiagonalPolicy) but here ess_...
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:412
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:521
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:400
int * GetJ()
Definition: table.hpp:114
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1053
Data and methods for matrix-free bilinear forms.
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:874
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse: .
void SetConstraintIntegrator(BilinearFormIntegrator *c_integ)
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
Definition: sparsemat.hpp:468
Array< BilinearFormIntegrator * > fbfi
Set of interior face Integrators to be applied.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
bool ReducesTrueVSize() const
Definition: staticcond.cpp:118
virtual const Operator * GetProlongation() const
Get the input finite element space prolongation matrix.
virtual void EliminateTestDofs(const Array< int > &bdr_attr_is_ess)
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1640
void EliminateTrialDofs(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Return in A that is column-constrained.
void SyncMemory(const Vector &v)
Update the memory location of the vector to match v.
Definition: vector.hpp:190
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:460
BilinearFormExtension * ext
Extension for supporting Full Assembly (FA), Element Assembly (EA), Partial Assembly (PA)...
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition: table.cpp:200
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
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
Array< BilinearFormIntegrator * > dbfi
Set of Domain Integrators to be applied.
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:831
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:71
Array< Array< int > * > bbfi_marker
Entries are not owned.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:178
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
Definition: fespace.cpp:346
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:601
int batch
Element batch size used in the form action (1, 8, num_elems, etc.)
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
void ReduceSystem(Vector &x, Vector &b, Vector &X, Vector &B, int copy_interior=0) const
Set the reduced solution X and r.h.s B vectors from the full linear system solution x and r...
Definition: staticcond.cpp:418
virtual void AddMult(const Vector &x, Vector &y, const double c=1.0) const =0
int * GetI()
Return the array I.
Definition: sparsemat.hpp:173
AssemblyLevel assembly
The form assembly level (full, partial, etc.)
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
y += At * x (default) or y += a * At * x
Definition: sparsemat.cpp:728
void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag.
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Definition: staticcond.cpp:288
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:154
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
Array< BilinearFormIntegrator * > btfbfi
Boundary trace face (skeleton) integrators.
Array< BilinearFormIntegrator * > bbfi
Boundary integrators.
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th face in the ...
Definition: fespace.cpp:2073
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
Definition: staticcond.cpp:478
virtual const Operator * GetRestriction() const
Get the input finite element space restriction matrix.
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:869
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:183
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
virtual void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const =0
Data and methods for partially-assembled bilinear forms.
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Definition: sparsemat.cpp:1208
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:92
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)=0
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
Compute the boundary element matrix of the given boundary element.
void EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
Enable hybridization.
Adds a boundary trace face integrator. Assumes ownership of bfi.
void AssembleElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given element matrix.
void EliminateEssentialBCFromDofsDiag(const Array< int > &ess_dofs, double value)
Perform elimination and set the diagonal entry to the given value.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:4184
long sequence
Indicates the Mesh::sequence corresponding to the current state of the BilinearForm.
void AbsMultTranspose(double a, const Vector &x, double b, Vector &y) const
Computes y = a * |At| * x + b * y, using entry-wise absolute values of the transpose of matrix A...
Definition: hypre.cpp:1066
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
Adds a domain integrator. Assumes ownership of bfi.
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:881
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:427
void SetSize(int m, int n)
Definition: array.hpp:352
void Finalize()
Finalize the construction of the hybridized matrix.
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
virtual void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)=0
Data type sparse matrix.
Definition: sparsemat.hpp:46
virtual void AssembleDiagonal(Vector &diag) const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:65
SparseMatrix & GetMatrix()
Return the serial hybridized matrix.
StaticCondensation * static_cond
Owned.
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
SparseMatrix * mat
Sparse matrix to be associated with the form. Owned.
void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)
Form the linear system A X = B, corresponding to this mixed bilinear form and the linear form b(...
Data and methods for fully-assembled bilinear forms.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
double b
Definition: lissajous.cpp:42
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
Returns the transformation defining the given face element in a user-defined variable.
Definition: mesh.cpp:480
void EliminateVDofs(const Array< int > &vdofs, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the given vdofs. NOTE: here, vdofs is a list of DOFs.
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
Abstract data type matrix.
Definition: matrix.hpp:27
int extern_bfs
Indicates the BilinearFormIntegrators stored in dbfi, bbfi, fbfi, and bfbfi are owned by another Bili...
void AssembleElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given element matrix.
void EliminateEssentialBC(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate essential boundary DOFs from the system.
virtual void Assemble()=0
Assemble at the level given for the BilinearFormExtension subclass.
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:330
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:414
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition: fespace.hpp:469
double * GetData(int k)
Definition: densemat.hpp:816
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:403
bool Conforming() const
Definition: fespace.hpp:320
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T into diag, where A is this mixed bilinear form and D is a diagonal...
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:389
SparseMatrix * mat
Owned.
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1003
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
Add a trace face integrator. Assumes ownership of bfi.
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2436
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Dynamic 2D array using row-major layout.
Definition: array.hpp:333
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Definition: vector.cpp:656
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2512
int SizeI() const
Definition: densemat.hpp:765
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:466
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACYFULL.
SparseMatrix * mat_e
Sparse Matrix used to store the eliminations from the b.c. Owned. .
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:460
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:672
Data and methods for element-assembled bilinear forms.
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:31
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
Array< Array< int > * > btfbfi_marker
Entries are not owned.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:315
MixedBilinearFormExtension * ext
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:594
Array< BilinearFormIntegrator * > bfbfi
Set of boundary face Integrators to be applied.
Array< BilinearFormIntegrator * > bbfi
Set of Boundary Integrators to be applied.
DenseTensor * element_matrices
Owned.
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:235
void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi, Vector &x, Vector &b, Vector &X, Vector &B) const
Initializes memory for true vectors of linear system.
Definition: operator.cpp:22
void ClearExternalData()
Definition: densemat.hpp:80
virtual double & Elem(int i, int j)
Returns a reference to: .
DiagonalPolicy diag_policy
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:227
Array< Array< int > * > bbfi_marker
Entries are not owned.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition: operator.cpp:85
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1463
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Table * GetFaceToElementTable() const
Definition: mesh.cpp:4903
FiniteElementSpace * test_fes
Not owned.
int SizeJ() const
Definition: densemat.hpp:766
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:530
int extern_bfs
Indicates the BilinearFormIntegrators stored in dbfi, bbfi, tfbfi and btfbfi are owned by another Mix...
virtual const Operator * GetOutputProlongation() const
Get the test finite element space prolongation matrix.
A standard isoparametric element transformation.
Definition: eltrans.hpp:348
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AssembleMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:189
SparseMatrix * mat_e
Owned.
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
Use the stored eliminated part of the matrix (see EliminateVDofs(const Array&lt;int&gt; &amp;...
Keep the diagonal value.
Definition: operator.hpp:49
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
void AbsMultTranspose(const Vector &x, Vector &y) const
y = |At| * x, using entry-wise absolute values of the transpose of matrix A
Definition: sparsemat.cpp:933
Adds new Boundary Integrator. Assumes ownership of bfi.
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:45
Adds new interior Face Integrator. Assumes ownership of bfi.
void NewMemoryAndSize(const Memory< double > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
Definition: vector.hpp:508
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:721
FiniteElementSpace * fes
FE space on which the form lives. Not owned.
Adds new Domain Integrator. Assumes ownership of bfi.
void ComputeElementMatrices()
Compute and store internally all element matrices.
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
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
Definition: sparsemat.cpp:1551
DenseMatrix elemmat
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th face element (2D and 3D).
Definition: fespace.cpp:185
FiniteElementSpace * trial_fes
Not owned.
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
virtual double & Elem(int i, int j)
Returns a reference to: .
virtual ~BilinearForm()
Destroys bilinear form.
Data and methods for partially-assembled mixed bilinear forms.
Vector data type.
Definition: vector.hpp:51
Adds a boundary integrator. Assumes ownership of bfi.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
Definition: staticcond.hpp:152
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)=0
virtual MatrixInverse * Inverse() const
Returns a pointer to (an approximation) of the matrix inverse.
Hybridization * hybridization
Owned.
void ReduceRHS(const Vector &b, Vector &b_r) const
int * GetI()
Definition: table.hpp:113
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
Definition: fespace.hpp:543
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
Compute the boundary element matrix of the given boundary element.
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
Array< Array< int > * > bfbfi_marker
Entries are not owned.
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:2047
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:179
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
Definition: staticcond.hpp:128
Array< int > vdofs
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
Abstract operator.
Definition: operator.hpp:24
void FreeElementMatrices()
Free the memory used by the element matrices.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:334
Array< BilinearFormIntegrator * > dbfi
Domain integrators.
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:707
AssemblyLevel assembly
The assembly level of the form (full, partial, etc.)
void Init(bool symmetric, bool block_diagonal)
Definition: staticcond.cpp:134
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:728
void UseSparsity(int *I, int *J, bool isSorted)
Use the given CSR sparsity pattern to allocate the internal SparseMatrix.
Array< BilinearFormIntegrator * > tfbfi
Trace face (skeleton) integrators.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:137