MFEM v4.7.0 Finite element discretization library
Loading...
Searching...
No Matches
bilinearform.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// Implementation of class BilinearForm
13
14#include "fem.hpp"
15#include "../general/device.hpp"
16#include "../mesh/nurbs.hpp"
17#include <cmath>
18
19namespace mfem
20{
21
23{
24 if (static_cond) { return; }
25
26 if (precompute_sparsity == 0 || fes->GetVDim() > 1)
27 {
28 mat = new SparseMatrix(height);
29 return;
30 }
31
32 const Table &elem_dof = fes->GetElementToDofTable();
33 Table dof_dof;
34
35 if (interior_face_integs.Size() > 0)
36 {
37 // the sparsity pattern is defined from the map: face->element->dof
38 Table face_dof, dof_face;
39 {
40 Table *face_elem = fes->GetMesh()->GetFaceToElementTable();
41 mfem::Mult(*face_elem, elem_dof, face_dof);
42 delete face_elem;
43 }
44 Transpose(face_dof, dof_face, height);
45 mfem::Mult(dof_face, face_dof, dof_dof);
46 }
47 else
48 {
49 // the sparsity pattern is defined from the map: element->dof
50 Table dof_elem;
51 Transpose(elem_dof, dof_elem, height);
52 mfem::Mult(dof_elem, elem_dof, dof_dof);
53 }
54
55 dof_dof.SortRows();
56
57 int *I = dof_dof.GetI();
58 int *J = dof_dof.GetJ();
59 real_t *data = Memory<real_t>(I[height]);
60
61 mat = new SparseMatrix(I, J, data, height, height, true, true, true);
62 *mat = 0.0;
63
64 dof_dof.LoseData();
65}
66
68 : Matrix (f->GetVSize())
69{
70 fes = f;
72 mat = mat_e = NULL;
73 extern_bfs = 0;
74 element_matrices = NULL;
75 static_cond = NULL;
76 hybridization = NULL;
79
81 batch = 1;
82 ext = NULL;
83}
84
86 : Matrix (f->GetVSize())
87{
88 fes = f;
90 mat_e = NULL;
91 extern_bfs = 1;
92 element_matrices = NULL;
93 static_cond = NULL;
94 hybridization = NULL;
97
99 batch = 1;
100 ext = NULL;
101
102 // Copy the pointers to the integrators
105
108
110
113
114 AllocMat();
115}
116
118{
119 if (ext)
120 {
121 MFEM_ABORT("the assembly level has already been set!");
122 }
123 assembly = assembly_level;
124 switch (assembly)
125 {
127 break;
129 SetDiagonalPolicy( DIAG_ONE ); // Only diagonal policy supported on device
130 ext = new FABilinearFormExtension(this);
131 break;
133 ext = new EABilinearFormExtension(this);
134 break;
136 ext = new PABilinearFormExtension(this);
137 break;
139 ext = new MFBilinearFormExtension(this);
140 break;
141 default:
142 MFEM_ABORT("BilinearForm: unknown assembly level");
143 }
144}
145
147{
148 delete static_cond;
150 {
151 static_cond = NULL;
152 MFEM_WARNING("Static condensation not supported for this assembly level");
153 return;
154 }
157 {
158 bool symmetric = false; // TODO
159 bool block_diagonal = false; // TODO
160 static_cond->Init(symmetric, block_diagonal);
161 }
162 else
163 {
164 delete static_cond;
165 static_cond = NULL;
166 }
167}
168
170 BilinearFormIntegrator *constr_integ,
171 const Array<int> &ess_tdof_list)
172{
173 delete hybridization;
175 {
176 delete constr_integ;
177 hybridization = NULL;
178 MFEM_WARNING("Hybridization not supported for this assembly level");
179 return;
180 }
181 hybridization = new Hybridization(fes, constr_space);
183 hybridization->Init(ess_tdof_list);
184}
185
186void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
187{
188 if (static_cond) { return; }
189
190 if (mat)
191 {
192 if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
193 {
194 return; // mat is already using the given sparsity
195 }
196 delete mat;
197 }
198 height = width = fes->GetVSize();
199 mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
200}
201
203{
204 MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
205 "invalid matrix A dimensions: "
206 << A.Height() << " x " << A.Width());
207 MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
208
209 UseSparsity(A.GetI(), A.GetJ(), A.ColumnsAreSorted());
210}
211
213{
214 return mat -> Elem(i,j);
215}
216
217const real_t& BilinearForm::Elem (int i, int j) const
218{
219 return mat -> Elem(i,j);
220}
221
223{
224 return mat -> Inverse();
225}
226
227void BilinearForm::Finalize (int skip_zeros)
228{
230 {
231 if (!static_cond) { mat->Finalize(skip_zeros); }
232 if (mat_e) { mat_e->Finalize(skip_zeros); }
233 if (static_cond) { static_cond->Finalize(); }
235 }
236}
237
239{
240 domain_integs.Append(bfi);
241 domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
242}
243
245 Array<int> &elem_marker)
246{
247 domain_integs.Append(bfi);
248 domain_integs_marker.Append(&elem_marker);
249}
250
252{
253 boundary_integs.Append (bfi);
254 boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
255}
256
258 Array<int> &bdr_marker)
259{
260 boundary_integs.Append (bfi);
261 boundary_integs_marker.Append(&bdr_marker);
262}
263
265{
266 interior_face_integs.Append (bfi);
267}
268
270{
271 boundary_face_integs.Append(bfi);
272 // NULL marker means apply everywhere
273 boundary_face_integs_marker.Append(NULL);
274}
275
277 Array<int> &bdr_marker)
278{
279 boundary_face_integs.Append(bfi);
280 boundary_face_integs_marker.Append(&bdr_marker);
281}
282
284{
286 {
288 elmat = element_matrices->GetData(i);
289 return;
290 }
291
292 if (domain_integs.Size())
293 {
294 const FiniteElement &fe = *fes->GetFE(i);
296 domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
297 for (int k = 1; k < domain_integs.Size(); k++)
298 {
299 domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
300 elmat += elemmat;
301 }
302 }
303 else
304 {
306 elmat.SetSize(vdofs.Size());
307 elmat = 0.0;
308 }
309}
310
312{
313 if (boundary_integs.Size())
314 {
315 const FiniteElement &be = *fes->GetBE(i);
317 boundary_integs[0]->AssembleElementMatrix(be, *eltrans, elmat);
318 for (int k = 1; k < boundary_integs.Size(); k++)
319 {
320 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
321 elmat += elemmat;
322 }
323 }
324 else
325 {
327 elmat.SetSize(vdofs.Size());
328 elmat = 0.0;
329 }
330}
331
333 int i, const DenseMatrix &elmat, int skip_zeros)
334{
335 AssembleElementMatrix(i, elmat, vdofs, skip_zeros);
336}
337
339 int i, const DenseMatrix &elmat, Array<int> &vdofs_, int skip_zeros)
340{
341 fes->GetElementVDofs(i, vdofs_);
342 if (static_cond)
343 {
344 static_cond->AssembleMatrix(i, elmat);
345 }
346 else
347 {
348 if (mat == NULL)
349 {
350 AllocMat();
351 }
352 mat->AddSubMatrix(vdofs_, vdofs_, elmat, skip_zeros);
353 if (hybridization)
354 {
355 hybridization->AssembleMatrix(i, elmat);
356 }
357 }
358}
359
361 int i, const DenseMatrix &elmat, int skip_zeros)
362{
363 AssembleBdrElementMatrix(i, elmat, vdofs, skip_zeros);
364}
365
367 int i, const DenseMatrix &elmat, Array<int> &vdofs_, int skip_zeros)
368{
369 fes->GetBdrElementVDofs(i, vdofs_);
370 if (static_cond)
371 {
373 }
374 else
375 {
376 if (mat == NULL)
377 {
378 AllocMat();
379 }
380 mat->AddSubMatrix(vdofs_, vdofs_, elmat, skip_zeros);
381 if (hybridization)
382 {
384 }
385 }
386}
387
388void BilinearForm::Assemble(int skip_zeros)
389{
390 if (ext)
391 {
392 ext->Assemble();
393 return;
394 }
395
396 ElementTransformation *eltrans;
397 DofTransformation * doftrans;
398 Mesh *mesh = fes -> GetMesh();
399 DenseMatrix elmat, *elmat_p;
400
401 if (mat == NULL)
402 {
403 AllocMat();
404 }
405
406#ifdef MFEM_USE_LEGACY_OPENMP
407 int free_element_matrices = 0;
408 if (!element_matrices)
409 {
411 free_element_matrices = 1;
412 }
413#endif
414
415 if (domain_integs.Size())
416 {
417 for (int k = 0; k < domain_integs.Size(); k++)
418 {
419 if (domain_integs_marker[k] != NULL)
420 {
421 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
422 (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
423 "invalid element marker for domain integrator #"
424 << k << ", counting from zero");
425 }
426
427 if (domain_integs[k]->Patchwise())
428 {
429 MFEM_VERIFY(fes->GetNURBSext(), "Patchwise integration requires a "
430 << "NURBS FE space");
431 }
432 }
433
434 // Element-wise integration
435 for (int i = 0; i < fes -> GetNE(); i++)
436 {
437 // Set both doftrans (potentially needed to assemble the element
438 // matrix) and vdofs, which is also needed when the element matrices
439 // are pre-assembled.
440 doftrans = fes->GetElementVDofs(i, vdofs);
442 {
443 elmat_p = &(*element_matrices)(i);
444 }
445 else
446 {
447 const int elem_attr = fes->GetMesh()->GetAttribute(i);
448 eltrans = fes->GetElementTransformation(i);
449
450 elmat.SetSize(0);
451 for (int k = 0; k < domain_integs.Size(); k++)
452 {
453 if ((domain_integs_marker[k] == NULL ||
454 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
455 && !domain_integs[k]->Patchwise())
456 {
457 domain_integs[k]->AssembleElementMatrix(*fes->GetFE(i),
458 *eltrans, elemmat);
459 if (elmat.Size() == 0)
460 {
461 elmat = elemmat;
462 }
463 else
464 {
465 elmat += elemmat;
466 }
467 }
468 }
469 if (elmat.Size() == 0)
470 {
471 continue;
472 }
473 else
474 {
475 elmat_p = &elmat;
476 }
477 if (doftrans)
478 {
479 doftrans->TransformDual(elmat);
480 }
481 elmat_p = &elmat;
482 }
483 if (static_cond)
484 {
485 static_cond->AssembleMatrix(i, *elmat_p);
486 }
487 else
488 {
489 mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
490 if (hybridization)
491 {
492 hybridization->AssembleMatrix(i, *elmat_p);
493 }
494 }
495 }
496
497 // Patch-wise integration
498 if (fes->GetNURBSext())
499 {
500 for (int p=0; p<mesh->NURBSext->GetNP(); ++p)
501 {
502 bool vdofsSet = false;
503 for (int k = 0; k < domain_integs.Size(); k++)
504 {
505 if (domain_integs[k]->Patchwise())
506 {
507 if (!vdofsSet)
508 {
510 vdofsSet = true;
511 }
512
513 SparseMatrix* spmat = nullptr;
514 domain_integs[k]->AssemblePatchMatrix(p, *fes, spmat);
515 Array<int> cols;
516 Vector srow;
517
518 for (int r=0; r<spmat->Height(); ++r)
519 {
520 spmat->GetRow(r, cols, srow);
521 for (int i=0; i<cols.Size(); ++i)
522 {
523 cols[i] = vdofs[cols[i]];
524 }
525 mat->AddRow(vdofs[r], cols, srow);
526 }
527
528 delete spmat;
529 }
530 }
531 }
532 }
533 }
534
535 if (boundary_integs.Size())
536 {
537 // Which boundary attributes need to be processed?
538 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
539 mesh->bdr_attributes.Max() : 0);
540 bdr_attr_marker = 0;
541 for (int k = 0; k < boundary_integs.Size(); k++)
542 {
543 if (boundary_integs_marker[k] == NULL)
544 {
545 bdr_attr_marker = 1;
546 break;
547 }
548 Array<int> &bdr_marker = *boundary_integs_marker[k];
549 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
550 "invalid boundary marker for boundary integrator #"
551 << k << ", counting from zero");
552 for (int i = 0; i < bdr_attr_marker.Size(); i++)
553 {
554 bdr_attr_marker[i] |= bdr_marker[i];
555 }
556 }
557
558 for (int i = 0; i < fes -> GetNBE(); i++)
559 {
560 const int bdr_attr = mesh->GetBdrAttribute(i);
561 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
562
563 const FiniteElement &be = *fes->GetBE(i);
564 doftrans = fes -> GetBdrElementVDofs (i, vdofs);
565 eltrans = fes -> GetBdrElementTransformation (i);
566 int k = 0;
567 for (; k < boundary_integs.Size(); k++)
568 {
569 if (boundary_integs_marker[k] &&
570 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
571
572 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elmat);
573 k++;
574 break;
575 }
576 for (; k < boundary_integs.Size(); k++)
577 {
578 if (boundary_integs_marker[k] &&
579 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
580
581 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
582 elmat += elemmat;
583 }
584 if (doftrans)
585 {
586 doftrans->TransformDual(elmat);
587 }
588 elmat_p = &elmat;
589 if (!static_cond)
590 {
591 mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
592 if (hybridization)
593 {
594 hybridization->AssembleBdrMatrix(i, *elmat_p);
595 }
596 }
597 else
598 {
599 static_cond->AssembleBdrMatrix(i, *elmat_p);
600 }
601 }
602 }
603
604 if (interior_face_integs.Size())
605 {
607 Array<int> vdofs2;
608
609 int nfaces = mesh->GetNumFaces();
610 for (int i = 0; i < nfaces; i++)
611 {
612 tr = mesh -> GetInteriorFaceTransformations (i);
613 if (tr != NULL)
614 {
615 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
616 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
617 vdofs.Append (vdofs2);
618 for (int k = 0; k < interior_face_integs.Size(); k++)
619 {
621 AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
622 *fes->GetFE(tr->Elem2No),
623 *tr, elemmat);
624 mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
625 }
626 }
627 }
628 }
629
630 if (boundary_face_integs.Size())
631 {
633 const FiniteElement *fe1, *fe2;
634
635 // Which boundary attributes need to be processed?
636 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
637 mesh->bdr_attributes.Max() : 0);
638 bdr_attr_marker = 0;
639 for (int k = 0; k < boundary_face_integs.Size(); k++)
640 {
641 if (boundary_face_integs_marker[k] == NULL)
642 {
643 bdr_attr_marker = 1;
644 break;
645 }
646 Array<int> &bdr_marker = *boundary_face_integs_marker[k];
647 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
648 "invalid boundary marker for boundary face integrator #"
649 << k << ", counting from zero");
650 for (int i = 0; i < bdr_attr_marker.Size(); i++)
651 {
652 bdr_attr_marker[i] |= bdr_marker[i];
653 }
654 }
655
656 for (int i = 0; i < fes -> GetNBE(); i++)
657 {
658 const int bdr_attr = mesh->GetBdrAttribute(i);
659 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
660
661 tr = mesh -> GetBdrFaceTransformations (i);
662 if (tr != NULL)
663 {
664 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
665 fe1 = fes -> GetFE (tr -> Elem1No);
666 // The fe2 object is really a dummy and not used on the boundaries,
667 // but we can't dereference a NULL pointer, and we don't want to
668 // actually make a fake element.
669 fe2 = fe1;
670 for (int k = 0; k < boundary_face_integs.Size(); k++)
671 {
673 (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
674 { continue; }
675
676 boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
677 elemmat);
678 mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
679 }
680 }
681 }
682 }
683
684#ifdef MFEM_USE_LEGACY_OPENMP
685 if (free_element_matrices)
686 {
688 }
689#endif
690}
691
693{
694 // Do not remove zero entries to preserve the symmetric structure of the
695 // matrix which in turn will give rise to symmetric structure in the new
696 // matrix. This ensures that subsequent calls to EliminateRowCol will work
697 // correctly.
698 Finalize(0);
699 MFEM_ASSERT(mat, "the BilinearForm is not assembled");
700
702 if (!P) { return; } // conforming mesh
703
704 SparseMatrix *R = Transpose(*P);
705 SparseMatrix *RA = mfem::Mult(*R, *mat);
706 delete mat;
707 if (mat_e)
708 {
709 SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
710 delete mat_e;
711 mat_e = RAe;
712 }
713 delete R;
714 mat = mfem::Mult(*RA, *P);
715 delete RA;
716 if (mat_e)
717 {
718 SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
719 delete mat_e;
720 mat_e = RAeP;
721 }
722
723 height = mat->Height();
724 width = mat->Width();
725}
726
728{
729 MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
730 "Vector for holding diagonal has wrong size!");
732 if (!ext)
733 {
734 MFEM_ASSERT(mat, "the BilinearForm is not assembled!");
735 MFEM_ASSERT(cP == nullptr || mat->Height() == cP->Width(),
736 "BilinearForm::ConformingAssemble() is not called!");
737 mat->GetDiag(diag);
738 return;
739 }
740 // Here, we have extension, ext.
741 if (!cP)
742 {
743 ext->AssembleDiagonal(diag);
744 return;
745 }
746 // Here, we have extension, ext, and conforming prolongation, cP.
747
748 // For an AMR mesh, a convergent diagonal is assembled with |P^T| d_l,
749 // where |P^T| has the entry-wise absolute values of the conforming
750 // prolongation transpose operator.
751 Vector local_diag(cP->Height());
752 ext->AssembleDiagonal(local_diag);
753 cP->AbsMultTranspose(local_diag, diag);
754}
755
757 Vector &b, OperatorHandle &A, Vector &X,
758 Vector &B, int copy_interior)
759{
760 if (ext)
761 {
762 ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
763 return;
764 }
766 FormSystemMatrix(ess_tdof_list, A);
767
768 // Transform the system and perform the elimination in B, based on the
769 // essential BC values from x. Restrict the BC part of x in X, and set the
770 // non-BC part to zero. Since there is no good initial guess for the Lagrange
771 // multipliers, set X = 0.0 for hybridization.
772 if (static_cond)
773 {
774 // Schur complement reduction to the exposed dofs
775 static_cond->ReduceSystem(x, b, X, B, copy_interior);
776 }
777 else if (!P) // conforming space
778 {
779 if (hybridization)
780 {
781 // Reduction to the Lagrange multipliers system
782 EliminateVDofsInRHS(ess_tdof_list, x, b);
784 X.SetSize(B.Size());
785 X = 0.0;
786 }
787 else
788 {
789 // A, X and B point to the same data as mat, x and b
790 EliminateVDofsInRHS(ess_tdof_list, x, b);
791 X.MakeRef(x, 0, x.Size());
792 B.MakeRef(b, 0, b.Size());
793 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
794 }
795 }
796 else // non-conforming space
797 {
798 if (hybridization)
799 {
800 // Reduction to the Lagrange multipliers system
802 Vector conf_b(P->Width()), conf_x(P->Width());
803 P->MultTranspose(b, conf_b);
804 R->Mult(x, conf_x);
805 EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
806 R->MultTranspose(conf_b, b); // store eliminated rhs in b
807 hybridization->ReduceRHS(conf_b, B);
808 X.SetSize(B.Size());
809 X = 0.0;
810 }
811 else
812 {
813 // Variational restriction with P
815 B.SetSize(P->Width());
816 P->MultTranspose(b, B);
817 X.SetSize(R->Height());
818 R->Mult(x, X);
819 EliminateVDofsInRHS(ess_tdof_list, X, B);
820 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
821 }
822 }
823}
824
827{
828 if (ext)
829 {
830 ext->FormSystemMatrix(ess_tdof_list, A);
831 return;
832 }
833
834 // Finish the matrix assembly and perform BC elimination, storing the
835 // eliminated part of the matrix.
836 if (static_cond)
837 {
839 {
840 static_cond->SetEssentialTrueDofs(ess_tdof_list);
841 static_cond->Finalize(); // finalize Schur complement (to true dofs)
843 static_cond->Finalize(); // finalize eliminated part
844 }
845 A.Reset(&static_cond->GetMatrix(), false);
846 }
847 else
848 {
849 if (!mat_e)
850 {
852 if (P) { ConformingAssemble(); }
853 EliminateVDofs(ess_tdof_list, diag_policy);
854 const int remove_zeros = 0;
855 Finalize(remove_zeros);
856 }
857 if (hybridization)
858 {
859 A.Reset(&hybridization->GetMatrix(), false);
860 }
861 else
862 {
863 A.Reset(mat, false);
864 }
865 }
866}
867
869 const Vector &b, Vector &x)
870{
871 if (ext)
872 {
873 ext->RecoverFEMSolution(X, b, x);
874 return;
875 }
876
878 if (!P) // conforming space
879 {
880 if (static_cond)
881 {
882 // Private dofs back solve
884 }
885 else if (hybridization)
886 {
887 // Primal unknowns recovery
889 }
890 else
891 {
892 // X and x point to the same data
893
894 // If the validity flags of X's Memory were changed (e.g. if it was
895 // moved to device memory) then we need to tell x about that.
896 x.SyncMemory(X);
897 }
898 }
899 else // non-conforming space
900 {
901 if (static_cond)
902 {
903 // Private dofs back solve
905 }
906 else if (hybridization)
907 {
908 // Primal unknowns recovery
909 Vector conf_b(P->Width()), conf_x(P->Width());
910 P->MultTranspose(b, conf_b);
912 R->Mult(x, conf_x); // get essential b.c. from x
913 hybridization->ComputeSolution(conf_b, X, conf_x);
914 x.SetSize(P->Height());
915 P->Mult(conf_x, x);
916 }
917 else
918 {
919 // Apply conforming prolongation
920 x.SetSize(P->Height());
921 P->Mult(X, x);
922 }
923 }
924}
925
927{
928 if (element_matrices || domain_integs.Size() == 0 || fes->GetNE() == 0)
929 {
930 return;
931 }
932
933 int num_elements = fes->GetNE();
934 int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
935
936 element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
937 num_elements);
938
939 DenseMatrix tmp;
941
942#ifdef MFEM_USE_LEGACY_OPENMP
943 #pragma omp parallel for private(tmp,eltrans)
944#endif
945 for (int i = 0; i < num_elements; i++)
946 {
948 num_dofs_per_el, num_dofs_per_el);
949 const FiniteElement &fe = *fes->GetFE(i);
950#ifdef MFEM_DEBUG
951 if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
952 mfem_error("BilinearForm::ComputeElementMatrices:"
953 " all elements must have same number of dofs");
954#endif
955 fes->GetElementTransformation(i, &eltrans);
956
957 domain_integs[0]->AssembleElementMatrix(fe, eltrans, elmat);
958 for (int k = 1; k < domain_integs.Size(); k++)
959 {
960 // note: some integrators may not be thread-safe
961 domain_integs[k]->AssembleElementMatrix(fe, eltrans, tmp);
962 elmat += tmp;
963 }
964 elmat.ClearExternalData();
965 }
966}
967
969 const Vector &sol, Vector &rhs,
970 DiagonalPolicy dpolicy)
971{
972 Array<int> ess_dofs, conf_ess_dofs;
973 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
974
975 if (fes->GetVSize() == height)
976 {
977 EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
978 }
979 else
980 {
981 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
982 EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
983 }
984}
985
987 DiagonalPolicy dpolicy)
988{
989 Array<int> ess_dofs, conf_ess_dofs;
990 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
991
992 if (fes->GetVSize() == height)
993 {
994 EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
995 }
996 else
997 {
998 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
999 EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
1000 }
1001}
1002
1004 real_t value)
1005{
1006 Array<int> ess_dofs, conf_ess_dofs;
1007 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
1008
1009 if (fes->GetVSize() == height)
1010 {
1011 EliminateEssentialBCFromDofsDiag(ess_dofs, value);
1012 }
1013 else
1014 {
1015 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
1016 EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
1017 }
1018}
1019
1021 const Vector &sol, Vector &rhs,
1022 DiagonalPolicy dpolicy)
1023{
1024 vdofs_.HostRead();
1025 for (int i = 0; i < vdofs_.Size(); i++)
1026 {
1027 int vdof = vdofs_[i];
1028 if ( vdof >= 0 )
1029 {
1030 mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
1031 }
1032 else
1033 {
1034 mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
1035 }
1036 }
1037}
1038
1040 DiagonalPolicy dpolicy)
1041{
1042 if (mat_e == NULL)
1043 {
1044 mat_e = new SparseMatrix(height);
1045 }
1046
1047 vdofs_.HostRead();
1048 for (int i = 0; i < vdofs_.Size(); i++)
1049 {
1050 int vdof = vdofs_[i];
1051 if ( vdof >= 0 )
1052 {
1053 mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
1054 }
1055 else
1056 {
1057 mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
1058 }
1059 }
1060}
1061
1063 const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
1064 DiagonalPolicy dpolicy)
1065{
1066 MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
1067 MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
1068 MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
1069
1070 for (int i = 0; i < ess_dofs.Size(); i++)
1071 if (ess_dofs[i] < 0)
1072 {
1073 mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
1074 }
1075}
1076
1078 DiagonalPolicy dpolicy)
1079{
1080 MFEM_ASSERT(ess_dofs.Size() == height,
1081 "incorrect dof Array size: " << ess_dofs.Size() << ' ' << height);
1082
1083 for (int i = 0; i < ess_dofs.Size(); i++)
1084 if (ess_dofs[i] < 0)
1085 {
1086 mat -> EliminateRowCol (i, dpolicy);
1087 }
1088}
1089
1091 real_t value)
1092{
1093 MFEM_ASSERT(ess_dofs.Size() == height,
1094 "incorrect dof Array size: " << ess_dofs.Size() << ' ' << height);
1095
1096 for (int i = 0; i < ess_dofs.Size(); i++)
1097 if (ess_dofs[i] < 0)
1098 {
1099 mat -> EliminateRowColDiag (i, value);
1100 }
1101}
1102
1104 const Array<int> &vdofs_, const Vector &x, Vector &b)
1105{
1106 mat_e->AddMult(x, b, -1.);
1107 mat->PartMult(vdofs_, x, b);
1108}
1109
1110void BilinearForm::Mult(const Vector &x, Vector &y) const
1111{
1112 if (ext)
1113 {
1114 ext->Mult(x, y);
1115 }
1116 else
1117 {
1118 mat->Mult(x, y);
1119 }
1120}
1121
1123{
1124 if (ext)
1125 {
1126 ext->MultTranspose(x, y);
1127 }
1128 else
1129 {
1130 y = 0.0;
1131 AddMultTranspose (x, y);
1132 }
1133}
1134
1136{
1137 bool full_update;
1138
1139 if (nfes && nfes != fes)
1140 {
1141 full_update = true;
1142 fes = nfes;
1143 }
1144 else
1145 {
1146 // Check for different size (e.g. assembled form on non-conforming space)
1147 // or different sequence number.
1148 full_update = (fes->GetVSize() != Height() ||
1149 sequence < fes->GetSequence());
1150 }
1151
1152 delete mat_e;
1153 mat_e = NULL;
1155 delete static_cond;
1156 static_cond = NULL;
1157
1158 if (full_update)
1159 {
1160 delete mat;
1161 mat = NULL;
1162 delete hybridization;
1163 hybridization = NULL;
1165 }
1166 else
1167 {
1168 if (mat) { *mat = 0.0; }
1169 if (hybridization) { hybridization->Reset(); }
1170 }
1171
1172 height = width = fes->GetVSize();
1173
1174 if (ext) { ext->Update(); }
1175}
1176
1178{
1179 diag_policy = policy;
1180}
1181
1183{
1184 delete mat_e;
1185 delete mat;
1186 delete element_matrices;
1187 delete static_cond;
1188 delete hybridization;
1189
1190 if (!extern_bfs)
1191 {
1192 int k;
1193 for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
1194 for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
1195 for (k=0; k < interior_face_integs.Size(); k++)
1196 { delete interior_face_integs[k]; }
1197 for (k=0; k < boundary_face_integs.Size(); k++)
1198 { delete boundary_face_integs[k]; }
1199 }
1200
1201 delete ext;
1202}
1203
1204
1205MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1206 FiniteElementSpace *te_fes)
1207 : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1208{
1209 trial_fes = tr_fes;
1210 test_fes = te_fes;
1211 mat = NULL;
1212 mat_e = NULL;
1213 extern_bfs = 0;
1215 ext = NULL;
1216}
1217
1218MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1219 FiniteElementSpace *te_fes,
1220 MixedBilinearForm * mbf)
1221 : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1222{
1223 trial_fes = tr_fes;
1224 test_fes = te_fes;
1225 mat = NULL;
1226 mat_e = NULL;
1227 extern_bfs = 1;
1228 ext = NULL;
1229
1230 // Copy the pointers to the integrators
1233
1236
1238
1241
1243 ext = NULL;
1244}
1245
1247{
1248 if (ext)
1249 {
1250 MFEM_ABORT("the assembly level has already been set!");
1251 }
1252 assembly = assembly_level;
1253 switch (assembly)
1254 {
1256 break;
1258 // ext = new FAMixedBilinearFormExtension(this);
1259 // Use the original BilinearForm implementation for now
1260 break;
1262 mfem_error("Element assembly not supported yet... stay tuned!");
1263 // ext = new EAMixedBilinearFormExtension(this);
1264 break;
1267 break;
1269 mfem_error("Matrix-free action not supported yet... stay tuned!");
1270 // ext = new MFMixedBilinearFormExtension(this);
1271 break;
1272 default:
1273 mfem_error("Unknown assembly level");
1274 }
1275}
1276
1278{
1279 return (*mat)(i, j);
1280}
1281
1282const real_t & MixedBilinearForm::Elem (int i, int j) const
1283{
1284 return (*mat)(i, j);
1285}
1286
1287void MixedBilinearForm::Mult(const Vector & x, Vector & y) const
1288{
1289 y = 0.0;
1290 AddMult(x, y);
1291}
1292
1294 const real_t a) const
1295{
1296 if (ext)
1297 {
1298 ext->AddMult(x, y, a);
1299 }
1300 else
1301 {
1302 mat->AddMult(x, y, a);
1303 }
1304}
1305
1307{
1308 y = 0.0;
1309 AddMultTranspose(x, y);
1310}
1311
1313 const real_t a) const
1314{
1315 if (ext)
1316 {
1317 ext->AddMultTranspose(x, y, a);
1318 }
1319 else
1320 {
1321 mat->AddMultTranspose(x, y, a);
1322 }
1323}
1324
1326{
1328 {
1329 MFEM_WARNING("MixedBilinearForm::Inverse not possible with this "
1330 "assembly level!");
1331 return NULL;
1332 }
1333 else
1334 {
1335 return mat -> Inverse ();
1336 }
1337}
1338
1339void MixedBilinearForm::Finalize (int skip_zeros)
1340{
1342 {
1343 mat -> Finalize (skip_zeros);
1344 }
1345}
1346
1348{
1349 MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
1351 "MixedBilinearForm::GetBlocks: both trial and test spaces "
1352 "must use Ordering::byNODES!");
1353
1354 blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
1355
1356 mat->GetBlocks(blocks);
1357}
1358
1360{
1361 domain_integs.Append (bfi);
1362 domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
1363}
1364
1366 Array<int> &elem_marker)
1367{
1368 domain_integs.Append (bfi);
1369 domain_integs_marker.Append(&elem_marker);
1370}
1371
1373{
1374 boundary_integs.Append (bfi);
1375 boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
1376}
1377
1379 Array<int> &bdr_marker)
1380{
1381 boundary_integs.Append (bfi);
1382 boundary_integs_marker.Append(&bdr_marker);
1383}
1384
1386{
1387 trace_face_integs.Append (bfi);
1388}
1389
1391{
1392 boundary_trace_face_integs.Append(bfi);
1393 // NULL marker means apply everywhere
1395}
1396
1398 Array<int> &bdr_marker)
1399{
1400 boundary_trace_face_integs.Append(bfi);
1401 boundary_trace_face_integs_marker.Append(&bdr_marker);
1402}
1403
1405{
1406 if (ext)
1407 {
1408 ext->Assemble();
1409 return;
1410 }
1411
1412 ElementTransformation *eltrans;
1413 DofTransformation * dom_dof_trans;
1414 DofTransformation * ran_dof_trans;
1415 DenseMatrix elmat;
1416
1417 Mesh *mesh = test_fes -> GetMesh();
1418
1419 if (mat == NULL)
1420 {
1421 mat = new SparseMatrix(height, width);
1422 }
1423
1424 if (domain_integs.Size())
1425 {
1426 for (int k = 0; k < domain_integs.Size(); k++)
1427 {
1428 if (domain_integs_marker[k] != NULL)
1429 {
1430 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
1431 (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
1432 "invalid element marker for domain integrator #"
1433 << k << ", counting from zero");
1434 }
1435 }
1436
1437 for (int i = 0; i < test_fes -> GetNE(); i++)
1438 {
1439 const int elem_attr = mesh->GetAttribute(i);
1440 dom_dof_trans = trial_fes -> GetElementVDofs (i, trial_vdofs);
1441 ran_dof_trans = test_fes -> GetElementVDofs (i, test_vdofs);
1442 eltrans = test_fes -> GetElementTransformation (i);
1443
1445 elmat = 0.0;
1446 for (int k = 0; k < domain_integs.Size(); k++)
1447 {
1448 if (domain_integs_marker[k] == NULL ||
1449 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
1450 {
1451 domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1452 *test_fes -> GetFE(i),
1453 *eltrans, elemmat);
1454 elmat += elemmat;
1455 }
1456 }
1457 if (ran_dof_trans || dom_dof_trans)
1458 {
1459 TransformDual(ran_dof_trans, dom_dof_trans, elmat);
1460 }
1461 mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
1462 }
1463 }
1464
1465 if (boundary_integs.Size())
1466 {
1467 // Which boundary attributes need to be processed?
1468 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1469 mesh->bdr_attributes.Max() : 0);
1470 bdr_attr_marker = 0;
1471 for (int k = 0; k < boundary_integs.Size(); k++)
1472 {
1473 if (boundary_integs_marker[k] == NULL)
1474 {
1475 bdr_attr_marker = 1;
1476 break;
1477 }
1478 Array<int> &bdr_marker = *boundary_integs_marker[k];
1479 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1480 "invalid boundary marker for boundary integrator #"
1481 << k << ", counting from zero");
1482 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1483 {
1484 bdr_attr_marker[i] |= bdr_marker[i];
1485 }
1486 }
1487
1488 for (int i = 0; i < test_fes -> GetNBE(); i++)
1489 {
1490 const int bdr_attr = mesh->GetBdrAttribute(i);
1491 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1492
1493 dom_dof_trans = trial_fes -> GetBdrElementVDofs (i, trial_vdofs);
1494 ran_dof_trans = test_fes -> GetBdrElementVDofs (i, test_vdofs);
1495 eltrans = test_fes -> GetBdrElementTransformation (i);
1496
1498 elmat = 0.0;
1499 for (int k = 0; k < boundary_integs.Size(); k++)
1500 {
1501 if (boundary_integs_marker[k] &&
1502 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
1503
1504 boundary_integs[k]->AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1505 *test_fes -> GetBE(i),
1506 *eltrans, elemmat);
1507 elmat += elemmat;
1508 }
1509 if (ran_dof_trans || dom_dof_trans)
1510 {
1511 TransformDual(ran_dof_trans, dom_dof_trans, elmat);
1512 }
1513 mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
1514 }
1515 }
1516
1517 if (trace_face_integs.Size())
1518 {
1520 Array<int> test_vdofs2;
1521 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1522
1523 int nfaces = mesh->GetNumFaces();
1524 for (int i = 0; i < nfaces; i++)
1525 {
1526 ftr = mesh->GetFaceElementTransformations(i);
1529 trial_face_fe = trial_fes->GetFaceElement(i);
1530 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1531 if (ftr->Elem2No >= 0)
1532 {
1533 test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
1534 test_vdofs.Append(test_vdofs2);
1535 test_fe2 = test_fes->GetFE(ftr->Elem2No);
1536 }
1537 else
1538 {
1539 // The test_fe2 object is really a dummy and not used on the
1540 // boundaries, but we can't dereference a NULL pointer, and we don't
1541 // want to actually make a fake element.
1542 test_fe2 = test_fe1;
1543 }
1544 for (int k = 0; k < trace_face_integs.Size(); k++)
1545 {
1546 trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
1547 *test_fe2, *ftr, elemmat);
1549 }
1550 }
1551 }
1552
1553 if (boundary_trace_face_integs.Size())
1554 {
1556 Array<int> te_vdofs2;
1557 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1558
1559 // Which boundary attributes need to be processed?
1560 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1561 mesh->bdr_attributes.Max() : 0);
1562 bdr_attr_marker = 0;
1563 for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
1564 {
1565 if (boundary_trace_face_integs_marker[k] == NULL)
1566 {
1567 bdr_attr_marker = 1;
1568 break;
1569 }
1571 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1572 "invalid boundary marker for boundary trace face"
1573 "integrator #" << k << ", counting from zero");
1574 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1575 {
1576 bdr_attr_marker[i] |= bdr_marker[i];
1577 }
1578 }
1579
1580 for (int i = 0; i < trial_fes -> GetNBE(); i++)
1581 {
1582 const int bdr_attr = mesh->GetBdrAttribute(i);
1583 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1584
1585 ftr = mesh->GetBdrFaceTransformations(i);
1586 if (ftr)
1587 {
1588 const int iface = mesh->GetBdrElementFaceIndex(i);
1591 trial_face_fe = trial_fes->GetFaceElement(iface);
1592 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1593 // The test_fe2 object is really a dummy and not used on the
1594 // boundaries, but we can't dereference a NULL pointer, and we don't
1595 // want to actually make a fake element.
1596 test_fe2 = test_fe1;
1597 for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
1598 {
1600 (*boundary_trace_face_integs_marker[k])[bdr_attr-1] == 0)
1601 { continue; }
1602
1603 boundary_trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe,
1604 *test_fe1,
1605 *test_fe2,
1606 *ftr, elemmat);
1608 }
1609 }
1610 }
1611 }
1612}
1613
1615 Vector &diag) const
1616{
1617 if (ext)
1618 {
1619 MFEM_ASSERT(diag.Size() == test_fes->GetTrueVSize(),
1620 "Vector for holding diagonal has wrong size!");
1621 MFEM_ASSERT(D.Size() == trial_fes->GetTrueVSize(),
1622 "Vector for holding diagonal has wrong size!");
1623 const Operator *P_trial = trial_fes->GetProlongationMatrix();
1624 const Operator *P_test = test_fes->GetProlongationMatrix();
1625 if (!IsIdentityProlongation(P_trial))
1626 {
1627 Vector local_D(P_trial->Height());
1628 P_trial->Mult(D, local_D);
1629
1630 if (!IsIdentityProlongation(P_test))
1631 {
1632 Vector local_diag(P_test->Height());
1633 ext->AssembleDiagonal_ADAt(local_D, local_diag);
1634 P_test->MultTranspose(local_diag, diag);
1635 }
1636 else
1637 {
1638 ext->AssembleDiagonal_ADAt(local_D, diag);
1639 }
1640 }
1641 else
1642 {
1643 if (!IsIdentityProlongation(P_test))
1644 {
1645 Vector local_diag(P_test->Height());
1646 ext->AssembleDiagonal_ADAt(D, local_diag);
1647 P_test->MultTranspose(local_diag, diag);
1648 }
1649 else
1650 {
1651 ext->AssembleDiagonal_ADAt(D, diag);
1652 }
1653 }
1654 }
1655 else
1656 {
1657 MFEM_ABORT("Not implemented. Maybe assemble your bilinear form into a "
1658 "matrix and use SparseMatrix functions?");
1659 }
1660}
1661
1663{
1665 {
1666 MFEM_WARNING("Conforming assemble not supported for this assembly level!");
1667 return;
1668 }
1669
1670 Finalize();
1671
1673 if (P2)
1674 {
1675 SparseMatrix *R = Transpose(*P2);
1676 SparseMatrix *RA = mfem::Mult(*R, *mat);
1677 delete R;
1678 delete mat;
1679 mat = RA;
1680 }
1681
1683 if (P1)
1684 {
1685 SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1686 delete mat;
1687 mat = RAP;
1688 }
1689
1690 height = mat->Height();
1691 width = mat->Width();
1692}
1693
1694
1696{
1697 if (domain_integs.Size())
1698 {
1699 const FiniteElement &trial_fe = *trial_fes->GetFE(i);
1700 const FiniteElement &test_fe = *test_fes->GetFE(i);
1702 domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1703 elmat);
1704 for (int k = 1; k < domain_integs.Size(); k++)
1705 {
1706 domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1707 elemmat);
1708 elmat += elemmat;
1709 }
1710 }
1711 else
1712 {
1716 elmat = 0.0;
1717 }
1718}
1719
1721{
1722 if (boundary_integs.Size())
1723 {
1724 const FiniteElement &trial_be = *trial_fes->GetBE(i);
1725 const FiniteElement &test_be = *test_fes->GetBE(i);
1727 boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1728 elmat);
1729 for (int k = 1; k < boundary_integs.Size(); k++)
1730 {
1731 boundary_integs[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1732 elemmat);
1733 elmat += elemmat;
1734 }
1735 }
1736 else
1737 {
1741 elmat = 0.0;
1742 }
1743}
1744
1746 int i, const DenseMatrix &elmat, int skip_zeros)
1747{
1748 AssembleElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1749}
1750
1752 int i, const DenseMatrix &elmat, Array<int> &trial_vdofs_,
1753 Array<int> &test_vdofs_, int skip_zeros)
1754{
1755 trial_fes->GetElementVDofs(i, trial_vdofs_);
1756 test_fes->GetElementVDofs(i, test_vdofs_);
1757 if (mat == NULL)
1758 {
1759 mat = new SparseMatrix(height, width);
1760 }
1761 mat->AddSubMatrix(test_vdofs_, trial_vdofs_, elmat, skip_zeros);
1762}
1763
1765 int i, const DenseMatrix &elmat, int skip_zeros)
1766{
1767 AssembleBdrElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1768}
1769
1771 int i, const DenseMatrix &elmat, Array<int> &trial_vdofs_,
1772 Array<int> &test_vdofs_, int skip_zeros)
1773{
1774 trial_fes->GetBdrElementVDofs(i, trial_vdofs_);
1775 test_fes->GetBdrElementVDofs(i, test_vdofs_);
1776 if (mat == NULL)
1777 {
1778 mat = new SparseMatrix(height, width);
1779 }
1780 mat->AddSubMatrix(test_vdofs_, trial_vdofs_, elmat, skip_zeros);
1781}
1782
1784 const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
1785{
1786 int i, j, k;
1787 Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1788
1789 cols_marker = 0;
1790 for (i = 0; i < trial_fes -> GetNBE(); i++)
1791 if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1792 {
1793 trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1794 for (j = 0; j < tr_vdofs.Size(); j++)
1795 {
1796 if ( (k = tr_vdofs[j]) < 0 )
1797 {
1798 k = -1-k;
1799 }
1800 cols_marker[k] = 1;
1801 }
1802 }
1803 mat -> EliminateCols (cols_marker, &sol, &rhs);
1804}
1805
1807 const Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
1808{
1809 mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1810}
1811
1813{
1814 int i, j, k;
1815 Array<int> te_vdofs;
1816
1817 for (i = 0; i < test_fes -> GetNBE(); i++)
1818 if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1819 {
1820 test_fes -> GetBdrElementVDofs (i, te_vdofs);
1821 for (j = 0; j < te_vdofs.Size(); j++)
1822 {
1823 if ( (k = te_vdofs[j]) < 0 )
1824 {
1825 k = -1-k;
1826 }
1827 mat -> EliminateRow (k);
1828 }
1829 }
1830}
1831
1833 const Array<int> &trial_tdof_list,
1834 const Array<int> &test_tdof_list,
1835 OperatorHandle &A)
1836
1837{
1838 if (ext)
1839 {
1840 ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
1841 return;
1842 }
1843
1846
1847 mat->Finalize();
1848
1849 if (test_P && trial_P)
1850 {
1851 SparseMatrix *m = RAP(*test_P, *mat, *trial_P);
1852 delete mat;
1853 mat = m;
1854 }
1855 else if (test_P)
1856 {
1857 SparseMatrix *m = TransposeMult(*test_P, *mat);
1858 delete mat;
1859 mat = m;
1860 }
1861 else if (trial_P)
1862 {
1863 SparseMatrix *m = mfem::Mult(*mat, *trial_P);
1864 delete mat;
1865 mat = m;
1866 }
1867
1868 Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
1870 ess_trial_tdof_marker);
1872 ess_test_tdof_marker);
1873
1874 mat_e = new SparseMatrix(mat->Height(), mat->Width());
1875 mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
1876
1877 for (int i=0; i<test_tdof_list.Size(); ++i)
1878 {
1879 mat->EliminateRow(test_tdof_list[i]);
1880 }
1881 mat_e->Finalize();
1882 A.Reset(mat, false);
1883}
1884
1886 const Array<int> &trial_tdof_list,
1887 const Array<int> &test_tdof_list,
1888 Vector &x, Vector &b,
1889 OperatorHandle &A,
1890 Vector &X, Vector &B)
1891{
1892 if (ext)
1893 {
1894 ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
1895 x, b, A, X, B);
1896 return;
1897 }
1898
1899 const Operator *Pi = this->GetProlongation();
1900 const Operator *Po = this->GetOutputProlongation();
1901 const Operator *Ri = this->GetRestriction();
1902 InitTVectors(Po, Ri, Pi, x, b, X, B);
1903
1904 if (!mat_e)
1905 {
1906 FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list,
1907 A); // Set A = mat_e
1908 }
1909 // Eliminate essential BCs with B -= Ab xb
1910 mat_e->AddMult(X, B, -1.0);
1911
1912 B.SetSubVector(test_tdof_list, 0.0);
1913}
1914
1916{
1917 delete mat;
1918 mat = NULL;
1919 delete mat_e;
1920 mat_e = NULL;
1923 if (ext) { ext->Update(); }
1924}
1925
1927{
1928 if (mat) { delete mat; }
1929 if (mat_e) { delete mat_e; }
1930 if (!extern_bfs)
1931 {
1932 int i;
1933 for (i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
1934 for (i = 0; i < boundary_integs.Size(); i++)
1935 { delete boundary_integs[i]; }
1936 for (i = 0; i < trace_face_integs.Size(); i++)
1937 { delete trace_face_integs[i]; }
1938 for (i = 0; i < boundary_trace_face_integs.Size(); i++)
1939 { delete boundary_trace_face_integs[i]; }
1940 }
1941 delete ext;
1942}
1943
1945{
1946 if (ext)
1947 {
1948 MFEM_ABORT("the assembly level has already been set!");
1949 }
1950 assembly = assembly_level;
1951 switch (assembly)
1952 {
1955 // Use the original implementation for now
1956 break;
1958 mfem_error("Element assembly not supported yet... stay tuned!");
1959 break;
1962 break;
1964 mfem_error("Matrix-free action not supported yet... stay tuned!");
1965 break;
1966 default:
1967 mfem_error("Unknown assembly level");
1968 }
1969}
1970
1972{
1973 if (ext)
1974 {
1975 ext->Assemble();
1976 return;
1977 }
1978
1979 ElementTransformation *eltrans;
1980 DofTransformation * dom_dof_trans;
1981 DofTransformation * ran_dof_trans;
1982 DenseMatrix elmat;
1983
1984 Mesh *mesh = test_fes->GetMesh();
1985
1986 if (mat == NULL)
1987 {
1988 mat = new SparseMatrix(height, width);
1989 }
1990
1991 if (domain_integs.Size())
1992 {
1993 for (int k = 0; k < domain_integs.Size(); k++)
1994 {
1995 if (domain_integs_marker[k] != NULL)
1996 {
1997 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
1998 (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
1999 "invalid element marker for domain integrator #"
2000 << k << ", counting from zero");
2001 }
2002 }
2003
2004 for (int i = 0; i < test_fes->GetNE(); i++)
2005 {
2006 const int elem_attr = mesh->GetAttribute(i);
2007 dom_dof_trans = trial_fes->GetElementVDofs(i, trial_vdofs);
2008 ran_dof_trans = test_fes->GetElementVDofs(i, test_vdofs);
2009 eltrans = test_fes->GetElementTransformation(i);
2010
2012 elmat = 0.0;
2013 for (int k = 0; k < domain_integs.Size(); k++)
2014 {
2015 if (domain_integs_marker[k] == NULL ||
2016 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
2017 {
2018 domain_integs[k]->AssembleElementMatrix2(*trial_fes->GetFE(i),
2019 *test_fes->GetFE(i),
2020 *eltrans, elemmat);
2021 elmat += elemmat;
2022 }
2023 }
2024 if (ran_dof_trans || dom_dof_trans)
2025 {
2026 TransformPrimal(ran_dof_trans, dom_dof_trans, elemmat);
2027 }
2029 }
2030 }
2031
2032 if (trace_face_integs.Size())
2033 {
2034 const int nfaces = test_fes->GetMesh()->GetNumFaces();
2035 for (int i = 0; i < nfaces; i++)
2036 {
2039 eltrans = test_fes->GetMesh()->GetFaceTransformation(i);
2040
2042 elmat = 0.0;
2043 for (int k = 0; k < trace_face_integs.Size(); k++)
2044 {
2045 trace_face_integs[k]->AssembleElementMatrix2(*trial_fes->GetFaceElement(i),
2047 *eltrans, elemmat);
2048 elmat += elemmat;
2049 }
2050 mat->SetSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
2051 }
2052 }
2053}
2054
2055}
Dynamic 2D array using row-major layout.
Definition array.hpp:372
void SetSize(int m, int n)
Definition array.hpp:385
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition array.hpp:321
int Size() const
Return the logical size of the array.
Definition array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
virtual void Assemble()=0
Assemble at the level given for the BilinearFormExtension subclass.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)=0
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)=0
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Abstract base class BilinearFormIntegrator.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets Operator::DiagonalPolicy used upon construction of the linear system. Policies include:
Array< BilinearFormIntegrator * > domain_integs
Set of Domain Integrators to be applied.
Array< Array< int > * > domain_integs_marker
Entries are not owned.
BilinearFormExtension * ext
Extension for supporting Full Assembly (FA), Element Assembly (EA),Partial Assembly (PA),...
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a tdof Vector.
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
Array< Array< int > * > boundary_face_integs_marker
Entries are not owned.
void UseSparsity(int *I, int *J, bool isSorted)
Use the given CSR sparsity pattern to allocate the internal SparseMatrix.
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
Use the stored eliminated part of the matrix (see EliminateVDofs(const Array<int> &,...
virtual void MultTranspose(const Vector &x, Vector &y) const
Matrix transpose vector multiplication: .
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void ComputeElementMatrices()
Compute and store internally all element matrices.
Array< BilinearFormIntegrator * > boundary_face_integs
Set of boundary face Integrators to be applied.
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
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 FreeElementMatrices()
Free the memory used by the element matrices.
virtual ~BilinearForm()
Deletes internal matrices, bilinear integrators, and the BilinearFormExtension.
BilinearForm()
may be used in the construction of derived classes
Hybridization * hybridization
Owned.
DenseTensor * element_matrices
Owned.
void EliminateEssentialBCFromDofsDiag(const Array< int > &ess_dofs, real_t value)
Perform elimination and set the diagonal entry to the given value.
DiagonalPolicy diag_policy
This data member allows one to specify what should be done to the diagonal matrix entries and corresp...
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator. Assumes ownership of bfi.
void AllocMat()
Allocate appropriate SparseMatrix and assign it to mat.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse: (currently returns NULL)
void EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
Enable hybridization.
int extern_bfs
Indicates the BilinearFormIntegrators stored in domain_integs, boundary_integs, interior_face_integs,...
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
void AssembleElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given element matrix.
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
virtual real_t & Elem(int i, int j)
Returns a reference to: .
long sequence
Indicates the Mesh::sequence corresponding to the current state of the BilinearForm.
Array< BilinearFormIntegrator * > interior_face_integs
Set of interior face Integrators to be applied.
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix transpose vector multiplication: .
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
Compute the boundary element matrix of the given boundary element.
void EliminateEssentialBCFromDofs(const Array< int > &ess_dofs, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateVDofs(const Array<int> &, const Vector &, Vector &,...
SparseMatrix * mat
Sparse matrix to be associated with the form. Owned.
AssemblyLevel assembly
The AssemblyLevel of the form (AssemblyLevel::LEGACY, AssemblyLevel::FULL, AssemblyLevel::ELEMENT,...
FiniteElementSpace * fes
FE space on which the form lives. Not owned.
int batch
Element batch size used in the form action (1, 8, num_elems, etc.)
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(....
int Size() const
Get the size of the BilinearForm as a square matrix.
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
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.
void ConformingAssemble()
For partially conforming trial and/or test FE spaces, complete the assembly process by performing wh...
Array< BilinearFormIntegrator * > boundary_integs
Set of Boundary Integrators to be applied.
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, real_t value)
Perform elimination and set the diagonal entry to the given value.
StaticCondensation * static_cond
Owned.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
SparseMatrix * mat_e
Sparse Matrix used to store the eliminations from the b.c. Owned. .
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void ClearExternalData()
Definition densemat.hpp:95
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition densemat.hpp:102
Rank 3 tensor (array of matrices)
real_t * GetData(int k)
int SizeJ() const
int SizeI() const
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
void TransformDual(real_t *v) const
Definition doftrans.cpp:81
Data and methods for element-assembled bilinear forms.
Data and methods for fully-assembled bilinear forms.
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:484
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1309
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element,...
Definition fespace.hpp:1136
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition fespace.cpp:3204
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:773
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:620
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
Definition fespace.cpp:667
const NURBSExtension * GetNURBSext() const
Definition fespace.hpp:561
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition fespace.hpp:597
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:287
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:725
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition fespace.hpp:782
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:740
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition fespace.cpp:316
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
Definition fespace.cpp:3237
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1302
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
Definition fespace.cpp:303
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
Definition fespace.cpp:514
void GetPatchVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom in vdofs for NURBS patch i.
Definition fespace.cpp:310
Abstract class for all finite elements.
Definition fe_base.hpp:239
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
Auxiliary class Hybridization, used to implement BilinearForm hybridization.
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
SparseMatrix & GetMatrix()
Return the serial hybridized matrix.
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
void Finalize()
Finalize the construction of the hybridized matrix.
void ReduceRHS(const Vector &b, Vector &b_r) const
void SetConstraintIntegrator(BilinearFormIntegrator *c_integ)
A standard isoparametric element transformation.
Definition eltrans.hpp:363
Data and methods for matrix-free bilinear forms.
Abstract data type for matrix inverse.
Definition matrix.hpp:63
Abstract data type matrix.
Definition matrix.hpp:28
Class used by MFEM to store pointers to host and/or device memory.
Mesh data type.
Definition mesh.hpp:56
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition mesh.cpp:980
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6250
int GetAttribute(int i) const
Return the attribute of element i.
Definition mesh.hpp:1333
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition mesh.hpp:1339
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
Definition mesh.hpp:1518
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.
Definition mesh.cpp:583
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition mesh.cpp:1099
Table * GetFaceToElementTable() const
Definition mesh.cpp:7167
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
virtual void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const =0
virtual void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)=0
virtual void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)=0
virtual void MultTranspose(const Vector &x, Vector &y) const
Matrix transpose vector multiplication: .
void ConformingAssemble()
For partially conforming trial and/or test FE spaces, complete the assembly process by performing wh...
MixedBilinearFormExtension * ext
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
virtual const Operator * GetProlongation() const
Get the input finite element space prolongation matrix.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds a boundary integrator. Assumes ownership of bfi.
virtual void EliminateTestDofs(const Array< int > &bdr_attr_is_ess)
Eliminate essential boundary DOFs from the rows of the system.
virtual void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Return in A that is column-constrained.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
int extern_bfs
Indicates the BilinearFormIntegrators stored in MixedBilinearForm::domain_integs, MixedBilinearForm::...
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T into diag, where A is this mixed bilinear form and D is a diagonal.
void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi)
Adds a boundary trace face integrator. Assumes ownership of bfi.
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
Array< BilinearFormIntegrator * > trace_face_integs
Trace face (skeleton) integrators.
virtual const Operator * GetRestriction() const
Get the input finite element space restriction matrix.
Array< BilinearFormIntegrator * > boundary_trace_face_integs
Boundary trace face (skeleton) integrators.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
Array< BilinearFormIntegrator * > domain_integs
Domain integrators.
SparseMatrix * mat
Owned.
void Update()
Must be called after making changes to trial_fes or test_fes.
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix vector multiple to a vector: .
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
Array< Array< int > * > domain_integs_marker
Entries are not owned.
Array< BilinearFormIntegrator * > boundary_integs
Boundary integrators.
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Eliminate the list of DOFs from the columns of the system.
void EliminateTrialDofs(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
Eliminate essential boundary DOFs from the columns of the system.
void AssembleElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given element matrix.
Array< Array< int > * > boundary_trace_face_integs_marker
Entries are not owned.
virtual void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)
Form the linear system A X = B, corresponding to this mixed bilinear form and the linear form b(....
virtual real_t & Elem(int i, int j)
Returns a reference to: .
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
Compute the boundary element matrix of the given boundary element.
SparseMatrix * mat_e
Owned.
FiniteElementSpace * trial_fes
Not owned.
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Extract the associated matrix as SparseMatrix blocks. The number of block rows and columns is given b...
void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi)
Add a trace face integrator. Assumes ownership of bfi.
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix transpose vector multiplication: .
virtual ~MixedBilinearForm()
Deletes internal matrices, bilinear integrators, and the BilinearFormExtension.
FiniteElementSpace * test_fes
Not owned.
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse: (currently unimplemented and returns NUL...
virtual const Operator * GetOutputProlongation() const
Get the test finite element space prolongation matrix.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
AssemblyLevel assembly
The form assembly level (full, partial, etc.)
int GetNP() const
Definition nurbs.hpp:482
Pointer to an Operator of a specified type.
Definition handle.hpp:34
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition handle.hpp:145
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition operator.hpp:48
@ DIAG_ONE
Set the diagonal value to one.
Definition operator.hpp:50
@ DIAG_KEEP
Keep the diagonal value.
Definition operator.hpp:51
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition operator.cpp:148
void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi, Vector &x, Vector &b, Vector &X, Vector &B) const
Initializes memory for true vectors of linear system.
Definition operator.cpp:22
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
Operator transpose application: y+=A^t(x) (default) or y+=a*A^t(x).
Definition operator.cpp:58
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition operator.hpp:93
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
Operator application: y+=A(x) (default) or y+=a*A(x).
Definition operator.cpp:51
Data and methods for partially-assembled bilinear forms.
Partial assembly extension for DiscreteLinearOperator.
Data and methods for partially-assembled mixed bilinear forms.
Data type sparse matrix.
Definition sparsemat.hpp:51
void GetDiag(Vector &d) const
Returns the Diagonal of A.
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Extract all column indices and values from a given row.
bool Finalized() const
Returns whether or not CSR format has been finalized.
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
y += At * x (default) or y += a * At * x
void AbsMultTranspose(const Vector &x, Vector &y) const
y = |At| * x, using entry-wise absolute values of the transpose of matrix A
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
void EliminateRow(int row, const real_t sol, Vector &rhs)
Eliminates a column from the transpose matrix.
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR.
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
y += A * x (default) or y += a * A * x
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
int * GetJ()
Return the array J.
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
int * GetI()
Return the array I.
void AssembleMatrix(int el, const DenseMatrix &elmat)
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
void Finalize()
Finalize the construction of the Schur complement matrix.
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
bool ReducesTrueVSize() const
void Init(bool symmetric, bool block_diagonal)
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
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....
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix.
void LoseData()
Call this if data has been stolen.
Definition table.hpp:180
int * GetJ()
Definition table.hpp:114
int * GetI()
Definition table.hpp:113
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition table.cpp:199
Vector data type.
Definition vector.hpp:80
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition vector.hpp:256
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
Definition vector.cpp:739
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
Mesh * GetMesh(int type)
Definition ex29.cpp:218
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void mfem_error(const char *msg)
Definition error.cpp:154
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
void TransformDual(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition doftrans.cpp:160
void TransformPrimal(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition doftrans.cpp:145
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
bool IsIdentityProlongation(const Operator *P)
Definition operator.hpp:720
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
real_t p(const Vector &x, real_t t)