MFEM  v4.1.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  {
124  case AssemblyLevel::FULL:
125  // ext = new FABilinearFormExtension(this);
126  // Use the original BilinearForm implementation for now
127  break;
129  mfem_error("Element assembly not supported yet... stay tuned!");
130  // ext = new EABilinearFormExtension(this);
131  break;
133  ext = new PABilinearFormExtension(this);
134  break;
135  case AssemblyLevel::NONE:
136  mfem_error("Matrix-free action not supported yet... stay tuned!");
137  // ext = new MFBilinearFormExtension(this);
138  break;
139  default:
140  mfem_error("Unknown assembly level");
141  }
142 }
143
145 {
146  delete static_cond;
148  {
149  static_cond = NULL;
150  MFEM_WARNING("Static condensation not supported for this assembly level");
151  return;
152  }
155  {
156  bool symmetric = false; // TODO
157  bool block_diagonal = false; // TODO
158  static_cond->Init(symmetric, block_diagonal);
159  }
160  else
161  {
162  delete static_cond;
163  static_cond = NULL;
164  }
165 }
166
168  BilinearFormIntegrator *constr_integ,
169  const Array<int> &ess_tdof_list)
170 {
171  delete hybridization;
173  {
174  delete constr_integ;
175  hybridization = NULL;
176  MFEM_WARNING("Hybridization not supported for this assembly level");
177  return;
178  }
179  hybridization = new Hybridization(fes, constr_space);
180  hybridization->SetConstraintIntegrator(constr_integ);
181  hybridization->Init(ess_tdof_list);
182 }
183
184 void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
185 {
186  if (static_cond) { return; }
187
188  if (mat)
189  {
190  if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
191  {
192  return; // mat is already using the given sparsity
193  }
194  delete mat;
195  }
196  height = width = fes->GetVSize();
197  mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
198 }
199
201 {
202  MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
203  "invalid matrix A dimensions: "
204  << A.Height() << " x " << A.Width());
205  MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
206
207  UseSparsity(A.GetI(), A.GetJ(), A.ColumnsAreSorted());
208 }
209
210 double& BilinearForm::Elem (int i, int j)
211 {
212  return mat -> Elem(i,j);
213 }
214
215 const double& BilinearForm::Elem (int i, int j) const
216 {
217  return mat -> Elem(i,j);
218 }
219
221 {
222  return mat -> Inverse();
223 }
224
225 void BilinearForm::Finalize (int skip_zeros)
226 {
228  {
229  if (!static_cond) { mat->Finalize(skip_zeros); }
230  if (mat_e) { mat_e->Finalize(skip_zeros); }
231  if (static_cond) { static_cond->Finalize(); }
233  }
234 }
235
237 {
238  dbfi.Append(bfi);
239 }
240
242 {
243  bbfi.Append (bfi);
244  bbfi_marker.Append(NULL); // NULL marker means apply everywhere
245 }
246
248  Array<int> &bdr_marker)
249 {
250  bbfi.Append (bfi);
251  bbfi_marker.Append(&bdr_marker);
252 }
253
255 {
256  fbfi.Append (bfi);
257 }
258
260 {
261  bfbfi.Append(bfi);
262  bfbfi_marker.Append(NULL); // NULL marker means apply everywhere
263 }
264
266  Array<int> &bdr_marker)
267 {
268  bfbfi.Append(bfi);
269  bfbfi_marker.Append(&bdr_marker);
270 }
271
273 {
274  if (element_matrices)
275  {
277  elmat = element_matrices->GetData(i);
278  return;
279  }
280
281  if (dbfi.Size())
282  {
283  const FiniteElement &fe = *fes->GetFE(i);
285  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
286  for (int k = 1; k < dbfi.Size(); k++)
287  {
288  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
289  elmat += elemmat;
290  }
291  }
292  else
293  {
295  elmat.SetSize(vdofs.Size());
296  elmat = 0.0;
297  }
298 }
299
301 {
302  if (bbfi.Size())
303  {
304  const FiniteElement &be = *fes->GetBE(i);
306  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
307  for (int k = 1; k < bbfi.Size(); k++)
308  {
309  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
310  elmat += elemmat;
311  }
312  }
313  else
314  {
316  elmat.SetSize(vdofs.Size());
317  elmat = 0.0;
318  }
319 }
320
322  int i, const DenseMatrix &elmat, int skip_zeros)
323 {
324  AssembleElementMatrix(i, elmat, vdofs, skip_zeros);
325 }
326
328  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
329 {
330  fes->GetElementVDofs(i, vdofs);
331  if (static_cond)
332  {
333  static_cond->AssembleMatrix(i, elmat);
334  }
335  else
336  {
337  if (mat == NULL)
338  {
339  AllocMat();
340  }
342  if (hybridization)
343  {
344  hybridization->AssembleMatrix(i, elmat);
345  }
346  }
347 }
348
350  int i, const DenseMatrix &elmat, int skip_zeros)
351 {
352  AssembleBdrElementMatrix(i, elmat, vdofs, skip_zeros);
353 }
354
356  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
357 {
358  fes->GetBdrElementVDofs(i, vdofs);
359  if (static_cond)
360  {
361  static_cond->AssembleBdrMatrix(i, elmat);
362  }
363  else
364  {
365  if (mat == NULL)
366  {
367  AllocMat();
368  }
370  if (hybridization)
371  {
372  hybridization->AssembleBdrMatrix(i, elmat);
373  }
374  }
375 }
376
377 void BilinearForm::Assemble(int skip_zeros)
378 {
379  if (ext)
380  {
381  ext->Assemble();
382  return;
383  }
384
385  ElementTransformation *eltrans;
386  Mesh *mesh = fes -> GetMesh();
387  DenseMatrix elmat, *elmat_p;
388
389  if (mat == NULL)
390  {
391  AllocMat();
392  }
393
394 #ifdef MFEM_USE_LEGACY_OPENMP
395  int free_element_matrices = 0;
396  if (!element_matrices)
397  {
399  free_element_matrices = 1;
400  }
401 #endif
402
403  if (dbfi.Size())
404  {
405  for (int i = 0; i < fes -> GetNE(); i++)
406  {
408  if (element_matrices)
409  {
410  elmat_p = &(*element_matrices)(i);
411  }
412  else
413  {
414  const FiniteElement &fe = *fes->GetFE(i);
415  eltrans = fes->GetElementTransformation(i);
416  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
417  for (int k = 1; k < dbfi.Size(); k++)
418  {
419  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
420  elmat += elemmat;
421  }
422  elmat_p = &elmat;
423  }
424  if (static_cond)
425  {
426  static_cond->AssembleMatrix(i, *elmat_p);
427  }
428  else
429  {
431  if (hybridization)
432  {
433  hybridization->AssembleMatrix(i, *elmat_p);
434  }
435  }
436  }
437  }
438
439  if (bbfi.Size())
440  {
441  // Which boundary attributes need to be processed?
442  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
443  mesh->bdr_attributes.Max() : 0);
444  bdr_attr_marker = 0;
445  for (int k = 0; k < bbfi.Size(); k++)
446  {
447  if (bbfi_marker[k] == NULL)
448  {
449  bdr_attr_marker = 1;
450  break;
451  }
452  Array<int> &bdr_marker = *bbfi_marker[k];
453  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
454  "invalid boundary marker for boundary integrator #"
455  << k << ", counting from zero");
456  for (int i = 0; i < bdr_attr_marker.Size(); i++)
457  {
458  bdr_attr_marker[i] |= bdr_marker[i];
459  }
460  }
461
462  for (int i = 0; i < fes -> GetNBE(); i++)
463  {
464  const int bdr_attr = mesh->GetBdrAttribute(i);
465  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
466
467  const FiniteElement &be = *fes->GetBE(i);
468  fes -> GetBdrElementVDofs (i, vdofs);
469  eltrans = fes -> GetBdrElementTransformation (i);
470  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
471  for (int k = 1; k < bbfi.Size(); k++)
472  {
473  if (bbfi_marker[k] &&
474  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
475
476  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
477  elmat += elemmat;
478  }
479  if (!static_cond)
480  {
482  if (hybridization)
483  {
484  hybridization->AssembleBdrMatrix(i, elmat);
485  }
486  }
487  else
488  {
489  static_cond->AssembleBdrMatrix(i, elmat);
490  }
491  }
492  }
493
494  if (fbfi.Size())
495  {
497  Array<int> vdofs2;
498
499  int nfaces = mesh->GetNumFaces();
500  for (int i = 0; i < nfaces; i++)
501  {
502  tr = mesh -> GetInteriorFaceTransformations (i);
503  if (tr != NULL)
504  {
505  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
506  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
507  vdofs.Append (vdofs2);
508  for (int k = 0; k < fbfi.Size(); k++)
509  {
510  fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
511  *fes -> GetFE (tr -> Elem2No),
512  *tr, elemmat);
513  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
514  }
515  }
516  }
517  }
518
519  if (bfbfi.Size())
520  {
522  const FiniteElement *fe1, *fe2;
523
524  // Which boundary attributes need to be processed?
525  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
526  mesh->bdr_attributes.Max() : 0);
527  bdr_attr_marker = 0;
528  for (int k = 0; k < bfbfi.Size(); k++)
529  {
530  if (bfbfi_marker[k] == NULL)
531  {
532  bdr_attr_marker = 1;
533  break;
534  }
535  Array<int> &bdr_marker = *bfbfi_marker[k];
536  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
537  "invalid boundary marker for boundary face integrator #"
538  << k << ", counting from zero");
539  for (int i = 0; i < bdr_attr_marker.Size(); i++)
540  {
541  bdr_attr_marker[i] |= bdr_marker[i];
542  }
543  }
544
545  for (int i = 0; i < fes -> GetNBE(); i++)
546  {
547  const int bdr_attr = mesh->GetBdrAttribute(i);
548  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
549
550  tr = mesh -> GetBdrFaceTransformations (i);
551  if (tr != NULL)
552  {
553  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
554  fe1 = fes -> GetFE (tr -> Elem1No);
555  // The fe2 object is really a dummy and not used on the boundaries,
556  // but we can't dereference a NULL pointer, and we don't want to
557  // actually make a fake element.
558  fe2 = fe1;
559  for (int k = 0; k < bfbfi.Size(); k++)
560  {
561  if (bfbfi_marker[k] &&
562  (*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
563
564  bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
565  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
566  }
567  }
568  }
569  }
570
571 #ifdef MFEM_USE_LEGACY_OPENMP
572  if (free_element_matrices)
573  {
575  }
576 #endif
577 }
578
580 {
581  // Do not remove zero entries to preserve the symmetric structure of the
582  // matrix which in turn will give rise to symmetric structure in the new
583  // matrix. This ensures that subsequent calls to EliminateRowCol will work
584  // correctly.
585  Finalize(0);
586  MFEM_ASSERT(mat, "the BilinearForm is not assembled");
587
589  if (!P) { return; } // conforming mesh
590
591  SparseMatrix *R = Transpose(*P);
592  SparseMatrix *RA = mfem::Mult(*R, *mat);
593  delete mat;
594  if (mat_e)
595  {
596  SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
597  delete mat_e;
598  mat_e = RAe;
599  }
600  delete R;
601  mat = mfem::Mult(*RA, *P);
602  delete RA;
603  if (mat_e)
604  {
605  SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
606  delete mat_e;
607  mat_e = RAeP;
608  }
609
610  height = mat->Height();
611  width = mat->Width();
612 }
613
615 {
616  if (ext)
617  {
618  MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
619  "Vector for holding diagonal has wrong size!");
620  const Operator *P = fes->GetProlongationMatrix();
621  if (!IsIdentityProlongation(P))
622  {
623  Vector local_diag(P->Height());
624  ext->AssembleDiagonal(local_diag);
625  P->MultTranspose(local_diag, diag);
626  }
627  else
628  {
629  ext->AssembleDiagonal(diag);
630  }
631  }
632  else
633  {
634  MFEM_ABORT("Not implemented. Maybe assemble your bilinear form into a "
635  "matrix and use SparseMatrix::GetDiag?");
636  }
637 }
638
639 void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
640  Vector &b, OperatorHandle &A, Vector &X,
641  Vector &B, int copy_interior)
642 {
643  if (ext)
644  {
645  ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
646  return;
647  }
649  FormSystemMatrix(ess_tdof_list, A);
650
651  // Transform the system and perform the elimination in B, based on the
652  // essential BC values from x. Restrict the BC part of x in X, and set the
653  // non-BC part to zero. Since there is no good initial guess for the Lagrange
654  // multipliers, set X = 0.0 for hybridization.
655  if (static_cond)
656  {
657  // Schur complement reduction to the exposed dofs
658  static_cond->ReduceSystem(x, b, X, B, copy_interior);
659  }
660  else if (!P) // conforming space
661  {
662  if (hybridization)
663  {
664  // Reduction to the Lagrange multipliers system
665  EliminateVDofsInRHS(ess_tdof_list, x, b);
666  hybridization->ReduceRHS(b, B);
667  X.SetSize(B.Size());
668  X = 0.0;
669  }
670  else
671  {
672  // A, X and B point to the same data as mat, x and b
673  EliminateVDofsInRHS(ess_tdof_list, x, b);
674  X.NewMemoryAndSize(x.GetMemory(), x.Size(), false);
675  B.NewMemoryAndSize(b.GetMemory(), b.Size(), false);
676  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
677  }
678  }
679  else // non-conforming space
680  {
681  if (hybridization)
682  {
683  // Reduction to the Lagrange multipliers system
685  Vector conf_b(P->Width()), conf_x(P->Width());
686  P->MultTranspose(b, conf_b);
687  R->Mult(x, conf_x);
688  EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
689  R->MultTranspose(conf_b, b); // store eliminated rhs in b
690  hybridization->ReduceRHS(conf_b, B);
691  X.SetSize(B.Size());
692  X = 0.0;
693  }
694  else
695  {
696  // Variational restriction with P
698  B.SetSize(P->Width());
699  P->MultTranspose(b, B);
700  X.SetSize(R->Height());
701  R->Mult(x, X);
702  EliminateVDofsInRHS(ess_tdof_list, X, B);
703  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
704  }
705  }
706 }
707
708 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
709  OperatorHandle &A)
710 {
711  if (ext)
712  {
713  ext->FormSystemMatrix(ess_tdof_list, A);
714  return;
715  }
716
717  // Finish the matrix assembly and perform BC elimination, storing the
718  // eliminated part of the matrix.
719  if (static_cond)
720  {
722  {
723  static_cond->SetEssentialTrueDofs(ess_tdof_list);
724  static_cond->Finalize(); // finalize Schur complement (to true dofs)
726  static_cond->Finalize(); // finalize eliminated part
727  }
728  A.Reset(&static_cond->GetMatrix(), false);
729  }
730  else
731  {
732  if (!mat_e)
733  {
735  if (P) { ConformingAssemble(); }
736  EliminateVDofs(ess_tdof_list, diag_policy);
737  const int remove_zeros = 0;
738  Finalize(remove_zeros);
739  }
740  if (hybridization)
741  {
742  A.Reset(&hybridization->GetMatrix(), false);
743  }
744  else
745  {
746  A.Reset(mat, false);
747  }
748  }
749 }
750
752  const Vector &b, Vector &x)
753 {
754  if (ext)
755  {
756  ext->RecoverFEMSolution(X, b, x);
757  return;
758  }
759
761  if (!P) // conforming space
762  {
763  if (static_cond)
764  {
765  // Private dofs back solve
766  static_cond->ComputeSolution(b, X, x);
767  }
768  else if (hybridization)
769  {
770  // Primal unknowns recovery
771  hybridization->ComputeSolution(b, X, x);
772  }
773  else
774  {
775  // X and x point to the same data
776
777  // If the validity flags of X's Memory were changed (e.g. if it was
778  // moved to device memory) then we need to tell x about that.
779  x.SyncMemory(X);
780  }
781  }
782  else // non-conforming space
783  {
784  if (static_cond)
785  {
786  // Private dofs back solve
787  static_cond->ComputeSolution(b, X, x);
788  }
789  else if (hybridization)
790  {
791  // Primal unknowns recovery
792  Vector conf_b(P->Width()), conf_x(P->Width());
793  P->MultTranspose(b, conf_b);
795  R->Mult(x, conf_x); // get essential b.c. from x
796  hybridization->ComputeSolution(conf_b, X, conf_x);
797  x.SetSize(P->Height());
798  P->Mult(conf_x, x);
799  }
800  else
801  {
802  // Apply conforming prolongation
803  x.SetSize(P->Height());
804  P->Mult(X, x);
805  }
806  }
807 }
808
810 {
811  if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
812  {
813  return;
814  }
815
816  int num_elements = fes->GetNE();
817  int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
818
819  element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
820  num_elements);
821
822  DenseMatrix tmp;
824
825 #ifdef MFEM_USE_LEGACY_OPENMP
826  #pragma omp parallel for private(tmp,eltrans)
827 #endif
828  for (int i = 0; i < num_elements; i++)
829  {
831  num_dofs_per_el, num_dofs_per_el);
832  const FiniteElement &fe = *fes->GetFE(i);
833 #ifdef MFEM_DEBUG
834  if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
835  mfem_error("BilinearForm::ComputeElementMatrices:"
836  " all elements must have same number of dofs");
837 #endif
838  fes->GetElementTransformation(i, &eltrans);
839
840  dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
841  for (int k = 1; k < dbfi.Size(); k++)
842  {
843  // note: some integrators may not be thread-safe
844  dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
845  elmat += tmp;
846  }
847  elmat.ClearExternalData();
848  }
849 }
850
851 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
852  const Vector &sol, Vector &rhs,
853  DiagonalPolicy dpolicy)
854 {
855  Array<int> ess_dofs, conf_ess_dofs;
856  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
857
858  if (fes->GetVSize() == height)
859  {
860  EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
861  }
862  else
863  {
864  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
865  EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
866  }
867 }
868
869 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
870  DiagonalPolicy dpolicy)
871 {
872  Array<int> ess_dofs, conf_ess_dofs;
873  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
874
875  if (fes->GetVSize() == height)
876  {
877  EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
878  }
879  else
880  {
881  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
882  EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
883  }
884 }
885
887  double value)
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  EliminateEssentialBCFromDofsDiag(ess_dofs, value);
895  }
896  else
897  {
898  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
899  EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
900  }
901 }
902
904  const Vector &sol, Vector &rhs,
905  DiagonalPolicy dpolicy)
906 {
907  for (int i = 0; i < vdofs.Size(); i++)
908  {
909  int vdof = vdofs[i];
910  if ( vdof >= 0 )
911  {
912  mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
913  }
914  else
915  {
916  mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
917  }
918  }
919 }
920
922  DiagonalPolicy dpolicy)
923 {
924  if (mat_e == NULL)
925  {
926  mat_e = new SparseMatrix(height);
927  }
928
929  for (int i = 0; i < vdofs.Size(); i++)
930  {
931  int vdof = vdofs[i];
932  if ( vdof >= 0 )
933  {
934  mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
935  }
936  else
937  {
938  mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
939  }
940  }
941 }
942
944  const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
945  DiagonalPolicy dpolicy)
946 {
947  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
948  MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
949  MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
950
951  for (int i = 0; i < ess_dofs.Size(); i++)
952  if (ess_dofs[i] < 0)
953  {
954  mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
955  }
956 }
957
959  DiagonalPolicy dpolicy)
960 {
961  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
962
963  for (int i = 0; i < ess_dofs.Size(); i++)
964  if (ess_dofs[i] < 0)
965  {
966  mat -> EliminateRowCol (i, dpolicy);
967  }
968 }
969
971  double value)
972 {
973  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
974
975  for (int i = 0; i < ess_dofs.Size(); i++)
976  if (ess_dofs[i] < 0)
977  {
978  mat -> EliminateRowColDiag (i, value);
979  }
980 }
981
983  const Array<int> &vdofs, const Vector &x, Vector &b)
984 {
986  mat->PartMult(vdofs, x, b);
987 }
988
989 void BilinearForm::Mult(const Vector &x, Vector &y) const
990 {
991  if (ext)
992  {
993  ext->Mult(x, y);
994  }
995  else
996  {
997  mat->Mult(x, y);
998  }
999 }
1000
1002 {
1003  bool full_update;
1004
1005  if (nfes && nfes != fes)
1006  {
1007  full_update = true;
1008  fes = nfes;
1009  }
1010  else
1011  {
1012  // Check for different size (e.g. assembled form on non-conforming space)
1013  // or different sequence number.
1014  full_update = (fes->GetVSize() != Height() ||
1015  sequence < fes->GetSequence());
1016  }
1017
1018  delete mat_e;
1019  mat_e = NULL;
1021  delete static_cond;
1022  static_cond = NULL;
1023
1024  if (full_update)
1025  {
1026  delete mat;
1027  mat = NULL;
1028  delete hybridization;
1029  hybridization = NULL;
1030  sequence = fes->GetSequence();
1031  }
1032  else
1033  {
1034  if (mat) { *mat = 0.0; }
1035  if (hybridization) { hybridization->Reset(); }
1036  }
1037
1038  height = width = fes->GetVSize();
1039
1040  if (ext) { ext->Update(); }
1041 }
1042
1044 {
1045  diag_policy = policy;
1046 }
1047
1049 {
1050  delete mat_e;
1051  delete mat;
1052  delete element_matrices;
1053  delete static_cond;
1054  delete hybridization;
1055
1056  if (!extern_bfs)
1057  {
1058  int k;
1059  for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
1060  for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
1061  for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
1062  for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
1063  }
1064
1065  delete ext;
1066 }
1067
1068
1069 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1070  FiniteElementSpace *te_fes)
1071  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1072 {
1073  trial_fes = tr_fes;
1074  test_fes = te_fes;
1075  mat = NULL;
1076  mat_e = NULL;
1077  extern_bfs = 0;
1079  ext = NULL;
1080 }
1081
1082 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1083  FiniteElementSpace *te_fes,
1084  MixedBilinearForm * mbf)
1085  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1086 {
1087  trial_fes = tr_fes;
1088  test_fes = te_fes;
1089  mat = NULL;
1090  mat_e = NULL;
1091  extern_bfs = 1;
1092  ext = NULL;
1093
1094  // Copy the pointers to the integrators
1095  dbfi = mbf->dbfi;
1096  bbfi = mbf->bbfi;
1097  tfbfi = mbf->tfbfi;
1098  btfbfi = mbf->btfbfi;
1099
1100  bbfi_marker = mbf->bbfi_marker;
1102
1104  ext = NULL;
1105 }
1106
1108 {
1109  if (ext)
1110  {
1111  MFEM_ABORT("the assembly level has already been set!");
1112  }
1113  assembly = assembly_level;
1114  switch (assembly)
1115  {
1116  case AssemblyLevel::FULL:
1117  // ext = new FAMixedBilinearFormExtension(this);
1118  // Use the original BilinearForm implementation for now
1119  break;
1121  mfem_error("Element assembly not supported yet... stay tuned!");
1122  // ext = new EAMixedBilinearFormExtension(this);
1123  break;
1125  ext = new PAMixedBilinearFormExtension(this);
1126  break;
1127  case AssemblyLevel::NONE:
1128  mfem_error("Matrix-free action not supported yet... stay tuned!");
1129  // ext = new MFMixedBilinearFormExtension(this);
1130  break;
1131  default:
1132  mfem_error("Unknown assembly level");
1133  }
1134 }
1135
1136 double & MixedBilinearForm::Elem (int i, int j)
1137 {
1138  return (*mat)(i, j);
1139 }
1140
1141 const double & MixedBilinearForm::Elem (int i, int j) const
1142 {
1143  return (*mat)(i, j);
1144 }
1145
1146 void MixedBilinearForm::Mult(const Vector & x, Vector & y) const
1147 {
1148  y = 0.0;
1150 }
1151
1153  const double a) const
1154 {
1155  if (ext)
1156  {
1158  }
1159  else
1160  {
1162  }
1163 }
1164
1166 {
1167  y = 0.0;
1169 }
1170
1172  const double a) const
1173 {
1174  if (ext)
1175  {
1177  }
1178  else
1179  {
1181  }
1182 }
1183
1185 {
1187  {
1188  MFEM_WARNING("MixedBilinearForm::Inverse not possible with this assembly level!");
1189  return NULL;
1190  }
1191  else
1192  {
1193  return mat -> Inverse ();
1194  }
1195 }
1196
1197 void MixedBilinearForm::Finalize (int skip_zeros)
1198 {
1200  {
1201  mat -> Finalize (skip_zeros);
1202  }
1203 }
1204
1206 {
1207  MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
1209  "MixedBilinearForm::GetBlocks: both trial and test spaces "
1210  "must use Ordering::byNODES!");
1211
1212  blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
1213
1214  mat->GetBlocks(blocks);
1215 }
1216
1218 {
1219  dbfi.Append (bfi);
1220 }
1221
1223 {
1224  bbfi.Append (bfi);
1225  bbfi_marker.Append(NULL); // NULL marker means apply everywhere
1226 }
1227
1229  Array<int> &bdr_marker)
1230 {
1231  bbfi.Append (bfi);
1232  bbfi_marker.Append(&bdr_marker);
1233 }
1234
1236 {
1237  tfbfi.Append (bfi);
1238 }
1239
1241 {
1242  btfbfi.Append(bfi);
1243  btfbfi_marker.Append(NULL); // NULL marker means apply everywhere
1244 }
1245
1247  Array<int> &bdr_marker)
1248 {
1249  btfbfi.Append(bfi);
1250  btfbfi_marker.Append(&bdr_marker);
1251 }
1252
1253 void MixedBilinearForm::Assemble (int skip_zeros)
1254 {
1255  if (ext)
1256  {
1257  ext->Assemble();
1258  return;
1259  }
1260
1261  Array<int> tr_vdofs, te_vdofs;
1262  ElementTransformation *eltrans;
1264
1265  Mesh *mesh = test_fes -> GetMesh();
1266
1267  if (mat == NULL)
1268  {
1269  mat = new SparseMatrix(height, width);
1270  }
1271
1272  if (dbfi.Size())
1273  {
1274  for (int i = 0; i < test_fes -> GetNE(); i++)
1275  {
1276  trial_fes -> GetElementVDofs (i, tr_vdofs);
1277  test_fes -> GetElementVDofs (i, te_vdofs);
1278  eltrans = test_fes -> GetElementTransformation (i);
1279  for (int k = 0; k < dbfi.Size(); k++)
1280  {
1281  dbfi[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1282  *test_fes -> GetFE(i),
1283  *eltrans, elemmat);
1284  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1285  }
1286  }
1287  }
1288
1289  if (bbfi.Size())
1290  {
1291  // Which boundary attributes need to be processed?
1292  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1293  mesh->bdr_attributes.Max() : 0);
1294  bdr_attr_marker = 0;
1295  for (int k = 0; k < bbfi.Size(); k++)
1296  {
1297  if (bbfi_marker[k] == NULL)
1298  {
1299  bdr_attr_marker = 1;
1300  break;
1301  }
1302  Array<int> &bdr_marker = *bbfi_marker[k];
1303  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1304  "invalid boundary marker for boundary integrator #"
1305  << k << ", counting from zero");
1306  for (int i = 0; i < bdr_attr_marker.Size(); i++)
1307  {
1308  bdr_attr_marker[i] |= bdr_marker[i];
1309  }
1310  }
1311
1312  for (int i = 0; i < test_fes -> GetNBE(); i++)
1313  {
1314  const int bdr_attr = mesh->GetBdrAttribute(i);
1315  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1316
1317  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1318  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1319  eltrans = test_fes -> GetBdrElementTransformation (i);
1320  for (int k = 0; k < bbfi.Size(); k++)
1321  {
1322  if (bbfi_marker[k] &&
1323  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
1324
1325  bbfi[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1326  *test_fes -> GetBE(i),
1327  *eltrans, elemmat);
1328  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1329  }
1330  }
1331  }
1332
1333  if (tfbfi.Size())
1334  {
1336  Array<int> te_vdofs2;
1337  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1338
1339  int nfaces = mesh->GetNumFaces();
1340  for (int i = 0; i < nfaces; i++)
1341  {
1342  ftr = mesh->GetFaceElementTransformations(i);
1343  trial_fes->GetFaceVDofs(i, tr_vdofs);
1344  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1345  trial_face_fe = trial_fes->GetFaceElement(i);
1346  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1347  if (ftr->Elem2No >= 0)
1348  {
1349  test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
1350  te_vdofs.Append(te_vdofs2);
1351  test_fe2 = test_fes->GetFE(ftr->Elem2No);
1352  }
1353  else
1354  {
1355  // The test_fe2 object is really a dummy and not used on the
1356  // boundaries, but we can't dereference a NULL pointer, and we don't
1357  // want to actually make a fake element.
1358  test_fe2 = test_fe1;
1359  }
1360  for (int k = 0; k < tfbfi.Size(); k++)
1361  {
1362  tfbfi[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1363  *ftr, elemmat);
1365  }
1366  }
1367  }
1368
1369  if (btfbfi.Size())
1370  {
1372  Array<int> te_vdofs2;
1373  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1374
1375  // Which boundary attributes need to be processed?
1376  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1377  mesh->bdr_attributes.Max() : 0);
1378  bdr_attr_marker = 0;
1379  for (int k = 0; k < btfbfi.Size(); k++)
1380  {
1381  if (btfbfi_marker[k] == NULL)
1382  {
1383  bdr_attr_marker = 1;
1384  break;
1385  }
1386  Array<int> &bdr_marker = *btfbfi_marker[k];
1387  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1388  "invalid boundary marker for boundary trace face integrator #"
1389  << k << ", counting from zero");
1390  for (int i = 0; i < bdr_attr_marker.Size(); i++)
1391  {
1392  bdr_attr_marker[i] |= bdr_marker[i];
1393  }
1394  }
1395
1396  for (int i = 0; i < trial_fes -> GetNBE(); i++)
1397  {
1398  const int bdr_attr = mesh->GetBdrAttribute(i);
1399  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1400
1401  ftr = mesh->GetBdrFaceTransformations(i);
1402  if (ftr)
1403  {
1404  trial_fes->GetFaceVDofs(i, tr_vdofs);
1405  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1406  trial_face_fe = trial_fes->GetFaceElement(i);
1407  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1408  // The test_fe2 object is really a dummy and not used on the
1409  // boundaries, but we can't dereference a NULL pointer, and we don't
1410  // want to actually make a fake element.
1411  test_fe2 = test_fe1;
1412  for (int k = 0; k < btfbfi.Size(); k++)
1413  {
1414  if (btfbfi_marker[k] &&
1415  (*btfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
1416
1417  btfbfi[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1418  *ftr, elemmat);
1420  }
1421  }
1422  }
1423  }
1424 }
1425
1427 {
1429  {
1430  MFEM_WARNING("Conforming assemble not supported for this assembly level!");
1431  return;
1432  }
1433
1434  Finalize();
1435
1437  if (P2)
1438  {
1439  SparseMatrix *R = Transpose(*P2);
1440  SparseMatrix *RA = mfem::Mult(*R, *mat);
1441  delete R;
1442  delete mat;
1443  mat = RA;
1444  }
1445
1447  if (P1)
1448  {
1449  SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1450  delete mat;
1451  mat = RAP;
1452  }
1453
1454  height = mat->Height();
1455  width = mat->Width();
1456 }
1457
1458
1460 {
1461  if (dbfi.Size())
1462  {
1463  const FiniteElement &trial_fe = *trial_fes->GetFE(i);
1464  const FiniteElement &test_fe = *test_fes->GetFE(i);
1466  dbfi[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans, elmat);
1467  for (int k = 1; k < dbfi.Size(); k++)
1468  {
1469  dbfi[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans, elemmat);
1470  elmat += elemmat;
1471  }
1472  }
1473  else
1474  {
1477  elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
1478  elmat = 0.0;
1479  }
1480 }
1481
1483 {
1484  if (bbfi.Size())
1485  {
1486  const FiniteElement &trial_be = *trial_fes->GetBE(i);
1487  const FiniteElement &test_be = *test_fes->GetBE(i);
1489  bbfi[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans, elmat);
1490  for (int k = 1; k < bbfi.Size(); k++)
1491  {
1492  bbfi[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans, elemmat);
1493  elmat += elemmat;
1494  }
1495  }
1496  else
1497  {
1500  elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
1501  elmat = 0.0;
1502  }
1503 }
1504
1506  int i, const DenseMatrix &elmat, int skip_zeros)
1507 {
1508  AssembleElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1509 }
1510
1512  int i, const DenseMatrix &elmat, Array<int> &trial_vdofs,
1513  Array<int> &test_vdofs, int skip_zeros)
1514 {
1515  trial_fes->GetElementVDofs(i, trial_vdofs);
1516  test_fes->GetElementVDofs(i, test_vdofs);
1517  if (mat == NULL)
1518  {
1519  mat = new SparseMatrix(height, width);
1520  }
1522 }
1523
1525  int i, const DenseMatrix &elmat, int skip_zeros)
1526 {
1527  AssembleBdrElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1528 }
1529
1531  int i, const DenseMatrix &elmat, Array<int> &trial_vdofs,
1532  Array<int> &test_vdofs, int skip_zeros)
1533 {
1534  trial_fes->GetBdrElementVDofs(i, trial_vdofs);
1535  test_fes->GetBdrElementVDofs(i, test_vdofs);
1536  if (mat == NULL)
1537  {
1538  mat = new SparseMatrix(height, width);
1539  }
1541 }
1542
1544  const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
1545 {
1546  int i, j, k;
1547  Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1548
1549  cols_marker = 0;
1550  for (i = 0; i < trial_fes -> GetNBE(); i++)
1551  if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1552  {
1553  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1554  for (j = 0; j < tr_vdofs.Size(); j++)
1555  {
1556  if ( (k = tr_vdofs[j]) < 0 )
1557  {
1558  k = -1-k;
1559  }
1560  cols_marker[k] = 1;
1561  }
1562  }
1563  mat -> EliminateCols (cols_marker, &sol, &rhs);
1564 }
1565
1567  const Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
1568 {
1569  mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1570 }
1571
1573 {
1574  int i, j, k;
1575  Array<int> te_vdofs;
1576
1577  for (i = 0; i < test_fes -> GetNBE(); i++)
1578  if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1579  {
1580  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1581  for (j = 0; j < te_vdofs.Size(); j++)
1582  {
1583  if ( (k = te_vdofs[j]) < 0 )
1584  {
1585  k = -1-k;
1586  }
1587  mat -> EliminateRow (k);
1588  }
1589  }
1590 }
1591
1593  &trial_tdof_list,
1594  const Array<int> &test_tdof_list,
1595  OperatorHandle &A)
1596
1597 {
1598  if (ext)
1599  {
1600  ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
1601  return;
1602  }
1603
1604  const SparseMatrix *test_P = test_fes->GetConformingProlongation();
1605  const SparseMatrix *trial_P = trial_fes->GetConformingProlongation();
1606
1607  mat->Finalize();
1608
1609  if (test_P) // TODO: Must actually check for trial_P too
1610  {
1611  SparseMatrix *m = RAP(*test_P, *mat, *trial_P);
1612  delete mat;
1613  mat = m;
1614  }
1615
1616  Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
1618  ess_trial_tdof_marker);
1620  ess_test_tdof_marker);
1621
1622  mat_e = new SparseMatrix(mat->Height(), mat->Width());
1623  mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
1624
1625  for (int i=0; i<test_tdof_list.Size(); ++i)
1626  {
1627  mat->EliminateRow(test_tdof_list[i]);
1628  }
1629  mat_e->Finalize();
1630  A.Reset(mat, false);
1631 }
1632
1634  &trial_tdof_list,
1635  const Array<int> &test_tdof_list,
1636  Vector &x, Vector &b,
1637  OperatorHandle &A,
1638  Vector &X, Vector &B)
1639 {
1640  if (ext)
1641  {
1642  ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list, x, b, A, X,
1643  B);
1644  return;
1645  }
1646
1647  const Operator *Pi = this->GetProlongation();
1648  const Operator *Po = this->GetOutputProlongation();
1649  const Operator *Ri = this->GetRestriction();
1650  InitTVectors(Po, Ri, Pi, x, b, X, B);
1651
1652  if (!mat_e)
1653  {
1654  FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list,
1655  A); // Set A = mat_e
1656  }
1657  // Eliminate essential BCs with B -= Ab xb
1659
1660  B.SetSubVector(test_tdof_list, 0.0);
1661 }
1662
1664 {
1665  delete mat;
1666  mat = NULL;
1667  delete mat_e;
1668  mat_e = NULL;
1669  height = test_fes->GetVSize();
1670  width = trial_fes->GetVSize();
1671  if (ext) { ext->Update(); }
1672 }
1673
1675 {
1676  if (mat) { delete mat; }
1677  if (mat_e) { delete mat_e; }
1678  if (!extern_bfs)
1679  {
1680  int i;
1681  for (i = 0; i < dbfi.Size(); i++) { delete dbfi[i]; }
1682  for (i = 0; i < bbfi.Size(); i++) { delete bbfi[i]; }
1683  for (i = 0; i < tfbfi.Size(); i++) { delete tfbfi[i]; }
1684  for (i = 0; i < btfbfi.Size(); i++) { delete btfbfi[i]; }
1685  }
1686  delete ext;
1687 }
1688
1689
1691 {
1692  Array<int> dom_vdofs, ran_vdofs;
1694  const FiniteElement *dom_fe, *ran_fe;
1695  DenseMatrix totelmat, elmat;
1696
1697  if (mat == NULL)
1698  {
1699  mat = new SparseMatrix(height, width);
1700  }
1701
1702  if (dbfi.Size() > 0)
1703  {
1704  for (int i = 0; i < test_fes->GetNE(); i++)
1705  {
1706  trial_fes->GetElementVDofs(i, dom_vdofs);
1707  test_fes->GetElementVDofs(i, ran_vdofs);
1709  dom_fe = trial_fes->GetFE(i);
1710  ran_fe = test_fes->GetFE(i);
1711
1712  dbfi[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1713  for (int j = 1; j < dbfi.Size(); j++)
1714  {
1715  dbfi[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1716  totelmat += elmat;
1717  }
1718  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1719  }
1720  }
1721
1722  if (tfbfi.Size())
1723  {
1724  const int nfaces = test_fes->GetMesh()->GetNumFaces();
1725  for (int i = 0; i < nfaces; i++)
1726  {
1727  trial_fes->GetFaceVDofs(i, dom_vdofs);
1728  test_fes->GetFaceVDofs(i, ran_vdofs);
1730  dom_fe = trial_fes->GetFaceElement(i);
1731  ran_fe = test_fes->GetFaceElement(i);
1732
1733  tfbfi[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1734  for (int j = 1; j < tfbfi.Size(); j++)
1735  {
1736  tfbfi[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1737  totelmat += elmat;
1738  }
1739  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1740  }
1741  }
1742 }
1743
1744 }
Abstract class for Finite Elements.
Definition: fe.hpp:232
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:393
int Size() const
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:498
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
int * GetJ()
Definition: table.hpp:114
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:1007
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:813
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:430
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:1636
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:187
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:422
BilinearFormExtension * ext
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:407
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:172
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:729
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:63
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:146
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
Definition: fespace.cpp:297
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:555
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:157
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:141
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:626
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:69
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. The default is AssemblyLevel::FULL.
Array< BilinearFormIntegrator * > btfbfi
Boundary trace face (skeleton) integrators.
Array< BilinearFormIntegrator * > bbfi
Boundary integrators.
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1908
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:864
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:180
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:1011
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)
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:3943
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:820
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
void SetSize(int m, int n)
Definition: array.hpp:335
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.
Keep the diagonal value.
Definition: matrix.hpp:36
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:40
virtual void AssembleDiagonal(Vector &diag) const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
SparseMatrix & GetMatrix()
Return the serial hybridized matrix.
StaticCondensation * static_cond
Owned.
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:707
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
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(...
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:478
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
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:311
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:450
double * GetData(int k)
Definition: densemat.hpp:813
virtual void Update(FiniteElementSpace *nfes=NULL)
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:384
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
SparseMatrix * mat
Owned.
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:948
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:2224
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Dynamic 2D array using row-major layout.
Definition: array.hpp:316
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the &#39;dofs&#39; array to the given &#39;val&#39;.
Definition: vector.cpp:633
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2300
int SizeI() const
Definition: densemat.hpp:762
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:428
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
SparseMatrix * mat_e
Matrix used to eliminate b.c. Owned.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:441
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:538
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
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:314
MixedBilinearFormExtension * ext
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:548
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 reference to a_{ij}.
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:1266
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Table * GetFaceToElementTable() const
Definition: mesh.cpp:4529
FiniteElementSpace * test_fes
Not owned.
int SizeJ() const
Definition: densemat.hpp:763
double a
Definition: lissajous.cpp:41
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.
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;...
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
Adds new Boundary Integrator. Assumes ownership of bfi.
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:456
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:619
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 associated with i&#39;th element.
Definition: fespace.cpp:1642
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:1354
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:184
FiniteElementSpace * trial_fes
Not owned.
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Definition: fespace.cpp:403
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual ~BilinearForm()
Destroys bilinear form.
Data and methods for partially-assembled mixed bilinear forms.
Vector data type.
Definition: vector.hpp:48
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
Definition: fespace.hpp:524
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 for the i&#39;th boundary element.
Definition: fespace.cpp:1882
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:178
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.
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:315
Array< BilinearFormIntegrator * > dbfi
Domain integrators.
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:657
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:725
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
Operator application: y=A(x).
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void EnableStaticCondensation()
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:137