MFEM  v3.1 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
31  const Table &elem_dof = fes->GetElementToDofTable();
32  Table dof_dof;
33
34  if (fbfi.Size() > 0)
35  {
36  // the sparsity pattern is defined from the map: face->element->dof
37  Table face_dof, dof_face;
38  {
39  Table *face_elem = fes->GetMesh()->GetFaceToElementTable();
40  mfem::Mult(*face_elem, elem_dof, face_dof);
41  delete face_elem;
42  }
43  Transpose(face_dof, dof_face, height);
44  mfem::Mult(dof_face, face_dof, dof_dof);
45  }
46  else
47  {
48  // the sparsity pattern is defined from the map: element->dof
49  Table dof_elem;
50  Transpose(elem_dof, dof_elem, height);
51  mfem::Mult(dof_elem, elem_dof, dof_dof);
52  }
53
54  int *I = dof_dof.GetI();
55  int *J = dof_dof.GetJ();
56  double *data = new double[I[height]];
57
58  mat = new SparseMatrix(I, J, data, height, height);
59  *mat = 0.0;
60
61  dof_dof.LoseData();
62 }
63
65  : Matrix (f->GetVSize())
66 {
67  fes = f;
68  mat = mat_e = NULL;
69  extern_bfs = 0;
70  element_matrices = NULL;
71  static_cond = NULL;
72  hybridization = NULL;
74 }
75
77  : Matrix (f->GetVSize())
78 {
79  int i;
81
82  fes = f;
83  mat_e = NULL;
84  extern_bfs = 1;
85  element_matrices = NULL;
86  static_cond = NULL;
87  hybridization = NULL;
89
90  bfi = bf->GetDBFI();
91  dbfi.SetSize (bfi->Size());
92  for (i = 0; i < bfi->Size(); i++)
93  {
94  dbfi[i] = (*bfi)[i];
95  }
96
97  bfi = bf->GetBBFI();
98  bbfi.SetSize (bfi->Size());
99  for (i = 0; i < bfi->Size(); i++)
100  {
101  bbfi[i] = (*bfi)[i];
102  }
103
104  bfi = bf->GetFBFI();
105  fbfi.SetSize (bfi->Size());
106  for (i = 0; i < bfi->Size(); i++)
107  {
108  fbfi[i] = (*bfi)[i];
109  }
110
111  bfi = bf->GetBFBFI();
112  bfbfi.SetSize (bfi->Size());
113  for (i = 0; i < bfi->Size(); i++)
114  {
115  bfbfi[i] = (*bfi)[i];
116  }
117
118  AllocMat();
119 }
120
122 {
123  delete static_cond;
126  {
127  bool symmetric = false; // TODO
128  bool block_diagonal = false; // TODO
129  static_cond->Init(symmetric, block_diagonal);
130  }
131  else
132  {
133  delete static_cond;
134  static_cond = NULL;
135  }
136 }
137
139  BilinearFormIntegrator *constr_integ,
140  const Array<int> &ess_tdof_list)
141 {
142  delete hybridization;
143  hybridization = new Hybridization(fes, constr_space);
144  hybridization->SetConstraintIntegrator(constr_integ);
145  hybridization->Init(ess_tdof_list);
146 }
147
148 double& BilinearForm::Elem (int i, int j)
149 {
150  return mat -> Elem(i,j);
151 }
152
153 const double& BilinearForm::Elem (int i, int j) const
154 {
155  return mat -> Elem(i,j);
156 }
157
158 void BilinearForm::Mult (const Vector & x, Vector & y) const
159 {
160  mat -> Mult (x, y);
161 }
162
164 {
165  return mat -> Inverse();
166 }
167
168 void BilinearForm::Finalize (int skip_zeros)
169 {
170  if (!static_cond) { mat->Finalize(skip_zeros); }
171  if (mat_e) { mat_e->Finalize(skip_zeros); }
172  if (static_cond) { static_cond->Finalize(); }
174 }
175
177 {
178  dbfi.Append (bfi);
179 }
180
182 {
183  bbfi.Append (bfi);
184 }
185
187 {
188  fbfi.Append (bfi);
189 }
190
192 {
193  bfbfi.Append (bfi);
194 }
195
197 {
198  if (element_matrices)
199  {
201  elmat = element_matrices->GetData(i);
202  return;
203  }
204
205  if (dbfi.Size())
206  {
207  const FiniteElement &fe = *fes->GetFE(i);
209  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
210  for (int k = 1; k < dbfi.Size(); k++)
211  {
212  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
213  elmat += elemmat;
214  }
215  }
216  else
217  {
219  elmat.SetSize(vdofs.Size());
220  elmat = 0.0;
221  }
222 }
223
225  int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
226 {
227  fes->GetElementVDofs(i, vdofs);
228  if (static_cond)
229  {
230  static_cond->AssembleMatrix(i, elmat);
231  }
232  else
233  {
234  if (mat == NULL)
235  {
236  AllocMat();
237  }
239  if (hybridization)
240  {
241  hybridization->AssembleMatrix(i, elmat);
242  }
243  }
244 }
245
246 void BilinearForm::Assemble (int skip_zeros)
247 {
248  ElementTransformation *eltrans;
249  Mesh *mesh = fes -> GetMesh();
250  DenseMatrix elmat, *elmat_p;
251
252  int i;
253
254  if (mat == NULL)
255  {
256  AllocMat();
257  }
258
259 #ifdef MFEM_USE_OPENMP
260  int free_element_matrices = 0;
261  if (!element_matrices)
262  {
264  free_element_matrices = 1;
265  }
266 #endif
267
268  if (dbfi.Size())
269  {
270  for (i = 0; i < fes -> GetNE(); i++)
271  {
273  if (element_matrices)
274  {
275  elmat_p = &(*element_matrices)(i);
276  }
277  else
278  {
279  const FiniteElement &fe = *fes->GetFE(i);
280  eltrans = fes->GetElementTransformation(i);
281  dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
282  for (int k = 1; k < dbfi.Size(); k++)
283  {
284  dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
285  elmat += elemmat;
286  }
287  elmat_p = &elmat;
288  }
289  if (static_cond)
290  {
291  static_cond->AssembleMatrix(i, *elmat_p);
292  }
293  else
294  {
296  if (hybridization)
297  {
298  hybridization->AssembleMatrix(i, *elmat_p);
299  }
300  }
301  }
302  }
303
304  if (bbfi.Size())
305  {
306  for (i = 0; i < fes -> GetNBE(); i++)
307  {
308  const FiniteElement &be = *fes->GetBE(i);
309  fes -> GetBdrElementVDofs (i, vdofs);
310  eltrans = fes -> GetBdrElementTransformation (i);
311  bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
312  for (int k = 1; k < bbfi.Size(); k++)
313  {
314  bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
315  elmat += elemmat;
316  }
317  if (!static_cond)
318  {
320  }
321  else
322  {
323  static_cond->AssembleBdrMatrix(i, elmat);
324  }
325  }
326  }
327
328  if (fbfi.Size())
329  {
331  Array<int> vdofs2;
332
333  int nfaces = mesh->GetNumFaces();
334  for (i = 0; i < nfaces; i++)
335  {
336  tr = mesh -> GetInteriorFaceTransformations (i);
337  if (tr != NULL)
338  {
339  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
340  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
341  vdofs.Append (vdofs2);
342  for (int k = 0; k < fbfi.Size(); k++)
343  {
344  fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
345  *fes -> GetFE (tr -> Elem2No),
346  *tr, elemmat);
347  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
348  }
349  }
350  }
351  }
352
353  if (bfbfi.Size())
354  {
356  const FiniteElement *fe1, *fe2;
357
358  for (i = 0; i < fes -> GetNBE(); i++)
359  {
360  tr = mesh -> GetBdrFaceTransformations (i);
361  if (tr != NULL)
362  {
363  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
364  fe1 = fes -> GetFE (tr -> Elem1No);
365  // The fe2 object is really a dummy and not used on the boundaries,
366  // but we can't dereference a NULL pointer, and we don't want to
367  // actually make a fake element.
368  fe2 = fe1;
369  for (int k = 0; k < bfbfi.Size(); k++)
370  {
371  bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
372  mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
373  }
374  }
375  }
376  }
377
378 #ifdef MFEM_USE_OPENMP
379  if (free_element_matrices)
380  {
382  }
383 #endif
384 }
385
387 {
388  // Do not remove zero entries to preserve the symmetric structure of the
389  // matrix which in turn will give rise to symmetric structure in the new
390  // matrix. This ensures that subsequent calls to EliminateRowCol will work
391  // correctly.
392  Finalize(0);
393  MFEM_ASSERT(mat, "the BilinearForm is not assembled")
394
396  if (!P) { return; } // conforming mesh
397
398  SparseMatrix *R = Transpose(*P);
399  SparseMatrix *RA = mfem::Mult(*R, *mat);
400  delete mat;
401  if (mat_e)
402  {
403  SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
404  delete mat_e;
405  mat_e = RAe;
406  }
407  delete R;
408  mat = mfem::Mult(*RA, *P);
409  delete RA;
410  if (mat_e)
411  {
412  SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
413  delete mat_e;
414  mat_e = RAeP;
415  }
416
417  height = mat->Height();
418  width = mat->Width();
419 }
420
422  Vector &x, Vector &b,
423  SparseMatrix &A, Vector &X, Vector &B,
424  int copy_interior)
425 {
427  Array<int> ess_rtdof_list;
428
429  // Finish the matrix assembly and perform BC elimination, storing the
430  // eliminated part of the matrix.
431  const int keep_diag = 1;
432  if (static_cond)
433  {
434  static_cond->ConvertListToReducedTrueDofs(ess_tdof_list, ess_rtdof_list);
436  {
437  static_cond->Finalize(); // finalize Schur complement (to true dofs)
438  static_cond->EliminateReducedTrueDofs(ess_rtdof_list, keep_diag);
439  static_cond->Finalize(); // finalize eliminated part
440  }
441  }
442  else if (!mat_e)
443  {
444  if (P) { ConformingAssemble(); }
445  EliminateVDofs(ess_tdof_list, keep_diag);
446  Finalize();
447  }
448
449  // Transform the system and perform the elimination in B, based on the
450  // essential BC values from x. Restrict the BC part of x in X, and set the
451  // non-BC part to zero. Since there is no good initial guess for the Lagrange
452  // multipliers, set X = 0.0 for hybridization.
453  if (static_cond)
454  {
455  // Schur complement reduction to the exposed dofs
456  static_cond->ReduceRHS(b, B);
459  static_cond->GetMatrix().PartMult(ess_rtdof_list, X, B);
460  if (!copy_interior) { X.SetSubVectorComplement(ess_rtdof_list, 0.0); }
462  }
463  else if (!P) // conforming space
464  {
465  if (hybridization)
466  {
467  // Reduction to the Lagrange multipliers system
468  EliminateVDofsInRHS(ess_tdof_list, x, b);
469  hybridization->ReduceRHS(b, B);
470  X.SetSize(B.Size());
471  X = 0.0;
473  }
474  else
475  {
476  // A, X and B point to the same data as mat, x and b
477  EliminateVDofsInRHS(ess_tdof_list, x, b);
478  X.NewDataAndSize(x.GetData(), x.Size());
479  B.NewDataAndSize(b.GetData(), b.Size());
480  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
481  A.MakeRef(*mat);
482  }
483  }
484  else // non-conforming space
485  {
486  if (hybridization)
487  {
488  // Reduction to the Lagrange multipliers system
490  Vector conf_b(P->Width()), conf_x(P->Width());
491  P->MultTranspose(b, conf_b);
492  R->Mult(x, conf_x);
493  EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
494  R->MultTranspose(conf_b, b); // store eliminated rhs in b
495  hybridization->ReduceRHS(conf_b, B);
496  X.SetSize(B.Size());
497  X = 0.0;
499  }
500  else
501  {
502  // Variational restriction with P
504  B.SetSize(P->Width());
505  P->MultTranspose(b, B);
506  X.SetSize(R->Height());
507  R->Mult(x, X);
508  EliminateVDofsInRHS(ess_tdof_list, X, B);
509  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
510  A.MakeRef(*mat);
511  }
512  }
513 }
514
516  const Vector &b, Vector &x)
517 {
519  if (!P) // conforming space
520  {
521  if (static_cond)
522  {
523  // Private dofs back solve
524  static_cond->ComputeSolution(b, X, x);
525  }
526  else if (hybridization)
527  {
528  // Primal unknowns recovery
529  hybridization->ComputeSolution(b, X, x);
530  }
531  else
532  {
533  // X and x point to the same data
534  }
535  }
536  else // non-conforming space
537  {
538  if (static_cond)
539  {
540  // Private dofs back solve
541  static_cond->ComputeSolution(b, X, x);
542  }
543  else if (hybridization)
544  {
545  // Primal unknowns recovery
546  Vector conf_b(P->Width()), conf_x(P->Width());
547  P->MultTranspose(b, conf_b);
549  R->Mult(x, conf_x); // get essential b.c. from x
550  hybridization->ComputeSolution(conf_b, X, conf_x);
551  x.SetSize(P->Height());
552  P->Mult(conf_x, x);
553  }
554  else
555  {
556  // Apply conforming prolongation
557  x.SetSize(P->Height());
558  P->Mult(X, x);
559  }
560  }
561 }
562
564 {
565  if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
566  {
567  return;
568  }
569
570  int num_elements = fes->GetNE();
571  int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
572
573  element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
574  num_elements);
575
576  DenseMatrix tmp;
578
579 #ifdef MFEM_USE_OPENMP
580  #pragma omp parallel for private(tmp,eltrans)
581 #endif
582  for (int i = 0; i < num_elements; i++)
583  {
585  num_dofs_per_el, num_dofs_per_el);
586  const FiniteElement &fe = *fes->GetFE(i);
587 #ifdef MFEM_DEBUG
588  if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
589  mfem_error("BilinearForm::ComputeElementMatrices:"
590  " all elements must have same number of dofs");
591 #endif
592  fes->GetElementTransformation(i, &eltrans);
593
594  dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
595  for (int k = 1; k < dbfi.Size(); k++)
596  {
597  // note: some integrators may not be thread-safe
598  dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
599  elmat += tmp;
600  }
601  elmat.ClearExternalData();
602  }
603 }
604
606  Vector &sol, Vector &rhs, int d)
607 {
608  Array<int> ess_dofs, conf_ess_dofs;
609  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
610
611  if (fes->GetConformingRestriction() == NULL)
612  {
613  EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, d);
614  }
615  else
616  {
617  fes->GetConformingRestriction()->BooleanMult(ess_dofs, conf_ess_dofs);
618  EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, d);
619  }
620 }
621
622 void BilinearForm::EliminateEssentialBC(Array<int> &bdr_attr_is_ess, int d)
623 {
624  Array<int> ess_dofs, conf_ess_dofs;
625  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
626
627  if (fes->GetConformingRestriction() == NULL)
628  {
629  EliminateEssentialBCFromDofs(ess_dofs, d);
630  }
631  else
632  {
633  fes->GetConformingRestriction()->BooleanMult(ess_dofs, conf_ess_dofs);
634  EliminateEssentialBCFromDofs(conf_ess_dofs, d);
635  }
636 }
637
639  double value)
640 {
641  Array<int> ess_dofs, conf_ess_dofs;
642  fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
643
644  if (fes->GetConformingRestriction() == NULL)
645  {
646  EliminateEssentialBCFromDofsDiag(ess_dofs, value);
647  }
648  else
649  {
650  fes->GetConformingRestriction()->BooleanMult(ess_dofs, conf_ess_dofs);
651  EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
652  }
653 }
654
656  Vector &sol, Vector &rhs, int d)
657 {
658  for (int i = 0; i < vdofs.Size(); i++)
659  {
660  int vdof = vdofs[i];
661  if ( vdof >= 0 )
662  {
663  mat -> EliminateRowCol (vdof, sol(vdof), rhs, d);
664  }
665  else
666  {
667  mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, d);
668  }
669  }
670 }
671
673 {
674  if (mat_e == NULL)
675  {
676  mat_e = new SparseMatrix(height);
677  }
678
679  for (int i = 0; i < vdofs.Size(); i++)
680  {
681  int vdof = vdofs[i];
682  if ( vdof >= 0 )
683  {
684  mat -> EliminateRowCol (vdof, *mat_e, d);
685  }
686  else
687  {
688  mat -> EliminateRowCol (-1-vdof, *mat_e, d);
689  }
690  }
691 }
692
694  Array<int> &ess_dofs, Vector &sol, Vector &rhs, int d)
695 {
696  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
697  MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
698  MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
699
700  for (int i = 0; i < ess_dofs.Size(); i++)
701  if (ess_dofs[i] < 0)
702  {
703  mat -> EliminateRowCol (i, sol(i), rhs, d);
704  }
705 }
706
708 {
709  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
710
711  for (int i = 0; i < ess_dofs.Size(); i++)
712  if (ess_dofs[i] < 0)
713  {
714  mat -> EliminateRowCol (i, d);
715  }
716 }
717
719  double value)
720 {
721  MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
722
723  for (int i = 0; i < ess_dofs.Size(); i++)
724  if (ess_dofs[i] < 0)
725  {
726  mat -> EliminateRowColDiag (i, value);
727  }
728 }
729
731  Array<int> &vdofs, const Vector &x, Vector &b)
732 {
734  mat->PartMult(vdofs, x, b);
735 }
736
738 {
739  if (nfes) { fes = nfes; }
740
741  delete mat_e;
742  delete mat;
744  delete static_cond;
745  static_cond = NULL;
746  delete hybridization;
747  hybridization = NULL;
748
749  height = width = fes->GetVSize();
750
751  mat = mat_e = NULL;
752 }
753
755 {
756  delete mat_e;
757  delete mat;
758  delete element_matrices;
759  delete static_cond;
760  delete hybridization;
761
762  if (!extern_bfs)
763  {
764  int k;
765  for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
766  for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
767  for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
768  for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
769  }
770 }
771
772
774  FiniteElementSpace *te_fes)
775  : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
776 {
777  trial_fes = tr_fes;
778  test_fes = te_fes;
779  mat = NULL;
780 }
781
782 double & MixedBilinearForm::Elem (int i, int j)
783 {
784  return (*mat)(i, j);
785 }
786
787 const double & MixedBilinearForm::Elem (int i, int j) const
788 {
789  return (*mat)(i, j);
790 }
791
792 void MixedBilinearForm::Mult (const Vector & x, Vector & y) const
793 {
794  mat -> Mult (x, y);
795 }
796
798  const double a) const
799 {
800  mat -> AddMult (x, y, a);
801 }
802
804  const double a) const
805 {
806  mat -> AddMultTranspose (x, y, a);
807 }
808
810 {
811  return mat -> Inverse ();
812 }
813
814 void MixedBilinearForm::Finalize (int skip_zeros)
815 {
816  mat -> Finalize (skip_zeros);
817 }
818
820 {
821  MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
823  "MixedBilinearForm::GetBlocks: both trial and test spaces "
824  "must use Ordering::byNODES!");
825
826  blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
827
828  mat->GetBlocks(blocks);
829 }
830
832 {
833  dom.Append (bfi);
834 }
835
837 {
838  bdr.Append (bfi);
839 }
840
842 {
843  skt.Append (bfi);
844 }
845
846 void MixedBilinearForm::Assemble (int skip_zeros)
847 {
848  int i, k;
849  Array<int> tr_vdofs, te_vdofs;
850  ElementTransformation *eltrans;
851  DenseMatrix elemmat;
852
853  Mesh *mesh = test_fes -> GetMesh();
854
855  if (mat == NULL)
856  {
857  mat = new SparseMatrix(height, width);
858  }
859
860  if (dom.Size())
861  {
862  for (i = 0; i < test_fes -> GetNE(); i++)
863  {
864  trial_fes -> GetElementVDofs (i, tr_vdofs);
865  test_fes -> GetElementVDofs (i, te_vdofs);
866  eltrans = test_fes -> GetElementTransformation (i);
867  for (k = 0; k < dom.Size(); k++)
868  {
869  dom[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
870  *test_fes -> GetFE(i),
871  *eltrans, elemmat);
872  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
873  }
874  }
875  }
876
877  if (bdr.Size())
878  {
879  for (i = 0; i < test_fes -> GetNBE(); i++)
880  {
881  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
882  test_fes -> GetBdrElementVDofs (i, te_vdofs);
883  eltrans = test_fes -> GetBdrElementTransformation (i);
884  for (k = 0; k < bdr.Size(); k++)
885  {
886  bdr[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
887  *test_fes -> GetBE(i),
888  *eltrans, elemmat);
889  mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
890  }
891  }
892  }
893
894  if (skt.Size())
895  {
897  Array<int> te_vdofs2;
898  const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
899
900  int nfaces = mesh->GetNumFaces();
901  for (i = 0; i < nfaces; i++)
902  {
903  ftr = mesh->GetFaceElementTransformations(i);
904  trial_fes->GetFaceVDofs(i, tr_vdofs);
905  test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
906  trial_face_fe = trial_fes->GetFaceElement(i);
907  test_fe1 = test_fes->GetFE(ftr->Elem1No);
908  if (ftr->Elem2No >= 0)
909  {
910  test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
911  te_vdofs.Append(te_vdofs2);
912  test_fe2 = test_fes->GetFE(ftr->Elem2No);
913  }
914  else
915  {
916  // The test_fe2 object is really a dummy and not used on the
917  // boundaries, but we can't dereference a NULL pointer, and we don't
918  // want to actually make a fake element.
919  test_fe2 = test_fe1;
920  }
921  for (int k = 0; k < skt.Size(); k++)
922  {
923  skt[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
924  *ftr, elemmat);
926  }
927  }
928  }
929 }
930
932 {
933  Finalize();
934
936  if (P2)
937  {
938  SparseMatrix *R = Transpose(*P2);
939  SparseMatrix *RA = mfem::Mult(*R, *mat);
940  delete R;
941  delete mat;
942  mat = RA;
943  }
944
946  if (P1)
947  {
948  SparseMatrix *RAP = mfem::Mult(*mat, *P1);
949  delete mat;
950  mat = RAP;
951  }
952
953  height = mat->Height();
954  width = mat->Width();
955 }
956
958  Array<int> &bdr_attr_is_ess, Vector &sol, Vector &rhs )
959 {
960  int i, j, k;
961  Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
962
963  cols_marker = 0;
964  for (i = 0; i < trial_fes -> GetNBE(); i++)
965  if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
966  {
967  trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
968  for (j = 0; j < tr_vdofs.Size(); j++)
969  {
970  if ( (k = tr_vdofs[j]) < 0 )
971  {
972  k = -1-k;
973  }
974  cols_marker[k] = 1;
975  }
976  }
977  mat -> EliminateCols (cols_marker, &sol, &rhs);
978 }
979
981  Array<int> &marked_vdofs, Vector &sol, Vector &rhs)
982 {
983  mat -> EliminateCols (marked_vdofs, &sol, &rhs);
984 }
985
987 {
988  int i, j, k;
989  Array<int> te_vdofs;
990
991  for (i = 0; i < test_fes -> GetNBE(); i++)
992  if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
993  {
994  test_fes -> GetBdrElementVDofs (i, te_vdofs);
995  for (j = 0; j < te_vdofs.Size(); j++)
996  {
997  if ( (k = te_vdofs[j]) < 0 )
998  {
999  k = -1-k;
1000  }
1001  mat -> EliminateRow (k);
1002  }
1003  }
1004 }
1005
1007 {
1008  delete mat;
1009  mat = NULL;
1010  height = test_fes->GetVSize();
1011  width = trial_fes->GetVSize();
1012 }
1013
1015 {
1016  int i;
1017
1018  if (mat) { delete mat; }
1019  for (i = 0; i < dom.Size(); i++) { delete dom[i]; }
1020  for (i = 0; i < bdr.Size(); i++) { delete bdr[i]; }
1021  for (i = 0; i < skt.Size(); i++) { delete skt[i]; }
1022 }
1023
1024
1026 {
1027  Array<int> dom_vdofs, ran_vdofs;
1029  const FiniteElement *dom_fe, *ran_fe;
1030  DenseMatrix totelmat, elmat;
1031
1032  if (mat == NULL)
1033  {
1034  mat = new SparseMatrix(height, width);
1035  }
1036
1037  if (dom.Size() > 0)
1038  {
1039  for (int i = 0; i < test_fes->GetNE(); i++)
1040  {
1041  trial_fes->GetElementVDofs(i, dom_vdofs);
1042  test_fes->GetElementVDofs(i, ran_vdofs);
1044  dom_fe = trial_fes->GetFE(i);
1045  ran_fe = test_fes->GetFE(i);
1046
1047  dom[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1048  for (int j = 1; j < dom.Size(); j++)
1049  {
1050  dom[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1051  totelmat += elmat;
1052  }
1053  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1054  }
1055  }
1056
1057  if (skt.Size())
1058  {
1059  const int nfaces = test_fes->GetMesh()->GetNumFaces();
1060  for (int i = 0; i < nfaces; i++)
1061  {
1062  trial_fes->GetFaceVDofs(i, dom_vdofs);
1063  test_fes->GetFaceVDofs(i, ran_vdofs);
1065  dom_fe = trial_fes->GetFaceElement(i);
1066  ran_fe = test_fes->GetFaceElement(i);
1067
1068  skt[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
1069  for (int j = 1; j < skt.Size(); j++)
1070  {
1071  skt[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
1072  totelmat += elmat;
1073  }
1074  mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1075  }
1076  }
1077 }
1078
1079 }
Abstract class for Finite Elements.
Definition: fe.hpp:44
void EliminateEssentialBCFromDofs(Array< int > &ess_dofs, Vector &sol, Vector &rhs, int d=0)
int Size() const
Logical size of the array.
Definition: array.hpp:109
int GetVSize() const
Definition: fespace.hpp:164
Array< BilinearFormIntegrator * > * GetBBFI()
int * GetJ()
Definition: table.hpp:108
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
void NewDataAndSize(double *d, int s)
Definition: vector.hpp:72
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Array< BilinearFormIntegrator * > skt
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
void ReduceRHS(const Vector &b, Vector &sc_b) const
Definition: staticcond.cpp:309
Array< BilinearFormIntegrator * > * GetDBFI()
void MakeRef(const SparseMatrix &master)
Definition: sparsemat.cpp:164
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
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:161
Array< BilinearFormIntegrator * > dbfi
Set of Domain Integrators to be applied.
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs) const
Definition: fespace.cpp:432
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:573
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:445
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
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:423
int Size() const
Returns the size of the vector.
Definition: vector.hpp:84
Abstract data type for matrix inverse.
Definition: matrix.hpp:58
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:138
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1445
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:452
Array< BilinearFormIntegrator * > * GetFBFI()
double * GetData() const
Definition: vector.hpp:88
MixedBilinearForm(FiniteElementSpace *tr_fes, FiniteElementSpace *te_fes)
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:557
void FormLinearSystem(Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Definition: sparsemat.cpp:752
void EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3755
void EliminateEssentialBCFromTrialDofs(Array< int > &marked_vdofs, Vector &sol, Vector &rhs)
Array< BilinearFormIntegrator * > bdr
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:185
void EliminateVDofs(Array< int > &vdofs, Vector &sol, Vector &rhs, int d=0)
Eliminate the given vdofs. NOTE: here, vdofs is a list of DOFs.
void SetSize(int m, int n)
Definition: array.hpp:256
void Finalize()
Finalize the construction of the hybridized matrix.
const SparseMatrix * GetConformingRestriction()
Definition: fespace.cpp:919
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:35
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:368
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:140
SparseMatrix * mat
Sparse matrix to be associated with the form.
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
Definition: mesh.cpp:280
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
Abstract data type matrix.
Definition: matrix.hpp:27
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1318
void Assemble(int skip_zeros=1)
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:385
double * GetData(int k)
Definition: densemat.hpp:600
virtual void Update(FiniteElementSpace *nfes=NULL)
void EliminateTrialDofs(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs)
const SparseMatrix * GetConformingProlongation()
Definition: fespace.cpp:912
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:154
virtual void Finalize(int skip_zeros=1)
Definition: sparsemat.cpp:680
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1733
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
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:575
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1809
int SizeI() const
Definition: densemat.hpp:579
void ComputeElementMatrix(int i, DenseMatrix &elmat)
void EliminateEssentialBCDiag(Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
int GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:176
SparseMatrix * mat_e
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i&#39;th element.
Definition: fespace.hpp:206
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:22
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.
Abstract finite element space.
Definition: fespace.hpp:62
int GetDof() const
Returns the degrees of freedom in the FE space.
Definition: fe.hpp:88
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:417
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:23
void ClearExternalData()
Definition: densemat.hpp:65
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:286
Table * GetFaceToElementTable() const
Definition: mesh.cpp:4247
FiniteElementSpace * test_fes
int SizeJ() const
Definition: densemat.hpp:580
virtual void EliminateTestDofs(Array< int > &bdr_attr_is_ess)
void EliminateVDofsInRHS(Array< int > &vdofs, const Vector &x, Vector &b)
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AssembleMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:189
bool HasEliminatedBC() const
Definition: staticcond.hpp:131
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:492
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:1199
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:173
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:33
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
Definition: staticcond.hpp:141
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:107
const Table & GetElementToDofTable() const
Definition: fespace.hpp:278
void ReduceSolution(const Vector &sol, Vector &sc_sol) const
Definition: staticcond.cpp:388
SparseMatrix & GetMatrixElim()
Return the eliminated part of the serial Schur complement matrix.
Definition: staticcond.hpp:144
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 for the i&#39;th boundary element.
Definition: fespace.cpp:1419
Array< int > vdofs
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:71
void FreeElementMatrices()
Free the memory used by the element matrices.
void EliminateEssentialBCFromDofsDiag(Array< int > &ess_dofs, double value)
Perform elimination and set the diagonal entry to the given value.
void Init(bool symmetric, bool block_diagonal)
Definition: staticcond.cpp:134
virtual void Assemble(int skip_zeros=1)
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:565
void ConvertListToReducedTrueDofs(const Array< int > &ess_tdof_list, Array< int > &ess_rtdof_list) const
Definition: staticcond.hpp:169
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
void EnableStaticCondensation()