MFEM  v3.3.2
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
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of class BilinearForm
13 
14 #include "fem.hpp"
15 #include <cmath>
16 
17 namespace mfem
18 {
19 
21 {
22  if (static_cond) { return; }
23 
24  if (precompute_sparsity == 0 || fes->GetVDim() > 1)
25  {
26  mat = new SparseMatrix(height);
27  return;
28  }
29 
30  const Table &elem_dof = fes->GetElementToDofTable();
31  Table dof_dof;
32 
33  if (fbfi.Size() > 0)
34  {
35  // the sparsity pattern is defined from the map: face->element->dof
36  Table face_dof, dof_face;
37  {
38  Table *face_elem = fes->GetMesh()->GetFaceToElementTable();
39  mfem::Mult(*face_elem, elem_dof, face_dof);
40  delete face_elem;
41  }
42  Transpose(face_dof, dof_face, height);
43  mfem::Mult(dof_face, face_dof, dof_dof);
44  }
45  else
46  {
47  // the sparsity pattern is defined from the map: element->dof
48  Table dof_elem;
49  Transpose(elem_dof, dof_elem, height);
50  mfem::Mult(dof_elem, elem_dof, dof_dof);
51  }
52 
53  dof_dof.SortRows();
54 
55  int *I = dof_dof.GetI();
56  int *J = dof_dof.GetJ();
57  double *data = new double[I[height]];
58 
59  mat = new SparseMatrix(I, J, data, height, height, true, true, true);
60  *mat = 0.0;
61 
62  dof_dof.LoseData();
63 }
64 
66  : Matrix (f->GetVSize())
67 {
68  fes = f;
69  sequence = f->GetSequence();
70  mat = mat_e = NULL;
71  extern_bfs = 0;
72  element_matrices = NULL;
73  static_cond = NULL;
74  hybridization = NULL;
76 }
77 
79  : Matrix (f->GetVSize())
80 {
81  int i;
83 
84  fes = f;
85  sequence = f->GetSequence();
86  mat_e = NULL;
87  extern_bfs = 1;
88  element_matrices = NULL;
89  static_cond = NULL;
90  hybridization = NULL;
92 
93  bfi = bf->GetDBFI();
94  dbfi.SetSize (bfi->Size());
95  for (i = 0; i < bfi->Size(); i++)
96  {
97  dbfi[i] = (*bfi)[i];
98  }
99 
100  bfi = bf->GetBBFI();
101  bbfi.SetSize (bfi->Size());
102  for (i = 0; i < bfi->Size(); i++)
103  {
104  bbfi[i] = (*bfi)[i];
105  }
106 
107  bfi = bf->GetFBFI();
108  fbfi.SetSize (bfi->Size());
109  for (i = 0; i < bfi->Size(); i++)
110  {
111  fbfi[i] = (*bfi)[i];
112  }
113 
114  bfi = bf->GetBFBFI();
115  bfbfi.SetSize (bfi->Size());
116  for (i = 0; i < bfi->Size(); i++)
117  {
118  bfbfi[i] = (*bfi)[i];
119  }
120 
121  AllocMat();
122 }
123 
125 {
126  delete static_cond;
129  {
130  bool symmetric = false; // TODO
131  bool block_diagonal = false; // TODO
132  static_cond->Init(symmetric, block_diagonal);
133  }
134  else
135  {
136  delete static_cond;
137  static_cond = NULL;
138  }
139 }
140 
142  BilinearFormIntegrator *constr_integ,
143  const Array<int> &ess_tdof_list)
144 {
145  delete hybridization;
146  hybridization = new Hybridization(fes, constr_space);
147  hybridization->SetConstraintIntegrator(constr_integ);
148  hybridization->Init(ess_tdof_list);
149 }
150 
151 void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
152 {
153  if (static_cond) { return; }
154 
155  if (mat)
156  {
157  if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
158  {
159  return; // mat is already using the given sparsity
160  }
161  delete mat;
162  }
163  height = width = fes->GetVSize();
164  mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
165 }
166 
168 {
169  MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
170  "invalid matrix A dimensions: "
171  << A.Height() << " x " << A.Width());
172  MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
173 
174  UseSparsity(A.GetI(), A.GetJ(), A.areColumnsSorted());
175 }
176 
177 double& BilinearForm::Elem (int i, int j)
178 {
179  return mat -> Elem(i,j);
180 }
181 
182 const double& BilinearForm::Elem (int i, int j) const
183 {
184  return mat -> Elem(i,j);
185 }
186 
188 {
189  return mat -> Inverse();
190 }
191 
192 void BilinearForm::Finalize (int skip_zeros)
193 {
194  if (!static_cond) { mat->Finalize(skip_zeros); }
195  if (mat_e) { mat_e->Finalize(skip_zeros); }
196  if (static_cond) { static_cond->Finalize(); }
198 }
199 
201 {
202  dbfi.Append (bfi);
203 }
204 
206 {
207  bbfi.Append (bfi);
208 }
209 
211 {
212  fbfi.Append (bfi);
213 }
214 
216 {
217  bfbfi.Append(bfi);
218  bfbfi_marker.Append(NULL); // NULL marker means apply everywhere
219 }
220 
222  Array<int> &bdr_marker)
223 {
224  bfbfi.Append(bfi);
225  bfbfi_marker.Append(&bdr_marker);
226 }
227 
229 {
230  if (element_matrices)
231  {
233  elmat = element_matrices->GetData(i);
234  return;
235  }
236 
237  if (dbfi.Size())
238  {
239  const FiniteElement &fe = *fes->GetFE(i);
241  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
242  for (int k = 1; k < dbfi.Size(); k++)
243  {
244  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
245  elmat += elemmat;
246  }
247  }
248  else
249  {
251  elmat.SetSize(vdofs.Size());
252  elmat = 0.0;
253  }
254 }
255 
257  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
258 {
259  fes->GetElementVDofs(i, vdofs);
260  if (static_cond)
261  {
262  static_cond->AssembleMatrix(i, elmat);
263  }
264  else
265  {
266  if (mat == NULL)
267  {
268  AllocMat();
269  }
270  mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
271  if (hybridization)
272  {
273  hybridization->AssembleMatrix(i, elmat);
274  }
275  }
276 }
277 
279  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
280 {
281  fes->GetBdrElementVDofs(i, vdofs);
282  if (static_cond)
283  {
284  static_cond->AssembleBdrMatrix(i, elmat);
285  }
286  else
287  {
288  if (mat == NULL)
289  {
290  AllocMat();
291  }
292  mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
293  if (hybridization)
294  {
295  hybridization->AssembleBdrMatrix(i, elmat);
296  }
297  }
298 }
299 
300 void BilinearForm::Assemble (int skip_zeros)
301 {
302  ElementTransformation *eltrans;
303  Mesh *mesh = fes -> GetMesh();
304  DenseMatrix elmat, *elmat_p;
305 
306  int i;
307 
308  if (mat == NULL)
309  {
310  AllocMat();
311  }
312 
313 #ifdef MFEM_USE_OPENMP
314  int free_element_matrices = 0;
315  if (!element_matrices)
316  {
318  free_element_matrices = 1;
319  }
320 #endif
321 
322  if (dbfi.Size())
323  {
324  for (i = 0; i < fes -> GetNE(); i++)
325  {
327  if (element_matrices)
328  {
329  elmat_p = &(*element_matrices)(i);
330  }
331  else
332  {
333  const FiniteElement &fe = *fes->GetFE(i);
334  eltrans = fes->GetElementTransformation(i);
335  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
336  for (int k = 1; k < dbfi.Size(); k++)
337  {
338  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
339  elmat += elemmat;
340  }
341  elmat_p = &elmat;
342  }
343  if (static_cond)
344  {
345  static_cond->AssembleMatrix(i, *elmat_p);
346  }
347  else
348  {
349  mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
350  if (hybridization)
351  {
352  hybridization->AssembleMatrix(i, *elmat_p);
353  }
354  }
355  }
356  }
357 
358  if (bbfi.Size())
359  {
360  for (i = 0; i < fes -> GetNBE(); i++)
361  {
362  const FiniteElement &be = *fes->GetBE(i);
363  fes -> GetBdrElementVDofs (i, vdofs);
364  eltrans = fes -> GetBdrElementTransformation (i);
365  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
366  for (int k = 1; k < bbfi.Size(); k++)
367  {
368  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
369  elmat += elemmat;
370  }
371  if (!static_cond)
372  {
373  mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
374  if (hybridization)
375  {
376  hybridization->AssembleBdrMatrix(i, elmat);
377  }
378  }
379  else
380  {
381  static_cond->AssembleBdrMatrix(i, elmat);
382  }
383  }
384  }
385 
386  if (fbfi.Size())
387  {
389  Array<int> vdofs2;
390 
391  int nfaces = mesh->GetNumFaces();
392  for (i = 0; i < nfaces; i++)
393  {
394  tr = mesh -> GetInteriorFaceTransformations (i);
395  if (tr != NULL)
396  {
397  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
398  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
399  vdofs.Append (vdofs2);
400  for (int k = 0; k < fbfi.Size(); k++)
401  {
402  fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
403  *fes -> GetFE (tr -> Elem2No),
404  *tr, elemmat);
405  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
406  }
407  }
408  }
409  }
410 
411  if (bfbfi.Size())
412  {
414  const FiniteElement *fe1, *fe2;
415 
416  // Which boundary attributes need to be processed?
417  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
418  mesh->bdr_attributes.Max() : 0);
419  bdr_attr_marker = 0;
420  for (int k = 0; k < bfbfi.Size(); k++)
421  {
422  if (bfbfi_marker[k] == NULL)
423  {
424  bdr_attr_marker = 1;
425  break;
426  }
427  Array<int> &bdr_marker = *bfbfi_marker[k];
428  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
429  "invalid boundary marker for boundary face integrator #"
430  << k << ", counting from zero");
431  for (int i = 0; i < bdr_attr_marker.Size(); i++)
432  {
433  bdr_attr_marker[i] |= bdr_marker[i];
434  }
435  }
436 
437  for (i = 0; i < fes -> GetNBE(); i++)
438  {
439  const int bdr_attr = mesh->GetBdrAttribute(i);
440  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
441 
442  tr = mesh -> GetBdrFaceTransformations (i);
443  if (tr != NULL)
444  {
445  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
446  fe1 = fes -> GetFE (tr -> Elem1No);
447  // The fe2 object is really a dummy and not used on the boundaries,
448  // but we can't dereference a NULL pointer, and we don't want to
449  // actually make a fake element.
450  fe2 = fe1;
451  for (int k = 0; k < bfbfi.Size(); k++)
452  {
453  if (bfbfi_marker[k] &&
454  (*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
455 
456  bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
457  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
458  }
459  }
460  }
461  }
462 
463 #ifdef MFEM_USE_OPENMP
464  if (free_element_matrices)
465  {
467  }
468 #endif
469 }
470 
472 {
473  // Do not remove zero entries to preserve the symmetric structure of the
474  // matrix which in turn will give rise to symmetric structure in the new
475  // matrix. This ensures that subsequent calls to EliminateRowCol will work
476  // correctly.
477  Finalize(0);
478  MFEM_ASSERT(mat, "the BilinearForm is not assembled");
479 
481  if (!P) { return; } // conforming mesh
482 
483  SparseMatrix *R = Transpose(*P);
484  SparseMatrix *RA = mfem::Mult(*R, *mat);
485  delete mat;
486  if (mat_e)
487  {
488  SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
489  delete mat_e;
490  mat_e = RAe;
491  }
492  delete R;
493  mat = mfem::Mult(*RA, *P);
494  delete RA;
495  if (mat_e)
496  {
497  SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
498  delete mat_e;
499  mat_e = RAeP;
500  }
501 
502  height = mat->Height();
503  width = mat->Width();
504 }
505 
506 void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
507  Vector &x, Vector &b,
508  SparseMatrix &A, Vector &X, Vector &B,
509  int copy_interior)
510 {
512 
513  FormSystemMatrix(ess_tdof_list, A);
514 
515  // Transform the system and perform the elimination in B, based on the
516  // essential BC values from x. Restrict the BC part of x in X, and set the
517  // non-BC part to zero. Since there is no good initial guess for the Lagrange
518  // multipliers, set X = 0.0 for hybridization.
519  if (static_cond)
520  {
521  // Schur complement reduction to the exposed dofs
522  static_cond->ReduceSystem(x, b, X, B, copy_interior);
523  }
524  else if (!P) // conforming space
525  {
526  if (hybridization)
527  {
528  // Reduction to the Lagrange multipliers system
529  EliminateVDofsInRHS(ess_tdof_list, x, b);
530  hybridization->ReduceRHS(b, B);
531  X.SetSize(B.Size());
532  X = 0.0;
533  }
534  else
535  {
536  // A, X and B point to the same data as mat, x and b
537  EliminateVDofsInRHS(ess_tdof_list, x, b);
538  X.NewDataAndSize(x.GetData(), x.Size());
539  B.NewDataAndSize(b.GetData(), b.Size());
540  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
541  }
542  }
543  else // non-conforming space
544  {
545  if (hybridization)
546  {
547  // Reduction to the Lagrange multipliers system
549  Vector conf_b(P->Width()), conf_x(P->Width());
550  P->MultTranspose(b, conf_b);
551  R->Mult(x, conf_x);
552  EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
553  R->MultTranspose(conf_b, b); // store eliminated rhs in b
554  hybridization->ReduceRHS(conf_b, B);
555  X.SetSize(B.Size());
556  X = 0.0;
557  }
558  else
559  {
560  // Variational restriction with P
562  B.SetSize(P->Width());
563  P->MultTranspose(b, B);
564  X.SetSize(R->Height());
565  R->Mult(x, X);
566  EliminateVDofsInRHS(ess_tdof_list, X, B);
567  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
568  }
569  }
570 }
571 
572 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
573  SparseMatrix &A)
574 {
575  // Finish the matrix assembly and perform BC elimination, storing the
576  // eliminated part of the matrix.
577  const int keep_diag = 1;
578  if (static_cond)
579  {
581  {
582  static_cond->SetEssentialTrueDofs(ess_tdof_list);
583  static_cond->Finalize(); // finalize Schur complement (to true dofs)
585  static_cond->Finalize(); // finalize eliminated part
586  }
588  }
589  else
590  {
591  if (!mat_e)
592  {
594  if (P) { ConformingAssemble(); }
595  EliminateVDofs(ess_tdof_list, keep_diag);
596  const int remove_zeros = 0;
597  Finalize(remove_zeros);
598  }
599  if (hybridization)
600  {
602  }
603  else
604  {
605  A.MakeRef(*mat);
606  }
607  }
608 }
609 
611  const Vector &b, Vector &x)
612 {
614  if (!P) // conforming space
615  {
616  if (static_cond)
617  {
618  // Private dofs back solve
619  static_cond->ComputeSolution(b, X, x);
620  }
621  else if (hybridization)
622  {
623  // Primal unknowns recovery
624  hybridization->ComputeSolution(b, X, x);
625  }
626  else
627  {
628  // X and x point to the same data
629  }
630  }
631  else // non-conforming space
632  {
633  if (static_cond)
634  {
635  // Private dofs back solve
636  static_cond->ComputeSolution(b, X, x);
637  }
638  else if (hybridization)
639  {
640  // Primal unknowns recovery
641  Vector conf_b(P->Width()), conf_x(P->Width());
642  P->MultTranspose(b, conf_b);
644  R->Mult(x, conf_x); // get essential b.c. from x
645  hybridization->ComputeSolution(conf_b, X, conf_x);
646  x.SetSize(P->Height());
647  P->Mult(conf_x, x);
648  }
649  else
650  {
651  // Apply conforming prolongation
652  x.SetSize(P->Height());
653  P->Mult(X, x);
654  }
655  }
656 }
657 
659 {
660  if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
661  {
662  return;
663  }
664 
665  int num_elements = fes->GetNE();
666  int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
667 
668  element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
669  num_elements);
670 
671  DenseMatrix tmp;
673 
674 #ifdef MFEM_USE_OPENMP
675  #pragma omp parallel for private(tmp,eltrans)
676 #endif
677  for (int i = 0; i < num_elements; i++)
678  {
680  num_dofs_per_el, num_dofs_per_el);
681  const FiniteElement &fe = *fes->GetFE(i);
682 #ifdef MFEM_DEBUG
683  if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
684  mfem_error("BilinearForm::ComputeElementMatrices:"
685  " all elements must have same number of dofs");
686 #endif
687  fes->GetElementTransformation(i, &eltrans);
688 
689  dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
690  for (int k = 1; k < dbfi.Size(); k++)
691  {
692  // note: some integrators may not be thread-safe
693  dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
694  elmat += tmp;
695  }
696  elmat.ClearExternalData();
697  }
698 }
699 
700 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
701  Vector &sol, Vector &rhs, int d)
702 {
703  Array<int> ess_dofs, conf_ess_dofs;
704  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
705 
706  if (fes->GetVSize() == height)
707  {
708  EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, d);
709  }
710  else
711  {
712  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
713  EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, d);
714  }
715 }
716 
717 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
718  int d)
719 {
720  Array<int> ess_dofs, conf_ess_dofs;
721  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
722 
723  if (fes->GetVSize() == height)
724  {
725  EliminateEssentialBCFromDofs(ess_dofs, d);
726  }
727  else
728  {
729  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
730  EliminateEssentialBCFromDofs(conf_ess_dofs, d);
731  }
732 }
733 
735  double value)
736 {
737  Array<int> ess_dofs, conf_ess_dofs;
738  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
739 
740  if (fes->GetVSize() == height)
741  {
742  EliminateEssentialBCFromDofsDiag(ess_dofs, value);
743  }
744  else
745  {
746  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
747  EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
748  }
749 }
750 
752  Vector &sol, Vector &rhs, int d)
753 {
754  for (int i = 0; i < vdofs.Size(); i++)
755  {
756  int vdof = vdofs[i];
757  if ( vdof >= 0 )
758  {
759  mat -> EliminateRowCol (vdof, sol(vdof), rhs, d);
760  }
761  else
762  {
763  mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, d);
764  }
765  }
766 }
767 
768 void BilinearForm::EliminateVDofs(const Array<int> &vdofs, int d)
769 {
770  if (mat_e == NULL)
771  {
772  mat_e = new SparseMatrix(height);
773  }
774 
775  for (int i = 0; i < vdofs.Size(); i++)
776  {
777  int vdof = vdofs[i];
778  if ( vdof >= 0 )
779  {
780  mat -> EliminateRowCol (vdof, *mat_e, d);
781  }
782  else
783  {
784  mat -> EliminateRowCol (-1-vdof, *mat_e, d);
785  }
786  }
787 }
788 
790  const Array<int> &ess_dofs, Vector &sol, Vector &rhs, int d)
791 {
792  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
793  MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
794  MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
795 
796  for (int i = 0; i < ess_dofs.Size(); i++)
797  if (ess_dofs[i] < 0)
798  {
799  mat -> EliminateRowCol (i, sol(i), rhs, d);
800  }
801 }
802 
804  int d)
805 {
806  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
807 
808  for (int i = 0; i < ess_dofs.Size(); i++)
809  if (ess_dofs[i] < 0)
810  {
811  mat -> EliminateRowCol (i, d);
812  }
813 }
814 
816  double value)
817 {
818  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
819 
820  for (int i = 0; i < ess_dofs.Size(); i++)
821  if (ess_dofs[i] < 0)
822  {
823  mat -> EliminateRowColDiag (i, value);
824  }
825 }
826 
828  const Array<int> &vdofs, const Vector &x, Vector &b)
829 {
830  mat_e->AddMult(x, b, -1.);
831  mat->PartMult(vdofs, x, b);
832 }
833 
835 {
836  bool full_update;
837 
838  if (nfes && nfes != fes)
839  {
840  full_update = true;
841  fes = nfes;
842  }
843  else
844  {
845  // Check for different size (e.g. assembled form on non-conforming space)
846  // or different sequence number.
847  full_update = (fes->GetVSize() != Height() ||
848  sequence < fes->GetSequence());
849  }
850 
851  delete mat_e;
852  mat_e = NULL;
854  delete static_cond;
855  static_cond = NULL;
856 
857  if (full_update)
858  {
859  delete mat;
860  mat = NULL;
861  delete hybridization;
862  hybridization = NULL;
863  sequence = fes->GetSequence();
864  }
865  else
866  {
867  if (mat) { *mat = 0.0; }
868  if (hybridization) { hybridization->Reset(); }
869  }
870 
871  height = width = fes->GetVSize();
872 }
873 
875 {
876  delete mat_e;
877  delete mat;
878  delete element_matrices;
879  delete static_cond;
880  delete hybridization;
881 
882  if (!extern_bfs)
883  {
884  int k;
885  for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
886  for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
887  for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
888  for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
889  }
890 }
891 
892 
894  FiniteElementSpace *te_fes)
895  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
896 {
897  trial_fes = tr_fes;
898  test_fes = te_fes;
899  mat = NULL;
900 }
901 
902 double & MixedBilinearForm::Elem (int i, int j)
903 {
904  return (*mat)(i, j);
905 }
906 
907 const double & MixedBilinearForm::Elem (int i, int j) const
908 {
909  return (*mat)(i, j);
910 }
911 
912 void MixedBilinearForm::Mult (const Vector & x, Vector & y) const
913 {
914  mat -> Mult (x, y);
915 }
916 
918  const double a) const
919 {
920  mat -> AddMult (x, y, a);
921 }
922 
924  const double a) const
925 {
926  mat -> AddMultTranspose (x, y, a);
927 }
928 
930 {
931  return mat -> Inverse ();
932 }
933 
934 void MixedBilinearForm::Finalize (int skip_zeros)
935 {
936  mat -> Finalize (skip_zeros);
937 }
938 
940 {
941  MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
943  "MixedBilinearForm::GetBlocks: both trial and test spaces "
944  "must use Ordering::byNODES!");
945 
946  blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
947 
948  mat->GetBlocks(blocks);
949 }
950 
952 {
953  dom.Append (bfi);
954 }
955 
957 {
958  bdr.Append (bfi);
959 }
960 
962 {
963  skt.Append (bfi);
964 }
965 
966 void MixedBilinearForm::Assemble (int skip_zeros)
967 {
968  int i, k;
969  Array<int> tr_vdofs, te_vdofs;
970  ElementTransformation *eltrans;
971  DenseMatrix elemmat;
972 
973  Mesh *mesh = test_fes -> GetMesh();
974 
975  if (mat == NULL)
976  {
977  mat = new SparseMatrix(height, width);
978  }
979 
980  if (dom.Size())
981  {
982  for (i = 0; i < test_fes -> GetNE(); i++)
983  {
984  trial_fes -> GetElementVDofs (i, tr_vdofs);
985  test_fes -> GetElementVDofs (i, te_vdofs);
986  eltrans = test_fes -> GetElementTransformation (i);
987  for (k = 0; k < dom.Size(); k++)
988  {
989  dom[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
990  *test_fes -> GetFE(i),
991  *eltrans, elemmat);
992  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
993  }
994  }
995  }
996 
997  if (bdr.Size())
998  {
999  for (i = 0; i < test_fes -> GetNBE(); i++)
1000  {
1001  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1002  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1003  eltrans = test_fes -> GetBdrElementTransformation (i);
1004  for (k = 0; k < bdr.Size(); k++)
1005  {
1006  bdr[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1007  *test_fes -> GetBE(i),
1008  *eltrans, elemmat);
1009  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1010  }
1011  }
1012  }
1013 
1014  if (skt.Size())
1015  {
1017  Array<int> te_vdofs2;
1018  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1019 
1020  int nfaces = mesh->GetNumFaces();
1021  for (i = 0; i < nfaces; i++)
1022  {
1023  ftr = mesh->GetFaceElementTransformations(i);
1024  trial_fes->GetFaceVDofs(i, tr_vdofs);
1025  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1026  trial_face_fe = trial_fes->GetFaceElement(i);
1027  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1028  if (ftr->Elem2No >= 0)
1029  {
1030  test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
1031  te_vdofs.Append(te_vdofs2);
1032  test_fe2 = test_fes->GetFE(ftr->Elem2No);
1033  }
1034  else
1035  {
1036  // The test_fe2 object is really a dummy and not used on the
1037  // boundaries, but we can't dereference a NULL pointer, and we don't
1038  // want to actually make a fake element.
1039  test_fe2 = test_fe1;
1040  }
1041  for (int k = 0; k < skt.Size(); k++)
1042  {
1043  skt[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1044  *ftr, elemmat);
1045  mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
1046  }
1047  }
1048  }
1049 }
1050 
1052 {
1053  Finalize();
1054 
1056  if (P2)
1057  {
1058  SparseMatrix *R = Transpose(*P2);
1059  SparseMatrix *RA = mfem::Mult(*R, *mat);
1060  delete R;
1061  delete mat;
1062  mat = RA;
1063  }
1064 
1066  if (P1)
1067  {
1068  SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1069  delete mat;
1070  mat = RAP;
1071  }
1072 
1073  height = mat->Height();
1074  width = mat->Width();
1075 }
1076 
1078  Array<int> &bdr_attr_is_ess, Vector &sol, Vector &rhs )
1079 {
1080  int i, j, k;
1081  Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1082 
1083  cols_marker = 0;
1084  for (i = 0; i < trial_fes -> GetNBE(); i++)
1085  if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1086  {
1087  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1088  for (j = 0; j < tr_vdofs.Size(); j++)
1089  {
1090  if ( (k = tr_vdofs[j]) < 0 )
1091  {
1092  k = -1-k;
1093  }
1094  cols_marker[k] = 1;
1095  }
1096  }
1097  mat -> EliminateCols (cols_marker, &sol, &rhs);
1098 }
1099 
1101  Array<int> &marked_vdofs, Vector &sol, Vector &rhs)
1102 {
1103  mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1104 }
1105 
1107 {
1108  int i, j, k;
1109  Array<int> te_vdofs;
1110 
1111  for (i = 0; i < test_fes -> GetNBE(); i++)
1112  if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1113  {
1114  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1115  for (j = 0; j < te_vdofs.Size(); j++)
1116  {
1117  if ( (k = te_vdofs[j]) < 0 )
1118  {
1119  k = -1-k;
1120  }
1121  mat -> EliminateRow (k);
1122  }
1123  }
1124 }
1125 
1127 {
1128  delete mat;
1129  mat = NULL;
1130  height = test_fes->GetVSize();
1131  width = trial_fes->GetVSize();
1132 }
1133 
1135 {
1136  int i;
1137 
1138  if (mat) { delete mat; }
1139  for (i = 0; i < dom.Size(); i++) { delete dom[i]; }
1140  for (i = 0; i < bdr.Size(); i++) { delete bdr[i]; }
1141  for (i = 0; i < skt.Size(); i++) { delete skt[i]; }
1142 }
1143 
1144 
1146 {
1147  Array<int> dom_vdofs, ran_vdofs;
1149  const FiniteElement *dom_fe, *ran_fe;
1150  DenseMatrix totelmat, elmat;
1151 
1152  if (mat == NULL)
1153  {
1154  mat = new SparseMatrix(height, width);
1155  }
1156 
1157  if (dom.Size() > 0)
1158  {
1159  for (int i = 0; i < test_fes->GetNE(); i++)
1160  {
1161  trial_fes->GetElementVDofs(i, dom_vdofs);
1162  test_fes->GetElementVDofs(i, ran_vdofs);
1164  dom_fe = trial_fes->GetFE(i);
1165  ran_fe = test_fes->GetFE(i);
1166 
1167  dom[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1168  for (int j = 1; j < dom.Size(); j++)
1169  {
1170  dom[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1171  totelmat += elmat;
1172  }
1173  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1174  }
1175  }
1176 
1177  if (skt.Size())
1178  {
1179  const int nfaces = test_fes->GetMesh()->GetNumFaces();
1180  for (int i = 0; i < nfaces; i++)
1181  {
1182  trial_fes->GetFaceVDofs(i, dom_vdofs);
1183  test_fes->GetFaceVDofs(i, ran_vdofs);
1185  dom_fe = trial_fes->GetFaceElement(i);
1186  ran_fe = test_fes->GetFaceElement(i);
1187 
1188  skt[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1189  for (int j = 1; j < skt.Size(); j++)
1190  {
1191  skt[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1192  totelmat += elmat;
1193  }
1194  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1195  }
1196  }
1197 }
1198 
1199 }
Abstract class for Finite Elements.
Definition: fe.hpp:47
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:175
int Size() const
Logical size of the array.
Definition: array.hpp:110
int GetVSize() const
Definition: fespace.hpp:163
Array< BilinearFormIntegrator * > * GetBBFI()
int * GetJ()
Definition: table.hpp:111
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:855
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:101
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Array< BilinearFormIntegrator * > skt
const SparseMatrix * GetConformingProlongation() const
Definition: fespace.cpp:725
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:1476
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:278
Array< BilinearFormIntegrator * > * GetDBFI()
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:187
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition: table.cpp:188
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Definition: fespace.cpp:133
Array< BilinearFormIntegrator * > dbfi
Set of Domain Integrators to be applied.
Array< BilinearFormIntegrator * > * GetBFBFI()
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:633
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Definition: fespace.cpp:258
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:481
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
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
Abstract data type for matrix inverse.
Definition: matrix.hpp:58
void EliminateVDofs(const Array< int > &vdofs, Vector &sol, Vector &rhs, int d=0)
Eliminate the given vdofs. NOTE: here, vdofs is a list of DOFs.
void FormSystemMatrix(const Array< int > &ess_tdof_list, SparseMatrix &A)
Form the linear system matrix A, see FormLinearSystem for details.
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:144
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1432
void EliminateEssentialBCFromDofs(const Array< int > &ess_dofs, Vector &sol, Vector &rhs, int d=0)
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
Array< BilinearFormIntegrator * > * GetFBFI()
double * GetData() const
Definition: vector.hpp:121
MixedBilinearForm(FiniteElementSpace *tr_fes, FiniteElementSpace *te_fes)
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:633
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Definition: sparsemat.cpp:825
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:3334
void EliminateEssentialBCFromTrialDofs(Array< int > &marked_vdofs, Vector &sol, Vector &rhs)
Array< BilinearFormIntegrator * > bdr
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
const SparseMatrix * GetConformingRestriction() const
Definition: fespace.cpp:732
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:184
void SetSize(int m, int n)
Definition: array.hpp:284
void Finalize()
Finalize the construction of the hybridized matrix.
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
Data type sparse matrix.
Definition: sparsemat.hpp:38
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
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:548
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:136
SparseMatrix * mat
Sparse matrix to be associated with the form.
bool areColumnsSorted() const
Definition: sparsemat.hpp:284
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
Definition: mesh.cpp:367
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
Abstract data type matrix.
Definition: matrix.hpp:27
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:108
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:408
double * GetData(int k)
Definition: densemat.hpp:692
virtual void Update(FiniteElementSpace *nfes=NULL)
void EliminateTrialDofs(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs)
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:153
virtual const SparseMatrix * GetRestrictionMatrix()
Definition: fespace.hpp:149
void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi)
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1850
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:131
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:601
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1926
int SizeI() const
Definition: densemat.hpp:658
void ComputeElementMatrix(int i, DenseMatrix &elmat)
bool Finalized() const
Definition: sparsemat.hpp:283
SparseMatrix * mat_e
Matrix used to eliminate b.c.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i'th element.
Definition: fespace.hpp:205
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:22
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:172
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.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:122
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:475
Array< BilinearFormIntegrator * > bfbfi
Set of boundary face Integrators to be applied.
Array< BilinearFormIntegrator * > bbfi
Set of Boundary Integrators to be applied.
DenseTensor * element_matrices
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:235
void mfem_error(const char *msg)
Definition: error.cpp:107
void ClearExternalData()
Definition: densemat.hpp:76
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Array< BilinearFormIntegrator * > dom
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:227
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, int keep_diagonal)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Definition: staticcond.cpp:288
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
Table * GetFaceToElementTable() const
Definition: mesh.cpp:3840
FiniteElementSpace * test_fes
int SizeJ() const
Definition: densemat.hpp:659
virtual void EliminateTestDofs(Array< int > &bdr_attr_is_ess)
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)
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator.
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator.
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:552
FiniteElementSpace * fes
FE space on which the form lives.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void ComputeElementMatrices()
Compute and store internally all element matrices.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i'th element.
Definition: fespace.cpp:1171
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'th face element (2D and 3D).
Definition: fespace.cpp:145
FiniteElementSpace * trial_fes
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual ~BilinearForm()
Destroys bilinear form.
Vector data type.
Definition: vector.hpp:41
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
Definition: staticcond.hpp:152
virtual MatrixInverse * Inverse() const
Returns a pointer to (an approximation) of the matrix inverse.
Hybridization * hybridization
void ReduceRHS(const Vector &b, Vector &b_r) const
int * GetI()
Definition: table.hpp:110
const Table & GetElementToDofTable() const
Definition: fespace.hpp:288
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
Array< Array< int > * > bfbfi_marker
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i'th boundary element.
Definition: fespace.cpp:1406
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
Definition: fespace.cpp:139
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:86
void FreeElementMatrices()
Free the memory used by the element matrices.
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:376
void EliminateEssentialBC(const Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
void Init(bool symmetric, bool block_diagonal)
Definition: staticcond.cpp:134
virtual void Assemble(int skip_zeros=1)
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:133
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:634
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 AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Definition: sparsemat.cpp:597