MFEM  v3.4 Finite element discretization library
bilinearform.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
11
12 // Implementation of class BilinearForm
13
14 #include "fem.hpp"
15 #include <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;
77 }
78
80  : Matrix (f->GetVSize())
81 {
82  int i;
84
85  fes = f;
86  sequence = f->GetSequence();
87  mat_e = NULL;
88  extern_bfs = 1;
89  element_matrices = NULL;
90  static_cond = NULL;
91  hybridization = NULL;
94
95  bfi = bf->GetDBFI();
96  dbfi.SetSize (bfi->Size());
97  for (i = 0; i < bfi->Size(); i++)
98  {
99  dbfi[i] = (*bfi)[i];
100  }
101
102  bfi = bf->GetBBFI();
103  bbfi.SetSize (bfi->Size());
104  for (i = 0; i < bfi->Size(); i++)
105  {
106  bbfi[i] = (*bfi)[i];
107  }
108
109  bfi = bf->GetFBFI();
110  fbfi.SetSize (bfi->Size());
111  for (i = 0; i < bfi->Size(); i++)
112  {
113  fbfi[i] = (*bfi)[i];
114  }
115
116  bfi = bf->GetBFBFI();
117  bfbfi.SetSize (bfi->Size());
118  for (i = 0; i < bfi->Size(); i++)
119  {
120  bfbfi[i] = (*bfi)[i];
121  }
122
123  AllocMat();
124 }
125
127 {
128  delete static_cond;
131  {
132  bool symmetric = false; // TODO
133  bool block_diagonal = false; // TODO
134  static_cond->Init(symmetric, block_diagonal);
135  }
136  else
137  {
138  delete static_cond;
139  static_cond = NULL;
140  }
141 }
142
144  BilinearFormIntegrator *constr_integ,
145  const Array<int> &ess_tdof_list)
146 {
147  delete hybridization;
148  hybridization = new Hybridization(fes, constr_space);
149  hybridization->SetConstraintIntegrator(constr_integ);
150  hybridization->Init(ess_tdof_list);
151 }
152
153 void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
154 {
155  if (static_cond) { return; }
156
157  if (mat)
158  {
159  if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
160  {
161  return; // mat is already using the given sparsity
162  }
163  delete mat;
164  }
165  height = width = fes->GetVSize();
166  mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
167 }
168
170 {
171  MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
172  "invalid matrix A dimensions: "
173  << A.Height() << " x " << A.Width());
174  MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
175
176  UseSparsity(A.GetI(), A.GetJ(), A.areColumnsSorted());
177 }
178
179 double& BilinearForm::Elem (int i, int j)
180 {
181  return mat -> Elem(i,j);
182 }
183
184 const double& BilinearForm::Elem (int i, int j) const
185 {
186  return mat -> Elem(i,j);
187 }
188
190 {
191  return mat -> Inverse();
192 }
193
194 void BilinearForm::Finalize (int skip_zeros)
195 {
196  if (!static_cond) { mat->Finalize(skip_zeros); }
197  if (mat_e) { mat_e->Finalize(skip_zeros); }
198  if (static_cond) { static_cond->Finalize(); }
200 }
201
203 {
204  dbfi.Append (bfi);
205 }
206
208 {
209  bbfi.Append (bfi);
210  bbfi_marker.Append(NULL); // NULL marker means apply everywhere
211 }
212
214  Array<int> &bdr_marker)
215 {
216  bbfi.Append (bfi);
217  bbfi_marker.Append(&bdr_marker);
218 }
219
221 {
222  fbfi.Append (bfi);
223 }
224
226 {
227  bfbfi.Append(bfi);
228  bfbfi_marker.Append(NULL); // NULL marker means apply everywhere
229 }
230
232  Array<int> &bdr_marker)
233 {
234  bfbfi.Append(bfi);
235  bfbfi_marker.Append(&bdr_marker);
236 }
237
239 {
240  if (element_matrices)
241  {
243  elmat = element_matrices->GetData(i);
244  return;
245  }
246
247  if (dbfi.Size())
248  {
249  const FiniteElement &fe = *fes->GetFE(i);
251  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
252  for (int k = 1; k < dbfi.Size(); k++)
253  {
254  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
255  elmat += elemmat;
256  }
257  }
258  else
259  {
261  elmat.SetSize(vdofs.Size());
262  elmat = 0.0;
263  }
264 }
265
267  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
268 {
269  fes->GetElementVDofs(i, vdofs);
270  if (static_cond)
271  {
272  static_cond->AssembleMatrix(i, elmat);
273  }
274  else
275  {
276  if (mat == NULL)
277  {
278  AllocMat();
279  }
281  if (hybridization)
282  {
283  hybridization->AssembleMatrix(i, elmat);
284  }
285  }
286 }
287
289  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
290 {
291  fes->GetBdrElementVDofs(i, vdofs);
292  if (static_cond)
293  {
294  static_cond->AssembleBdrMatrix(i, elmat);
295  }
296  else
297  {
298  if (mat == NULL)
299  {
300  AllocMat();
301  }
303  if (hybridization)
304  {
305  hybridization->AssembleBdrMatrix(i, elmat);
306  }
307  }
308 }
309
310 void BilinearForm::Assemble (int skip_zeros)
311 {
312  ElementTransformation *eltrans;
313  Mesh *mesh = fes -> GetMesh();
314  DenseMatrix elmat, *elmat_p;
315
316  int i;
317
318  if (mat == NULL)
319  {
320  AllocMat();
321  }
322
323 #ifdef MFEM_USE_OPENMP
324  int free_element_matrices = 0;
325  if (!element_matrices)
326  {
328  free_element_matrices = 1;
329  }
330 #endif
331
332  if (dbfi.Size())
333  {
334  for (i = 0; i < fes -> GetNE(); i++)
335  {
337  if (element_matrices)
338  {
339  elmat_p = &(*element_matrices)(i);
340  }
341  else
342  {
343  const FiniteElement &fe = *fes->GetFE(i);
344  eltrans = fes->GetElementTransformation(i);
345  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
346  for (int k = 1; k < dbfi.Size(); k++)
347  {
348  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
349  elmat += elemmat;
350  }
351  elmat_p = &elmat;
352  }
353  if (static_cond)
354  {
355  static_cond->AssembleMatrix(i, *elmat_p);
356  }
357  else
358  {
360  if (hybridization)
361  {
362  hybridization->AssembleMatrix(i, *elmat_p);
363  }
364  }
365  }
366  }
367
368  if (bbfi.Size())
369  {
370  // Which boundary attributes need to be processed?
371  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
372  mesh->bdr_attributes.Max() : 0);
373  bdr_attr_marker = 0;
374  for (int k = 0; k < bbfi.Size(); k++)
375  {
376  if (bbfi_marker[k] == NULL)
377  {
378  bdr_attr_marker = 1;
379  break;
380  }
381  Array<int> &bdr_marker = *bbfi_marker[k];
382  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
383  "invalid boundary marker for boundary integrator #"
384  << k << ", counting from zero");
385  for (int i = 0; i < bdr_attr_marker.Size(); i++)
386  {
387  bdr_attr_marker[i] |= bdr_marker[i];
388  }
389  }
390
391  for (i = 0; i < fes -> GetNBE(); i++)
392  {
393  const int bdr_attr = mesh->GetBdrAttribute(i);
394  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
395
396  const FiniteElement &be = *fes->GetBE(i);
397  fes -> GetBdrElementVDofs (i, vdofs);
398  eltrans = fes -> GetBdrElementTransformation (i);
399  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
400  for (int k = 1; k < bbfi.Size(); k++)
401  {
402  if (bbfi_marker[k] &&
403  (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
404
405  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
406  elmat += elemmat;
407  }
408  if (!static_cond)
409  {
411  if (hybridization)
412  {
413  hybridization->AssembleBdrMatrix(i, elmat);
414  }
415  }
416  else
417  {
418  static_cond->AssembleBdrMatrix(i, elmat);
419  }
420  }
421  }
422
423  if (fbfi.Size())
424  {
426  Array<int> vdofs2;
427
428  int nfaces = mesh->GetNumFaces();
429  for (i = 0; i < nfaces; i++)
430  {
431  tr = mesh -> GetInteriorFaceTransformations (i);
432  if (tr != NULL)
433  {
434  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
435  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
436  vdofs.Append (vdofs2);
437  for (int k = 0; k < fbfi.Size(); k++)
438  {
439  fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
440  *fes -> GetFE (tr -> Elem2No),
441  *tr, elemmat);
442  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
443  }
444  }
445  }
446  }
447
448  if (bfbfi.Size())
449  {
451  const FiniteElement *fe1, *fe2;
452
453  // Which boundary attributes need to be processed?
454  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
455  mesh->bdr_attributes.Max() : 0);
456  bdr_attr_marker = 0;
457  for (int k = 0; k < bfbfi.Size(); k++)
458  {
459  if (bfbfi_marker[k] == NULL)
460  {
461  bdr_attr_marker = 1;
462  break;
463  }
464  Array<int> &bdr_marker = *bfbfi_marker[k];
465  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
466  "invalid boundary marker for boundary face integrator #"
467  << k << ", counting from zero");
468  for (int i = 0; i < bdr_attr_marker.Size(); i++)
469  {
470  bdr_attr_marker[i] |= bdr_marker[i];
471  }
472  }
473
474  for (i = 0; i < fes -> GetNBE(); i++)
475  {
476  const int bdr_attr = mesh->GetBdrAttribute(i);
477  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
478
479  tr = mesh -> GetBdrFaceTransformations (i);
480  if (tr != NULL)
481  {
482  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
483  fe1 = fes -> GetFE (tr -> Elem1No);
484  // The fe2 object is really a dummy and not used on the boundaries,
485  // but we can't dereference a NULL pointer, and we don't want to
486  // actually make a fake element.
487  fe2 = fe1;
488  for (int k = 0; k < bfbfi.Size(); k++)
489  {
490  if (bfbfi_marker[k] &&
491  (*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
492
493  bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
494  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
495  }
496  }
497  }
498  }
499
500 #ifdef MFEM_USE_OPENMP
501  if (free_element_matrices)
502  {
504  }
505 #endif
506 }
507
509 {
510  // Do not remove zero entries to preserve the symmetric structure of the
511  // matrix which in turn will give rise to symmetric structure in the new
512  // matrix. This ensures that subsequent calls to EliminateRowCol will work
513  // correctly.
514  Finalize(0);
515  MFEM_ASSERT(mat, "the BilinearForm is not assembled");
516
518  if (!P) { return; } // conforming mesh
519
520  SparseMatrix *R = Transpose(*P);
521  SparseMatrix *RA = mfem::Mult(*R, *mat);
522  delete mat;
523  if (mat_e)
524  {
525  SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
526  delete mat_e;
527  mat_e = RAe;
528  }
529  delete R;
530  mat = mfem::Mult(*RA, *P);
531  delete RA;
532  if (mat_e)
533  {
534  SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
535  delete mat_e;
536  mat_e = RAeP;
537  }
538
539  height = mat->Height();
540  width = mat->Width();
541 }
542
543 void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
544  Vector &x, Vector &b,
545  SparseMatrix &A, Vector &X, Vector &B,
546  int copy_interior)
547 {
549
550  FormSystemMatrix(ess_tdof_list, A);
551
552  // Transform the system and perform the elimination in B, based on the
553  // essential BC values from x. Restrict the BC part of x in X, and set the
554  // non-BC part to zero. Since there is no good initial guess for the Lagrange
555  // multipliers, set X = 0.0 for hybridization.
556  if (static_cond)
557  {
558  // Schur complement reduction to the exposed dofs
559  static_cond->ReduceSystem(x, b, X, B, copy_interior);
560  }
561  else if (!P) // conforming space
562  {
563  if (hybridization)
564  {
565  // Reduction to the Lagrange multipliers system
566  EliminateVDofsInRHS(ess_tdof_list, x, b);
567  hybridization->ReduceRHS(b, B);
568  X.SetSize(B.Size());
569  X = 0.0;
570  }
571  else
572  {
573  // A, X and B point to the same data as mat, x and b
574  EliminateVDofsInRHS(ess_tdof_list, x, b);
575  X.NewDataAndSize(x.GetData(), x.Size());
576  B.NewDataAndSize(b.GetData(), b.Size());
577  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
578  }
579  }
580  else // non-conforming space
581  {
582  if (hybridization)
583  {
584  // Reduction to the Lagrange multipliers system
586  Vector conf_b(P->Width()), conf_x(P->Width());
587  P->MultTranspose(b, conf_b);
588  R->Mult(x, conf_x);
589  EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
590  R->MultTranspose(conf_b, b); // store eliminated rhs in b
591  hybridization->ReduceRHS(conf_b, B);
592  X.SetSize(B.Size());
593  X = 0.0;
594  }
595  else
596  {
597  // Variational restriction with P
599  B.SetSize(P->Width());
600  P->MultTranspose(b, B);
601  X.SetSize(R->Height());
602  R->Mult(x, X);
603  EliminateVDofsInRHS(ess_tdof_list, X, B);
604  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
605  }
606  }
607 }
608
609 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
610  SparseMatrix &A)
611 {
612  // Finish the matrix assembly and perform BC elimination, storing the
613  // eliminated part of the matrix.
614  if (static_cond)
615  {
617  {
618  static_cond->SetEssentialTrueDofs(ess_tdof_list);
619  static_cond->Finalize(); // finalize Schur complement (to true dofs)
621  static_cond->Finalize(); // finalize eliminated part
622  }
624  }
625  else
626  {
627  if (!mat_e)
628  {
630  if (P) { ConformingAssemble(); }
631  EliminateVDofs(ess_tdof_list, diag_policy);
632  const int remove_zeros = 0;
633  Finalize(remove_zeros);
634  }
635  if (hybridization)
636  {
638  }
639  else
640  {
641  A.MakeRef(*mat);
642  }
643  }
644 }
645
647  const Vector &b, Vector &x)
648 {
650  if (!P) // conforming space
651  {
652  if (static_cond)
653  {
654  // Private dofs back solve
655  static_cond->ComputeSolution(b, X, x);
656  }
657  else if (hybridization)
658  {
659  // Primal unknowns recovery
660  hybridization->ComputeSolution(b, X, x);
661  }
662  else
663  {
664  // X and x point to the same data
665  }
666  }
667  else // non-conforming space
668  {
669  if (static_cond)
670  {
671  // Private dofs back solve
672  static_cond->ComputeSolution(b, X, x);
673  }
674  else if (hybridization)
675  {
676  // Primal unknowns recovery
677  Vector conf_b(P->Width()), conf_x(P->Width());
678  P->MultTranspose(b, conf_b);
680  R->Mult(x, conf_x); // get essential b.c. from x
681  hybridization->ComputeSolution(conf_b, X, conf_x);
682  x.SetSize(P->Height());
683  P->Mult(conf_x, x);
684  }
685  else
686  {
687  // Apply conforming prolongation
688  x.SetSize(P->Height());
689  P->Mult(X, x);
690  }
691  }
692 }
693
695 {
696  if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
697  {
698  return;
699  }
700
701  int num_elements = fes->GetNE();
702  int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
703
704  element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
705  num_elements);
706
707  DenseMatrix tmp;
709
710 #ifdef MFEM_USE_OPENMP
711  #pragma omp parallel for private(tmp,eltrans)
712 #endif
713  for (int i = 0; i < num_elements; i++)
714  {
716  num_dofs_per_el, num_dofs_per_el);
717  const FiniteElement &fe = *fes->GetFE(i);
718 #ifdef MFEM_DEBUG
719  if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
720  mfem_error("BilinearForm::ComputeElementMatrices:"
721  " all elements must have same number of dofs");
722 #endif
723  fes->GetElementTransformation(i, &eltrans);
724
725  dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
726  for (int k = 1; k < dbfi.Size(); k++)
727  {
728  // note: some integrators may not be thread-safe
729  dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
730  elmat += tmp;
731  }
732  elmat.ClearExternalData();
733  }
734 }
735
736 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
737  const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy)
738 {
739  Array<int> ess_dofs, conf_ess_dofs;
740  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
741
742  if (fes->GetVSize() == height)
743  {
744  EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
745  }
746  else
747  {
748  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
749  EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
750  }
751 }
752
753 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
754  DiagonalPolicy dpolicy)
755 {
756  Array<int> ess_dofs, conf_ess_dofs;
757  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
758
759  if (fes->GetVSize() == height)
760  {
761  EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
762  }
763  else
764  {
765  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
766  EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
767  }
768 }
769
771  double value)
772 {
773  Array<int> ess_dofs, conf_ess_dofs;
774  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
775
776  if (fes->GetVSize() == height)
777  {
778  EliminateEssentialBCFromDofsDiag(ess_dofs, value);
779  }
780  else
781  {
782  fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
783  EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
784  }
785 }
786
788  const Vector &sol, Vector &rhs,
789  DiagonalPolicy dpolicy)
790 {
791  for (int i = 0; i < vdofs.Size(); i++)
792  {
793  int vdof = vdofs[i];
794  if ( vdof >= 0 )
795  {
796  mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
797  }
798  else
799  {
800  mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
801  }
802  }
803 }
804
806  DiagonalPolicy dpolicy)
807 {
808  if (mat_e == NULL)
809  {
810  mat_e = new SparseMatrix(height);
811  }
812
813  for (int i = 0; i < vdofs.Size(); i++)
814  {
815  int vdof = vdofs[i];
816  if ( vdof >= 0 )
817  {
818  mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
819  }
820  else
821  {
822  mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
823  }
824  }
825 }
826
828  const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
829  DiagonalPolicy dpolicy)
830 {
831  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
832  MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
833  MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
834
835  for (int i = 0; i < ess_dofs.Size(); i++)
836  if (ess_dofs[i] < 0)
837  {
838  mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
839  }
840 }
841
843  DiagonalPolicy dpolicy)
844 {
845  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
846
847  for (int i = 0; i < ess_dofs.Size(); i++)
848  if (ess_dofs[i] < 0)
849  {
850  mat -> EliminateRowCol (i, dpolicy);
851  }
852 }
853
855  double value)
856 {
857  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
858
859  for (int i = 0; i < ess_dofs.Size(); i++)
860  if (ess_dofs[i] < 0)
861  {
862  mat -> EliminateRowColDiag (i, value);
863  }
864 }
865
867  const Array<int> &vdofs, const Vector &x, Vector &b)
868 {
870  mat->PartMult(vdofs, x, b);
871 }
872
874 {
875  bool full_update;
876
877  if (nfes && nfes != fes)
878  {
879  full_update = true;
880  fes = nfes;
881  }
882  else
883  {
884  // Check for different size (e.g. assembled form on non-conforming space)
885  // or different sequence number.
886  full_update = (fes->GetVSize() != Height() ||
887  sequence < fes->GetSequence());
888  }
889
890  delete mat_e;
891  mat_e = NULL;
893  delete static_cond;
894  static_cond = NULL;
895
896  if (full_update)
897  {
898  delete mat;
899  mat = NULL;
900  delete hybridization;
901  hybridization = NULL;
902  sequence = fes->GetSequence();
903  }
904  else
905  {
906  if (mat) { *mat = 0.0; }
907  if (hybridization) { hybridization->Reset(); }
908  }
909
910  height = width = fes->GetVSize();
911 }
912
913 void BilinearForm::SetDiagonalPolicy(DiagonalPolicy policy)
914 {
915  diag_policy = policy;
916 }
917
919 {
920  delete mat_e;
921  delete mat;
922  delete element_matrices;
923  delete static_cond;
924  delete hybridization;
925
926  if (!extern_bfs)
927  {
928  int k;
929  for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
930  for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
931  for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
932  for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
933  }
934 }
935
936
938  FiniteElementSpace *te_fes)
939  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
940 {
941  trial_fes = tr_fes;
942  test_fes = te_fes;
943  mat = NULL;
944 }
945
946 double & MixedBilinearForm::Elem (int i, int j)
947 {
948  return (*mat)(i, j);
949 }
950
951 const double & MixedBilinearForm::Elem (int i, int j) const
952 {
953  return (*mat)(i, j);
954 }
955
956 void MixedBilinearForm::Mult (const Vector & x, Vector & y) const
957 {
958  mat -> Mult (x, y);
959 }
960
962  const double a) const
963 {
964  mat -> AddMult (x, y, a);
965 }
966
968  const double a) const
969 {
970  mat -> AddMultTranspose (x, y, a);
971 }
972
974 {
975  return mat -> Inverse ();
976 }
977
978 void MixedBilinearForm::Finalize (int skip_zeros)
979 {
980  mat -> Finalize (skip_zeros);
981 }
982
984 {
985  MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
987  "MixedBilinearForm::GetBlocks: both trial and test spaces "
988  "must use Ordering::byNODES!");
989
990  blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
991
992  mat->GetBlocks(blocks);
993 }
994
996 {
997  dom.Append (bfi);
998 }
999
1001 {
1002  bdr.Append (bfi);
1003 }
1004
1006 {
1007  skt.Append (bfi);
1008 }
1009
1010 void MixedBilinearForm::Assemble (int skip_zeros)
1011 {
1012  int i, k;
1013  Array<int> tr_vdofs, te_vdofs;
1014  ElementTransformation *eltrans;
1015  DenseMatrix elemmat;
1016
1017  Mesh *mesh = test_fes -> GetMesh();
1018
1019  if (mat == NULL)
1020  {
1021  mat = new SparseMatrix(height, width);
1022  }
1023
1024  if (dom.Size())
1025  {
1026  for (i = 0; i < test_fes -> GetNE(); i++)
1027  {
1028  trial_fes -> GetElementVDofs (i, tr_vdofs);
1029  test_fes -> GetElementVDofs (i, te_vdofs);
1030  eltrans = test_fes -> GetElementTransformation (i);
1031  for (k = 0; k < dom.Size(); k++)
1032  {
1033  dom[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1034  *test_fes -> GetFE(i),
1035  *eltrans, elemmat);
1036  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1037  }
1038  }
1039  }
1040
1041  if (bdr.Size())
1042  {
1043  for (i = 0; i < test_fes -> GetNBE(); i++)
1044  {
1045  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1046  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1047  eltrans = test_fes -> GetBdrElementTransformation (i);
1048  for (k = 0; k < bdr.Size(); k++)
1049  {
1050  bdr[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1051  *test_fes -> GetBE(i),
1052  *eltrans, elemmat);
1053  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1054  }
1055  }
1056  }
1057
1058  if (skt.Size())
1059  {
1061  Array<int> te_vdofs2;
1062  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1063
1064  int nfaces = mesh->GetNumFaces();
1065  for (i = 0; i < nfaces; i++)
1066  {
1067  ftr = mesh->GetFaceElementTransformations(i);
1068  trial_fes->GetFaceVDofs(i, tr_vdofs);
1069  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1070  trial_face_fe = trial_fes->GetFaceElement(i);
1071  test_fe1 = test_fes->GetFE(ftr->Elem1No);
1072  if (ftr->Elem2No >= 0)
1073  {
1074  test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
1075  te_vdofs.Append(te_vdofs2);
1076  test_fe2 = test_fes->GetFE(ftr->Elem2No);
1077  }
1078  else
1079  {
1080  // The test_fe2 object is really a dummy and not used on the
1081  // boundaries, but we can't dereference a NULL pointer, and we don't
1082  // want to actually make a fake element.
1083  test_fe2 = test_fe1;
1084  }
1085  for (int k = 0; k < skt.Size(); k++)
1086  {
1087  skt[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
1088  *ftr, elemmat);
1090  }
1091  }
1092  }
1093 }
1094
1096 {
1097  Finalize();
1098
1100  if (P2)
1101  {
1102  SparseMatrix *R = Transpose(*P2);
1103  SparseMatrix *RA = mfem::Mult(*R, *mat);
1104  delete R;
1105  delete mat;
1106  mat = RA;
1107  }
1108
1110  if (P1)
1111  {
1112  SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1113  delete mat;
1114  mat = RAP;
1115  }
1116
1117  height = mat->Height();
1118  width = mat->Width();
1119 }
1120
1122  Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
1123 {
1124  int i, j, k;
1125  Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1126
1127  cols_marker = 0;
1128  for (i = 0; i < trial_fes -> GetNBE(); i++)
1129  if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1130  {
1131  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1132  for (j = 0; j < tr_vdofs.Size(); j++)
1133  {
1134  if ( (k = tr_vdofs[j]) < 0 )
1135  {
1136  k = -1-k;
1137  }
1138  cols_marker[k] = 1;
1139  }
1140  }
1141  mat -> EliminateCols (cols_marker, &sol, &rhs);
1142 }
1143
1145  Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
1146 {
1147  mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1148 }
1149
1151 {
1152  int i, j, k;
1153  Array<int> te_vdofs;
1154
1155  for (i = 0; i < test_fes -> GetNBE(); i++)
1156  if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1157  {
1158  test_fes -> GetBdrElementVDofs (i, te_vdofs);
1159  for (j = 0; j < te_vdofs.Size(); j++)
1160  {
1161  if ( (k = te_vdofs[j]) < 0 )
1162  {
1163  k = -1-k;
1164  }
1165  mat -> EliminateRow (k);
1166  }
1167  }
1168 }
1169
1171 {
1172  delete mat;
1173  mat = NULL;
1174  height = test_fes->GetVSize();
1175  width = trial_fes->GetVSize();
1176 }
1177
1179 {
1180  int i;
1181
1182  if (mat) { delete mat; }
1183  for (i = 0; i < dom.Size(); i++) { delete dom[i]; }
1184  for (i = 0; i < bdr.Size(); i++) { delete bdr[i]; }
1185  for (i = 0; i < skt.Size(); i++) { delete skt[i]; }
1186 }
1187
1188
1190 {
1191  Array<int> dom_vdofs, ran_vdofs;
1193  const FiniteElement *dom_fe, *ran_fe;
1194  DenseMatrix totelmat, elmat;
1195
1196  if (mat == NULL)
1197  {
1198  mat = new SparseMatrix(height, width);
1199  }
1200
1201  if (dom.Size() > 0)
1202  {
1203  for (int i = 0; i < test_fes->GetNE(); i++)
1204  {
1205  trial_fes->GetElementVDofs(i, dom_vdofs);
1206  test_fes->GetElementVDofs(i, ran_vdofs);
1208  dom_fe = trial_fes->GetFE(i);
1209  ran_fe = test_fes->GetFE(i);
1210
1211  dom[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1212  for (int j = 1; j < dom.Size(); j++)
1213  {
1214  dom[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1215  totelmat += elmat;
1216  }
1217  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1218  }
1219  }
1220
1221  if (skt.Size())
1222  {
1223  const int nfaces = test_fes->GetMesh()->GetNumFaces();
1224  for (int i = 0; i < nfaces; i++)
1225  {
1226  trial_fes->GetFaceVDofs(i, dom_vdofs);
1227  test_fes->GetFaceVDofs(i, ran_vdofs);
1229  dom_fe = trial_fes->GetFaceElement(i);
1230  ran_fe = test_fes->GetFaceElement(i);
1231
1232  skt[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1233  for (int j = 1; j < skt.Size(); j++)
1234  {
1235  skt[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1236  totelmat += elmat;
1237  }
1238  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1239  }
1240  }
1241 }
1242
1243 }
Abstract class for Finite Elements.
Definition: fe.hpp:140
void EliminateEssentialBCFromDofs(const Array< int > &ess_dofs, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateVDofs(const Array&lt;int&gt; &amp;, const Vector &amp;, Vector &amp;, DiagonalPolicy) but here ess_...
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:265
int Size() const
Logical size of the array.
Definition: array.hpp:133
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:253
Array< BilinearFormIntegrator * > * GetBBFI()
int * GetJ()
Definition: table.hpp:114
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:868
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:108
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Array< BilinearFormIntegrator * > skt
const SparseMatrix * GetConformingProlongation() const
Definition: fespace.cpp:761
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:1558
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:305
Array< BilinearFormIntegrator * > * GetDBFI()
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:197
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition: table.cpp:198
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
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:171
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:644
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
void EliminateEssentialBCFromTrialDofs(Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Array< Array< int > * > bbfi_marker
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
void EliminateTrialDofs(Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Definition: fespace.cpp:296
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:492
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
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
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Definition: staticcond.cpp:288
Abstract data type for matrix inverse.
Definition: matrix.hpp:66
void 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:147
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1597
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
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
MixedBilinearForm(FiniteElementSpace *tr_fes, FiniteElementSpace *te_fes)
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:657
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:836
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
void EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
void EliminateEssentialBCFromDofsDiag(const Array< int > &ess_dofs, double value)
Perform elimination and set the diagonal entry to the given value.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3374
Array< BilinearFormIntegrator * > bdr
const SparseMatrix * GetConformingRestriction() const
Definition: fespace.cpp:768
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:277
void SetSize(int m, int n)
Definition: array.hpp:308
void Finalize()
Finalize the construction of the hybridized matrix.
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
Keep the diagonal value.
Definition: matrix.hpp:36
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
Data type sparse matrix.
Definition: sparsemat.hpp: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:614
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:224
SparseMatrix * mat
Sparse matrix to be associated with the form.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:146
bool areColumnsSorted() const
Definition: sparsemat.hpp:311
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
Returns the transformation defining the given face element in a user-defined variable.
Definition: mesh.cpp:391
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
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.
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:109
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:418
double * GetData(int k)
Definition: densemat.hpp:712
virtual void Update(FiniteElementSpace *nfes=NULL)
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:242
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1978
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:134
Dynamic 2D array using row-major layout.
Definition: array.hpp:289
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the &#39;dofs&#39; array to the given &#39;val&#39;.
Definition: vector.cpp:597
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2054
int SizeI() const
Definition: densemat.hpp:678
void ComputeElementMatrix(int i, DenseMatrix &elmat)
bool Finalized() const
Definition: sparsemat.hpp:310
SparseMatrix * mat_e
Matrix used to eliminate b.c.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:301
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:174
void AssembleElementMatrix(int i, const DenseMatrix &elmat, Array< int > &vdofs, int skip_zeros=1)
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:66
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:217
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:486
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 ClearExternalData()
Definition: densemat.hpp:76
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Array< BilinearFormIntegrator * > dom
DiagonalPolicy diag_policy
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:227
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
Table * GetFaceToElementTable() const
Definition: mesh.cpp:3880
FiniteElementSpace * test_fes
int SizeJ() const
Definition: densemat.hpp:679
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)
Use the stored eliminated part of the matrix (see EliminateVDofs(const Array&lt;int&gt; &amp;...
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:563
FiniteElementSpace * fes
FE space on which the form lives.
void ComputeElementMatrices()
Compute and store internally all element matrices.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1335
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, Array< int > &vdofs, int skip_zeros=1)
DenseMatrix elemmat
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th face element (2D and 3D).
Definition: fespace.cpp:183
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:48
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:113
const Table & GetElementToDofTable() const
Definition: fespace.hpp:384
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&#39;th boundary element.
Definition: fespace.cpp:1571
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:177
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.
virtual const SparseMatrix * GetRestrictionMatrix() const
Definition: fespace.hpp:238
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:517
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:136
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
Form a linear system, A X = B.
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:638
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()