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