MFEM  v4.0 Finite element discretization library
bilinearform.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
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 = new 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.areColumnsSorted());
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 {
227  if (!static_cond) { mat->Finalize(skip_zeros); }
228  if (mat_e) { mat_e->Finalize(skip_zeros); }
229  if (static_cond) { static_cond->Finalize(); }
231 }
232
234 {
235  dbfi.Append(bfi);
236 }
237
239 {
240  bbfi.Append (bfi);
241  bbfi_marker.Append(NULL); // NULL marker means apply everywhere
242 }
243
245  Array<int> &bdr_marker)
246 {
247  bbfi.Append (bfi);
248  bbfi_marker.Append(&bdr_marker);
249 }
250
252 {
253  fbfi.Append (bfi);
254 }
255
257 {
258  bfbfi.Append(bfi);
259  bfbfi_marker.Append(NULL); // NULL marker means apply everywhere
260 }
261
263  Array<int> &bdr_marker)
264 {
265  bfbfi.Append(bfi);
266  bfbfi_marker.Append(&bdr_marker);
267 }
268
270 {
271  if (element_matrices)
272  {
274  elmat = element_matrices->GetData(i);
275  return;
276  }
277
278  if (dbfi.Size())
279  {
280  const FiniteElement &fe = *fes->GetFE(i);
282  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
283  for (int k = 1; k < dbfi.Size(); k++)
284  {
285  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
286  elmat += elemmat;
287  }
288  }
289  else
290  {
292  elmat.SetSize(vdofs.Size());
293  elmat = 0.0;
294  }
295 }
296
298  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
299 {
300  fes->GetElementVDofs(i, vdofs);
301  if (static_cond)
302  {
303  static_cond->AssembleMatrix(i, elmat);
304  }
305  else
306  {
307  if (mat == NULL)
308  {
309  AllocMat();
310  }
312  if (hybridization)
313  {
314  hybridization->AssembleMatrix(i, elmat);
315  }
316  }
317 }
318
320  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
321 {
322  fes->GetBdrElementVDofs(i, vdofs);
323  if (static_cond)
324  {
325  static_cond->AssembleBdrMatrix(i, elmat);
326  }
327  else
328  {
329  if (mat == NULL)
330  {
331  AllocMat();
332  }
334  if (hybridization)
335  {
336  hybridization->AssembleBdrMatrix(i, elmat);
337  }
338  }
339 }
340
341 void BilinearForm::Assemble(int skip_zeros)
342 {
343  if (ext)
344  {
345  ext->Assemble();
346  return;
347  }
348
349  ElementTransformation *eltrans;
350  Mesh *mesh = fes -> GetMesh();
351  DenseMatrix elmat, *elmat_p;
352
353  if (mat == NULL)
354  {
355  AllocMat();
356  }
357
358 #ifdef MFEM_USE_LEGACY_OPENMP
359  int free_element_matrices = 0;
360  if (!element_matrices)
361  {
363  free_element_matrices = 1;
364  }
365 #endif
366
367  if (dbfi.Size())
368  {
369  for (int i = 0; i < fes -> GetNE(); i++)
370  {
372  if (element_matrices)
373  {
374  elmat_p = &(*element_matrices)(i);
375  }
376  else
377  {
378  const FiniteElement &fe = *fes->GetFE(i);
379  eltrans = fes->GetElementTransformation(i);
380  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
381  for (int k = 1; k < dbfi.Size(); k++)
382  {
383  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
384  elmat += elemmat;
385  }
386  elmat_p = &elmat;
387  }
388  if (static_cond)
389  {
390  static_cond->AssembleMatrix(i, *elmat_p);
391  }
392  else
393  {
395  if (hybridization)
396  {
397  hybridization->AssembleMatrix(i, *elmat_p);
398  }
399  }
400  }
401  }
402
403  if (bbfi.Size())
404  {
405  // Which boundary attributes need to be processed?
406  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
407  mesh->bdr_attributes.Max() : 0);
408  bdr_attr_marker = 0;
409  for (int k = 0; k < bbfi.Size(); k++)
410  {
411  if (bbfi_marker[k] == NULL)
412  {
413  bdr_attr_marker = 1;
414  break;
415  }
416  Array<int> &bdr_marker = *bbfi_marker[k];
417  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
418  "invalid boundary marker for boundary integrator #"
419  << k << ", counting from zero");
420  for (int i = 0; i < bdr_attr_marker.Size(); i++)
421  {
422  bdr_attr_marker[i] |= bdr_marker[i];
423  }
424  }
425
426  for (int i = 0; i < fes -> GetNBE(); i++)
427  {
428  const int bdr_attr = mesh->GetBdrAttribute(i);
429  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
430
431  const FiniteElement &be = *fes->GetBE(i);
432  fes -> GetBdrElementVDofs (i, vdofs);
433  eltrans = fes -> GetBdrElementTransformation (i);
434  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
435  for (int k = 1; k < bbfi.Size(); k++)
436  {
437  if (bbfi_marker[k] &&
438  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
439
440  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
441  elmat += elemmat;
442  }
443  if (!static_cond)
444  {
446  if (hybridization)
447  {
448  hybridization->AssembleBdrMatrix(i, elmat);
449  }
450  }
451  else
452  {
453  static_cond->AssembleBdrMatrix(i, elmat);
454  }
455  }
456  }
457
458  if (fbfi.Size())
459  {
461  Array<int> vdofs2;
462
463  int nfaces = mesh->GetNumFaces();
464  for (int i = 0; i < nfaces; i++)
465  {
466  tr = mesh -> GetInteriorFaceTransformations (i);
467  if (tr != NULL)
468  {
469  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
470  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
471  vdofs.Append (vdofs2);
472  for (int k = 0; k < fbfi.Size(); k++)
473  {
474  fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
475  *fes -> GetFE (tr -> Elem2No),
476  *tr, elemmat);
477  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
478  }
479  }
480  }
481  }
482
483  if (bfbfi.Size())
484  {
486  const FiniteElement *fe1, *fe2;
487
488  // Which boundary attributes need to be processed?
489  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
490  mesh->bdr_attributes.Max() : 0);
491  bdr_attr_marker = 0;
492  for (int k = 0; k < bfbfi.Size(); k++)
493  {
494  if (bfbfi_marker[k] == NULL)
495  {
496  bdr_attr_marker = 1;
497  break;
498  }
499  Array<int> &bdr_marker = *bfbfi_marker[k];
500  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
501  "invalid boundary marker for boundary face integrator #"
502  << k << ", counting from zero");
503  for (int i = 0; i < bdr_attr_marker.Size(); i++)
504  {
505  bdr_attr_marker[i] |= bdr_marker[i];
506  }
507  }
508
509  for (int i = 0; i < fes -> GetNBE(); i++)
510  {
511  const int bdr_attr = mesh->GetBdrAttribute(i);
512  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
513
514  tr = mesh -> GetBdrFaceTransformations (i);
515  if (tr != NULL)
516  {
517  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
518  fe1 = fes -> GetFE (tr -> Elem1No);
519  // The fe2 object is really a dummy and not used on the boundaries,
520  // but we can't dereference a NULL pointer, and we don't want to
521  // actually make a fake element.
522  fe2 = fe1;
523  for (int k = 0; k < bfbfi.Size(); k++)
524  {
525  if (bfbfi_marker[k] &&
526  (*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
527
528  bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
529  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
530  }
531  }
532  }
533  }
534
535 #ifdef MFEM_USE_LEGACY_OPENMP
536  if (free_element_matrices)
537  {
539  }
540 #endif
541 }
542
544 {
545  // Do not remove zero entries to preserve the symmetric structure of the
546  // matrix which in turn will give rise to symmetric structure in the new
547  // matrix. This ensures that subsequent calls to EliminateRowCol will work
548  // correctly.
549  Finalize(0);
550  MFEM_ASSERT(mat, "the BilinearForm is not assembled");
551
553  if (!P) { return; } // conforming mesh
554
555  SparseMatrix *R = Transpose(*P);
556  SparseMatrix *RA = mfem::Mult(*R, *mat);
557  delete mat;
558  if (mat_e)
559  {
560  SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
561  delete mat_e;
562  mat_e = RAe;
563  }
564  delete R;
565  mat = mfem::Mult(*RA, *P);
566  delete RA;
567  if (mat_e)
568  {
569  SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
570  delete mat_e;
571  mat_e = RAeP;
572  }
573
574  height = mat->Height();
575  width = mat->Width();
576 }
577
578 void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
579  Vector &b, OperatorHandle &A, Vector &X,
580  Vector &B, int copy_interior)
581 {
583
584  if (ext)
585  {
586  ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
587  return;
588  }
589
590  FormSystemMatrix(ess_tdof_list, A);
591
592  // Transform the system and perform the elimination in B, based on the
593  // essential BC values from x. Restrict the BC part of x in X, and set the
594  // non-BC part to zero. Since there is no good initial guess for the Lagrange
595  // multipliers, set X = 0.0 for hybridization.
596  if (static_cond)
597  {
598  // Schur complement reduction to the exposed dofs
599  static_cond->ReduceSystem(x, b, X, B, copy_interior);
600  }
601  else if (!P) // conforming space
602  {
603  if (hybridization)
604  {
605  // Reduction to the Lagrange multipliers system
606  EliminateVDofsInRHS(ess_tdof_list, x, b);
607  hybridization->ReduceRHS(b, B);
608  X.SetSize(B.Size());
609  X = 0.0;
610  }
611  else
612  {
613  // A, X and B point to the same data as mat, x and b
614  EliminateVDofsInRHS(ess_tdof_list, x, b);
615  X.NewMemoryAndSize(x.GetMemory(), x.Size(), false);
616  B.NewMemoryAndSize(b.GetMemory(), b.Size(), false);
617  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
618  }
619  }
620  else // non-conforming space
621  {
622  if (hybridization)
623  {
624  // Reduction to the Lagrange multipliers system
626  Vector conf_b(P->Width()), conf_x(P->Width());
627  P->MultTranspose(b, conf_b);
628  R->Mult(x, conf_x);
629  EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
630  R->MultTranspose(conf_b, b); // store eliminated rhs in b
631  hybridization->ReduceRHS(conf_b, B);
632  X.SetSize(B.Size());
633  X = 0.0;
634  }
635  else
636  {
637  // Variational restriction with P
639  B.SetSize(P->Width());
640  P->MultTranspose(b, B);
641  X.SetSize(R->Height());
642  R->Mult(x, X);
643  EliminateVDofsInRHS(ess_tdof_list, X, B);
644  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
645  }
646  }
647 }
648
649 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
650  OperatorHandle &A)
651 {
652  if (ext)
653  {
654  ext->FormSystemMatrix(ess_tdof_list, A);
655  return;
656  }
657
658  // Finish the matrix assembly and perform BC elimination, storing the
659  // eliminated part of the matrix.
660  if (static_cond)
661  {
663  {
664  static_cond->SetEssentialTrueDofs(ess_tdof_list);
665  static_cond->Finalize(); // finalize Schur complement (to true dofs)
667  static_cond->Finalize(); // finalize eliminated part
668  }
669  A.Reset(&static_cond->GetMatrix(), false);
670  }
671  else
672  {
673  if (!mat_e)
674  {
676  if (P) { ConformingAssemble(); }
677  EliminateVDofs(ess_tdof_list, diag_policy);
678  const int remove_zeros = 0;
679  Finalize(remove_zeros);
680  }
681  if (hybridization)
682  {
683  A.Reset(&hybridization->GetMatrix(), false);
684  }
685  else
686  {
687  A.Reset(mat, false);
688  }
689  }
690 }
691
693  const Vector &b, Vector &x)
694 {
695  if (ext)
696  {
697  ext->RecoverFEMSolution(X, b, x);
698  return;
699  }
700
702  if (!P) // conforming space
703  {
704  if (static_cond)
705  {
706  // Private dofs back solve
707  static_cond->ComputeSolution(b, X, x);
708  }
709  else if (hybridization)
710  {
711  // Primal unknowns recovery
712  hybridization->ComputeSolution(b, X, x);
713  }
714  else
715  {
716  // X and x point to the same data
717
718  // If the validity flags of X's Memory were changed (e.g. if it was
719  // moved to device memory) then we need to tell x about that.
720  x.SyncMemory(X);
721  }
722  }
723  else // non-conforming space
724  {
725  if (static_cond)
726  {
727  // Private dofs back solve
728  static_cond->ComputeSolution(b, X, x);
729  }
730  else if (hybridization)
731  {
732  // Primal unknowns recovery
733  Vector conf_b(P->Width()), conf_x(P->Width());
734  P->MultTranspose(b, conf_b);
736  R->Mult(x, conf_x); // get essential b.c. from x
737  hybridization->ComputeSolution(conf_b, X, conf_x);
738  x.SetSize(P->Height());
739  P->Mult(conf_x, x);
740  }
741  else
742  {
743  // Apply conforming prolongation
744  x.SetSize(P->Height());
745  P->Mult(X, x);
746  }
747  }
748 }
749
751 {
752  if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
753  {
754  return;
755  }
756
757  int num_elements = fes->GetNE();
758  int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
759
760  element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
761  num_elements);
762
763  DenseMatrix tmp;
765
766 #ifdef MFEM_USE_LEGACY_OPENMP
767  #pragma omp parallel for private(tmp,eltrans)
768 #endif
769  for (int i = 0; i < num_elements; i++)
770  {
772  num_dofs_per_el, num_dofs_per_el);
773  const FiniteElement &fe = *fes->GetFE(i);
774 #ifdef MFEM_DEBUG
775  if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
776  mfem_error("BilinearForm::ComputeElementMatrices:"
777  " all elements must have same number of dofs");
778 #endif
779  fes->GetElementTransformation(i, &eltrans);
780
781  dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
782  for (int k = 1; k < dbfi.Size(); k++)
783  {
784  // note: some integrators may not be thread-safe
785  dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
786  elmat += tmp;
787  }
788  elmat.ClearExternalData();
789  }
790 }
791
792 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
793  const Vector &sol, Vector &rhs,
794  DiagonalPolicy dpolicy)
795 {
796  Array<int> ess_dofs, conf_ess_dofs;
797  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
798
799  if (fes->GetVSize() == height)
800  {
801  EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
802  }
803  else
804  {
805  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
806  EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
807  }
808 }
809
810 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
811  DiagonalPolicy dpolicy)
812 {
813  Array<int> ess_dofs, conf_ess_dofs;
814  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
815
816  if (fes->GetVSize() == height)
817  {
818  EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
819  }
820  else
821  {
822  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
823  EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
824  }
825 }
826
828  double value)
829 {
830  Array<int> ess_dofs, conf_ess_dofs;
831  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
832
833  if (fes->GetVSize() == height)
834  {
835  EliminateEssentialBCFromDofsDiag(ess_dofs, value);
836  }
837  else
838  {
839  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
840  EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
841  }
842 }
843
845  const Vector &sol, Vector &rhs,
846  DiagonalPolicy dpolicy)
847 {
848  for (int i = 0; i < vdofs.Size(); i++)
849  {
850  int vdof = vdofs[i];
851  if ( vdof >= 0 )
852  {
853  mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
854  }
855  else
856  {
857  mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
858  }
859  }
860 }
861
863  DiagonalPolicy dpolicy)
864 {
865  if (mat_e == NULL)
866  {
867  mat_e = new SparseMatrix(height);
868  }
869
870  for (int i = 0; i < vdofs.Size(); i++)
871  {
872  int vdof = vdofs[i];
873  if ( vdof >= 0 )
874  {
875  mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
876  }
877  else
878  {
879  mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
880  }
881  }
882 }
883
885  const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
886  DiagonalPolicy dpolicy)
887 {
888  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
889  MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
890  MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
891
892  for (int i = 0; i < ess_dofs.Size(); i++)
893  if (ess_dofs[i] < 0)
894  {
895  mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
896  }
897 }
898
900  DiagonalPolicy dpolicy)
901 {
902  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
903
904  for (int i = 0; i < ess_dofs.Size(); i++)
905  if (ess_dofs[i] < 0)
906  {
907  mat -> EliminateRowCol (i, dpolicy);
908  }
909 }
910
912  double value)
913 {
914  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
915
916  for (int i = 0; i < ess_dofs.Size(); i++)
917  if (ess_dofs[i] < 0)
918  {
919  mat -> EliminateRowColDiag (i, value);
920  }
921 }
922
924  const Array<int> &vdofs, const Vector &x, Vector &b)
925 {
927  mat->PartMult(vdofs, x, b);
928 }
929
931 {
932  bool full_update;
933
934  if (nfes && nfes != fes)
935  {
936  full_update = true;
937  fes = nfes;
938  }
939  else
940  {
941  // Check for different size (e.g. assembled form on non-conforming space)
942  // or different sequence number.
943  full_update = (fes->GetVSize() != Height() ||
944  sequence < fes->GetSequence());
945  }
946
947  delete mat_e;
948  mat_e = NULL;
950  delete static_cond;
951  static_cond = NULL;
952
953  if (full_update)
954  {
955  delete mat;
956  mat = NULL;
957  delete hybridization;
958  hybridization = NULL;
959  sequence = fes->GetSequence();
960  }
961  else
962  {
963  if (mat) { *mat = 0.0; }
964  if (hybridization) { hybridization->Reset(); }
965  }
966
967  height = width = fes->GetVSize();
968
969  if (ext) { ext->Update(); }
970 }
971
973 {
974  diag_policy = policy;
975 }
976
978 {
979  delete mat_e;
980  delete mat;
981  delete element_matrices;
982  delete static_cond;
983  delete hybridization;
984
985  if (!extern_bfs)
986  {
987  int k;
988  for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
989  for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
990  for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
991  for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
992  }
993
994  delete ext;
995 }
996
997
998 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
999  FiniteElementSpace *te_fes)
1000  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1001 {
1002  trial_fes = tr_fes;
1003  test_fes = te_fes;
1004  mat = NULL;
1005  extern_bfs = 0;
1006 }
1007
1008 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1009  FiniteElementSpace *te_fes,
1010  MixedBilinearForm * mbf)
1011  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1012 {
1013  trial_fes = tr_fes;
1014  test_fes = te_fes;
1015  mat = NULL;
1016  extern_bfs = 1;
1017
1018  // Copy the pointers to the integrators
1019  dom = mbf->dom;
1020  bdr = mbf->bdr;
1021  skt = mbf->skt;
1022 }
1023
1024 double & MixedBilinearForm::Elem (int i, int j)
1025 {
1026  return (*mat)(i, j);
1027 }
1028
1029 const double & MixedBilinearForm::Elem (int i, int j) const
1030 {
1031  return (*mat)(i, j);
1032 }
1033
1034 void MixedBilinearForm::Mult (const Vector & x, Vector & y) const
1035 {
1036  mat -> Mult (x, y);
1037 }
1038
1040  const double a) const
1041 {
1042  mat -> AddMult (x, y, a);
1043 }
1044
1046  const double a) const
1047 {
1048  mat -> AddMultTranspose (x, y, a);
1049 }
1050
1052 {
1053  return mat -> Inverse ();
1054 }
1055
1056 void MixedBilinearForm::Finalize (int skip_zeros)
1057 {
1058  mat -> Finalize (skip_zeros);
1059 }
1060
1062 {
1063  MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
1065  "MixedBilinearForm::GetBlocks: both trial and test spaces "
1066  "must use Ordering::byNODES!");
1067
1068  blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
1069
1070  mat->GetBlocks(blocks);
1071 }
1072
1074 {
1075  dom.Append (bfi);
1076 }
1077
1079 {
1080  bdr.Append (bfi);
1081 }
1082
1084 {
1085  skt.Append (bfi);
1086 }
1087
1088 void MixedBilinearForm::Assemble (int skip_zeros)
1089 {
1090  int i, k;
1091  Array<int> tr_vdofs, te_vdofs;
1092  ElementTransformation *eltrans;
1093  DenseMatrix elemmat;
1094
1095  Mesh *mesh = test_fes -> GetMesh();
1096
1097  if (mat == NULL)
1098  {
1099  mat = new SparseMatrix(height, width);
1100  }
1101
1102  if (dom.Size())
1103  {
1104  for (i = 0; i < test_fes -> GetNE(); i++)
1105  {
1106  trial_fes -> GetElementVDofs (i, tr_vdofs);
1107  test_fes -> GetElementVDofs (i, te_vdofs);
1108  eltrans = test_fes -> GetElementTransformation (i);
1109  for (k = 0; k < dom.Size(); k++)
1110  {
1111  dom[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1112  *test_fes -> GetFE(i),
1113  *eltrans, elemmat);
1114  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1115  }
1116  }
1117  }
1118
1119  if (bdr.Size())
1120  {
1121  for (i = 0; i < test_fes -> GetNBE(); i++)
1122  {
1123  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1124  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1125  eltrans = test_fes -> GetBdrElementTransformation (i);
1126  for (k = 0; k < bdr.Size(); k++)
1127  {
1128  bdr[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1129  *test_fes -> GetBE(i),
1130  *eltrans, elemmat);
1131  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1132  }
1133  }
1134  }
1135
1136  if (skt.Size())
1137  {
1139  Array<int> te_vdofs2;
1140  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1141
1142  int nfaces = mesh->GetNumFaces();
1143  for (i = 0; i < nfaces; i++)
1144  {
1145  ftr = mesh->GetFaceElementTransformations(i);
1146  trial_fes->GetFaceVDofs(i, tr_vdofs);
1147  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1148  trial_face_fe = trial_fes->GetFaceElement(i);
1149  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1150  if (ftr->Elem2No >= 0)
1151  {
1152  test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
1153  te_vdofs.Append(te_vdofs2);
1154  test_fe2 = test_fes->GetFE(ftr->Elem2No);
1155  }
1156  else
1157  {
1158  // The test_fe2 object is really a dummy and not used on the
1159  // boundaries, but we can't dereference a NULL pointer, and we don't
1160  // want to actually make a fake element.
1161  test_fe2 = test_fe1;
1162  }
1163  for (int k = 0; k < skt.Size(); k++)
1164  {
1165  skt[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1166  *ftr, elemmat);
1168  }
1169  }
1170  }
1171 }
1172
1174 {
1175  Finalize();
1176
1178  if (P2)
1179  {
1180  SparseMatrix *R = Transpose(*P2);
1181  SparseMatrix *RA = mfem::Mult(*R, *mat);
1182  delete R;
1183  delete mat;
1184  mat = RA;
1185  }
1186
1188  if (P1)
1189  {
1190  SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1191  delete mat;
1192  mat = RAP;
1193  }
1194
1195  height = mat->Height();
1196  width = mat->Width();
1197 }
1198
1200  Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
1201 {
1202  int i, j, k;
1203  Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1204
1205  cols_marker = 0;
1206  for (i = 0; i < trial_fes -> GetNBE(); i++)
1207  if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1208  {
1209  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1210  for (j = 0; j < tr_vdofs.Size(); j++)
1211  {
1212  if ( (k = tr_vdofs[j]) < 0 )
1213  {
1214  k = -1-k;
1215  }
1216  cols_marker[k] = 1;
1217  }
1218  }
1219  mat -> EliminateCols (cols_marker, &sol, &rhs);
1220 }
1221
1223  Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
1224 {
1225  mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1226 }
1227
1229 {
1230  int i, j, k;
1231  Array<int> te_vdofs;
1232
1233  for (i = 0; i < test_fes -> GetNBE(); i++)
1234  if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1235  {
1236  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1237  for (j = 0; j < te_vdofs.Size(); j++)
1238  {
1239  if ( (k = te_vdofs[j]) < 0 )
1240  {
1241  k = -1-k;
1242  }
1243  mat -> EliminateRow (k);
1244  }
1245  }
1246 }
1247
1249 {
1250  delete mat;
1251  mat = NULL;
1252  height = test_fes->GetVSize();
1253  width = trial_fes->GetVSize();
1254 }
1255
1257 {
1258  if (mat) { delete mat; }
1259  if (!extern_bfs)
1260  {
1261  int i;
1262  for (i = 0; i < dom.Size(); i++) { delete dom[i]; }
1263  for (i = 0; i < bdr.Size(); i++) { delete bdr[i]; }
1264  for (i = 0; i < skt.Size(); i++) { delete skt[i]; }
1265  }
1266 }
1267
1268
1270 {
1271  Array<int> dom_vdofs, ran_vdofs;
1273  const FiniteElement *dom_fe, *ran_fe;
1274  DenseMatrix totelmat, elmat;
1275
1276  if (mat == NULL)
1277  {
1278  mat = new SparseMatrix(height, width);
1279  }
1280
1281  if (dom.Size() > 0)
1282  {
1283  for (int i = 0; i < test_fes->GetNE(); i++)
1284  {
1285  trial_fes->GetElementVDofs(i, dom_vdofs);
1286  test_fes->GetElementVDofs(i, ran_vdofs);
1288  dom_fe = trial_fes->GetFE(i);
1289  ran_fe = test_fes->GetFE(i);
1290
1291  dom[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1292  for (int j = 1; j < dom.Size(); j++)
1293  {
1294  dom[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1295  totelmat += elmat;
1296  }
1297  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1298  }
1299  }
1300
1301  if (skt.Size())
1302  {
1303  const int nfaces = test_fes->GetMesh()->GetNumFaces();
1304  for (int i = 0; i < nfaces; i++)
1305  {
1306  trial_fes->GetFaceVDofs(i, dom_vdofs);
1307  test_fes->GetFaceVDofs(i, ran_vdofs);
1309  dom_fe = trial_fes->GetFaceElement(i);
1310  ran_fe = test_fes->GetFaceElement(i);
1311
1312  skt[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1313  for (int j = 1; j < skt.Size(); j++)
1314  {
1315  skt[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1316  totelmat += elmat;
1317  }
1318  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1319  }
1320  }
1321 }
1322
1323 }
Abstract class for Finite Elements.
Definition: fe.hpp:229
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:359
int Size() const
Logical size of the array.
Definition: array.hpp:118
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:347
int * GetJ()
Definition: table.hpp:114
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:976
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Array< BilinearFormIntegrator * > skt
Trace face (skeleton) integrators.
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:769
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse.
void SetConstraintIntegrator(BilinearFormIntegrator *c_integ)
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
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1563
void SyncMemory(const Vector &v)
Update the memory location of the vector to match v.
Definition: vector.hpp:180
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:367
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:400
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:726
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:42
void EliminateEssentialBCFromTrialDofs(Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
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().
void EliminateTrialDofs(Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
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:552
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:150
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:417
int * GetI()
Return the array I.
Definition: sparsemat.hpp:141
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:66
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:154
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1807
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:477
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:831
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:173
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
Data and methods for partially-assembled bilinear forms.
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Definition: sparsemat.cpp:1008
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 EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
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:3663
Array< BilinearFormIntegrator * > bdr
Boundary integrators.
Adds a domain integrator. Assumes ownership of bfi.
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:776
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:374
void SetSize(int m, int n)
Definition: array.hpp:323
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
Data type sparse matrix.
Definition: sparsemat.hpp:40
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
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:690
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:272
SparseMatrix * mat
Sparse matrix to be associated with the form. Owned.
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:146
bool areColumnsSorted() const
Definition: sparsemat.hpp:373
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
Returns the transformation defining the given face element in a user-defined variable.
Definition: mesh.cpp:473
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 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)
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
double * GetData(int k)
Definition: densemat.hpp:728
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.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:336
SparseMatrix * mat
Owned.
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:2187
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Dynamic 2D array using row-major layout.
Definition: array.hpp:304
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:2263
int SizeI() const
Definition: densemat.hpp:693
void ComputeElementMatrix(int i, DenseMatrix &elmat)
bool Finalized() const
Definition: sparsemat.hpp:372
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:398
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:23
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:179
void AssembleElementMatrix(int i, const DenseMatrix &elmat, Array< int > &vdofs, int skip_zeros=1)
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:85
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:311
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:545
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 ClearExternalData()
Definition: densemat.hpp:78
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Array< BilinearFormIntegrator * > dom
Domain integrators.
DiagonalPolicy diag_policy
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:227
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:59
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
Table * GetFaceToElementTable() const
Definition: mesh.cpp:4224
FiniteElementSpace * test_fes
Not owned.
int SizeJ() const
Definition: densemat.hpp:694
virtual void EliminateTestDofs(Array< int > &bdr_attr_is_ess)
int extern_bfs
Indicates the BilinearFormIntegrators stored in dom, bdr, and skt are owned by another MixedBilinearF...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AssembleMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:189
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:449
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:616
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:1541
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, Array< int > &vdofs, int skip_zeros=1)
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.
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual ~BilinearForm()
Destroys bilinear form.
Vector data type.
Definition: vector.hpp:48
Adds a boundary integrator. Assumes ownership of bfi.
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:481
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:1781
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:88
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:292
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:614
AssemblyLevel assembly
The form assembly level (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:656
void UseSparsity(int *I, int *J, bool isSorted)
Use the given CSR sparsity pattern to allocate the internal SparseMatrix.
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:25
void EnableStaticCondensation()
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:137