MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
bilinearform.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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;
76
78 batch = 1;
79}
80
108
110{
111 if (ext)
112 {
113 MFEM_ABORT("the assembly level has already been set!");
114 }
115 assembly = assembly_level;
116 switch (assembly)
117 {
119 break;
121 SetDiagonalPolicy( DIAG_ONE ); // Only diagonal policy supported on device
122 ext.reset(new FABilinearFormExtension(this));
123 break;
125 ext.reset(new EABilinearFormExtension(this));
126 break;
128 ext.reset(new PABilinearFormExtension(this));
129 break;
131 ext.reset(new MFBilinearFormExtension(this));
132 break;
133 default:
134 MFEM_ABORT("BilinearForm: unknown assembly level");
135 }
136}
137
139{
141 {
142 static_cond.reset();
143 MFEM_WARNING("Static condensation not supported for this assembly level");
144 return;
145 }
147 if (static_cond->ReducesTrueVSize())
148 {
149 bool symmetric = false; // TODO
150 bool block_diagonal = false; // TODO
151 static_cond->Init(symmetric, block_diagonal);
152 }
153 else
154 {
155 static_cond.reset();
156 }
157}
158
160 BilinearFormIntegrator *constr_integ,
161 const Array<int> &ess_tdof_list)
162{
164 {
165 delete constr_integ;
166 hybridization.reset();
167 MFEM_WARNING("Hybridization not supported for this assembly level");
168 return;
169 }
170 hybridization.reset(new Hybridization(fes, constr_space));
172 {
173 hybridization->EnableDeviceExecution();
174 }
175 hybridization->SetConstraintIntegrator(constr_integ);
176 hybridization->Init(ess_tdof_list);
177}
178
179void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
180{
181 if (static_cond) { return; }
182
183 if (mat)
184 {
185 if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
186 {
187 return; // mat is already using the given sparsity
188 }
189 delete mat;
190 }
191 height = width = fes->GetVSize();
192 mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
193}
194
196{
197 MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
198 "invalid matrix A dimensions: "
199 << A.Height() << " x " << A.Width());
200 MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
201
202 UseSparsity(A.GetI(), A.GetJ(), A.ColumnsAreSorted());
203}
204
206{
207 return mat -> Elem(i,j);
208}
209
210const real_t& BilinearForm::Elem (int i, int j) const
211{
212 return mat -> Elem(i,j);
213}
214
216{
217 return mat -> Inverse();
218}
219
220void BilinearForm::Finalize (int skip_zeros)
221{
223 {
224 if (!static_cond) { mat->Finalize(skip_zeros); }
225 if (mat_e) { mat_e->Finalize(skip_zeros); }
226 if (static_cond) { static_cond->Finalize(); }
227 }
228 if (hybridization) { hybridization->Finalize(); }
229}
230
232{
233 domain_integs.Append(bfi);
234 domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
235}
236
238 Array<int> &elem_marker)
239{
240 domain_integs.Append(bfi);
241 domain_integs_marker.Append(&elem_marker);
242}
243
245{
246 boundary_integs.Append (bfi);
247 boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
248}
249
251 Array<int> &bdr_marker)
252{
253 boundary_integs.Append (bfi);
254 boundary_integs_marker.Append(&bdr_marker);
255}
256
261
263{
264 boundary_face_integs.Append(bfi);
265 // NULL marker means apply everywhere
266 boundary_face_integs_marker.Append(NULL);
267}
268
270 Array<int> &bdr_marker)
271{
272 boundary_face_integs.Append(bfi);
273 boundary_face_integs_marker.Append(&bdr_marker);
274}
275
277{
279 {
280 elmat.SetSize(element_matrices->SizeI(), element_matrices->SizeJ());
281 elmat = element_matrices->GetData(i);
282 return;
283 }
284
285 const FiniteElement &fe = *fes->GetFE(i);
286
287 if (domain_integs.Size())
288 {
290 domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
291 for (int k = 1; k < domain_integs.Size(); k++)
292 {
293 domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
294 elmat += elemmat;
295 }
296 }
297 else
298 {
299 const int ndof = fe.GetDof() * fes->GetVDim();
300 elmat.SetSize(ndof);
301 elmat = 0.0;
302 }
303}
304
306{
307 const FiniteElement &be = *fes->GetBE(i);
308
309 if (boundary_integs.Size())
310 {
312 boundary_integs[0]->AssembleElementMatrix(be, *eltrans, elmat);
313 for (int k = 1; k < boundary_integs.Size(); k++)
314 {
315 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
316 elmat += elemmat;
317 }
318 }
319 else
320 {
321 const int ndof = be.GetDof() * fes->GetVDim();
322 elmat.SetSize(ndof);
323 elmat = 0.0;
324 }
325}
326
328{
330 Mesh *mesh = fes -> GetMesh();
331 tr = mesh -> GetFaceElementTransformations (i);
332
333 const FiniteElement *fe1, *fe2;
334 fe1 = fes->GetFE(tr->Elem1No);
335 if (tr->Elem2No >= 0)
336 {
337 fe2 = fes->GetFE(tr->Elem2No);
338 }
339 else
340 {
341 // The fe2 object is really a dummy and not used on the
342 // boundaries, but we can't dereference a NULL pointer, and we don't
343 // want to actually make a fake element.
344 fe2 = fe1;
345 }
346
347 if (interior_face_integs.Size())
348 {
349 interior_face_integs[0] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elmat);
350 for (int k = 1; k < interior_face_integs.Size(); k++)
351 {
352 interior_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
353 elmat += elemmat;
354 }
355 }
356 else
357 {
358 int ndof = fe1->GetDof() * fes->GetVDim();
359 if (tr->Elem2No >= 0)
360 {
361 ndof += fe2->GetDof() * fes->GetVDim();
362 }
363
364 elmat.SetSize(ndof);
365 elmat = 0.0;
366 }
367}
368
370{
372 Mesh *mesh = fes -> GetMesh();
373 tr = mesh -> GetBdrFaceTransformations (i);
374
375 const FiniteElement *fe1, *fe2;
376
377 fe1 = fes -> GetFE (tr -> Elem1No);
378 // The fe2 object is really a dummy and not used on the boundaries,
379 // but we can't dereference a NULL pointer, and we don't want to
380 // actually make a fake element.
381 fe2 = fe1;
382
383 if (boundary_face_integs.Size())
384 {
385 boundary_face_integs[0] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elmat);
386 for (int k = 1; k < boundary_face_integs.Size(); k++)
387 {
388 boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
389 elmat += elemmat;
390 }
391 }
392 else
393 {
394 int ndof = fe1->GetDof() * fes->GetVDim();
395 elmat.SetSize(ndof);
396 elmat = 0.0;
397 }
398}
399
401 int i, const DenseMatrix &elmat, int skip_zeros)
402{
403 AssembleElementMatrix(i, elmat, vdofs, skip_zeros);
404}
405
407 int i, const DenseMatrix &elmat, Array<int> &vdofs_, int skip_zeros)
408{
409 fes->GetElementVDofs(i, vdofs_);
410 if (static_cond)
411 {
412 static_cond->AssembleMatrix(i, elmat);
413 }
414 else
415 {
416 if (mat == NULL)
417 {
418 AllocMat();
419 }
420 mat->AddSubMatrix(vdofs_, vdofs_, elmat, skip_zeros);
421 if (hybridization)
422 {
423 hybridization->AssembleMatrix(i, elmat);
424 }
425 }
426}
427
429 int i, const DenseMatrix &elmat, int skip_zeros)
430{
431 AssembleBdrElementMatrix(i, elmat, vdofs, skip_zeros);
432}
433
435 int i, const DenseMatrix &elmat, Array<int> &vdofs_, int skip_zeros)
436{
437 fes->GetBdrElementVDofs(i, vdofs_);
438 if (static_cond)
439 {
440 static_cond->AssembleBdrMatrix(i, elmat);
441 }
442 else
443 {
444 if (mat == NULL)
445 {
446 AllocMat();
447 }
448 mat->AddSubMatrix(vdofs_, vdofs_, elmat, skip_zeros);
449 if (hybridization)
450 {
451 hybridization->AssembleBdrMatrix(i, elmat);
452 }
453 }
454}
455
456void BilinearForm::Assemble(int skip_zeros)
457{
458 if (ext)
459 {
460 ext->Assemble();
461 if (hybridization)
462 {
463 hybridization->AssembleElementMatrices(GetElementMatrices());
464 }
465 return;
466 }
467
468 ElementTransformation *eltrans;
469 DofTransformation * doftrans;
470 Mesh *mesh = fes -> GetMesh();
471 DenseMatrix elmat, *elmat_p;
472
473 if (mat == NULL)
474 {
475 AllocMat();
476 }
477
478#ifdef MFEM_USE_LEGACY_OPENMP
479 int free_element_matrices = 0;
480 if (!element_matrices)
481 {
483 free_element_matrices = 1;
484 }
485#endif
486
487 if (domain_integs.Size())
488 {
489 for (int k = 0; k < domain_integs.Size(); k++)
490 {
491 if (domain_integs_marker[k] != NULL)
492 {
493 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
494 (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
495 "invalid element marker for domain integrator #"
496 << k << ", counting from zero");
497 }
498
499 if (domain_integs[k]->Patchwise())
500 {
501 MFEM_VERIFY(fes->GetNURBSext(), "Patchwise integration requires a "
502 << "NURBS FE space");
503 }
504 }
505
506 // Element-wise integration
507 for (int i = 0; i < fes -> GetNE(); i++)
508 {
509 // Set both doftrans (potentially needed to assemble the element
510 // matrix) and vdofs, which is also needed when the element matrices
511 // are pre-assembled.
512 doftrans = fes->GetElementVDofs(i, vdofs);
514 {
515 elmat_p = &(*element_matrices)(i);
516 }
517 else
518 {
519 const int elem_attr = fes->GetMesh()->GetAttribute(i);
520 eltrans = fes->GetElementTransformation(i);
521
522 elmat.SetSize(0);
523 for (int k = 0; k < domain_integs.Size(); k++)
524 {
525 if (domain_integs_marker[k]) { domain_integs_marker[k]->HostRead(); }
526 if ((domain_integs_marker[k] == NULL ||
527 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
528 && !domain_integs[k]->Patchwise())
529 {
530 domain_integs[k]->AssembleElementMatrix(*fes->GetFE(i),
531 *eltrans, elemmat);
532 if (elmat.Size() == 0)
533 {
534 elmat = elemmat;
535 }
536 else
537 {
538 elmat += elemmat;
539 }
540 }
541 }
542 if (elmat.Size() == 0)
543 {
544 continue;
545 }
546 else
547 {
548 elmat_p = &elmat;
549 }
550 if (doftrans)
551 {
552 doftrans->TransformDual(elmat);
553 }
554 elmat_p = &elmat;
555 }
556 if (static_cond)
557 {
558 static_cond->AssembleMatrix(i, *elmat_p);
559 }
560 else
561 {
562 mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
563 if (hybridization)
564 {
565 hybridization->AssembleMatrix(i, *elmat_p);
566 }
567 }
568 }
569
570 // Patch-wise integration
571 if (fes->GetNURBSext())
572 {
573 for (int p=0; p<mesh->NURBSext->GetNP(); ++p)
574 {
575 bool vdofsSet = false;
576 for (int k = 0; k < domain_integs.Size(); k++)
577 {
578 if (domain_integs[k]->Patchwise())
579 {
580 if (!vdofsSet)
581 {
583 vdofsSet = true;
584 }
585
586 SparseMatrix* spmat = nullptr;
587 domain_integs[k]->AssemblePatchMatrix(p, *fes, spmat);
588 Array<int> cols;
589 Vector srow;
590
591 for (int r=0; r<spmat->Height(); ++r)
592 {
593 spmat->GetRow(r, cols, srow);
594 for (int i=0; i<cols.Size(); ++i)
595 {
596 cols[i] = vdofs[cols[i]];
597 }
598 mat->AddRow(vdofs[r], cols, srow);
599 }
600
601 delete spmat;
602 }
603 }
604 }
605 }
606 }
607
608 if (boundary_integs.Size())
609 {
610 // Which boundary attributes need to be processed?
611 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
612 mesh->bdr_attributes.Max() : 0);
613 bdr_attr_marker = 0;
614 for (int k = 0; k < boundary_integs.Size(); k++)
615 {
616 if (boundary_integs_marker[k] == NULL)
617 {
618 bdr_attr_marker = 1;
619 break;
620 }
621 Array<int> &bdr_marker = *boundary_integs_marker[k];
622 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
623 "invalid boundary marker for boundary integrator #"
624 << k << ", counting from zero");
625 for (int i = 0; i < bdr_attr_marker.Size(); i++)
626 {
627 bdr_attr_marker[i] |= bdr_marker[i];
628 }
629 }
630
631 for (int i = 0; i < fes -> GetNBE(); i++)
632 {
633 const int bdr_attr = mesh->GetBdrAttribute(i);
634 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
635
636 const FiniteElement &be = *fes->GetBE(i);
637 doftrans = fes -> GetBdrElementVDofs (i, vdofs);
638 eltrans = fes -> GetBdrElementTransformation (i);
639 int k = 0;
640 for (; k < boundary_integs.Size(); k++)
641 {
642 if (boundary_integs_marker[k] &&
643 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
644
645 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elmat);
646 k++;
647 break;
648 }
649 for (; k < boundary_integs.Size(); k++)
650 {
651 if (boundary_integs_marker[k] &&
652 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
653
654 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
655 elmat += elemmat;
656 }
657 if (doftrans)
658 {
659 doftrans->TransformDual(elmat);
660 }
661 elmat_p = &elmat;
662 if (!static_cond)
663 {
664 mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
665 if (hybridization)
666 {
667 hybridization->AssembleBdrMatrix(i, *elmat_p);
668 }
669 }
670 else
671 {
672 static_cond->AssembleBdrMatrix(i, *elmat_p);
673 }
674 }
675 }
676
677 if (interior_face_integs.Size())
678 {
680 Array<int> vdofs2;
681
682 int nfaces = mesh->GetNumFaces();
683 for (int i = 0; i < nfaces; i++)
684 {
685 tr = mesh -> GetInteriorFaceTransformations (i);
686 if (tr != NULL)
687 {
688 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
689 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
690 vdofs.Append (vdofs2);
691 for (int k = 0; k < interior_face_integs.Size(); k++)
692 {
694 AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
695 *fes->GetFE(tr->Elem2No),
696 *tr, elemmat);
697 mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
698 }
699 }
700 }
701 }
702
703 if (boundary_face_integs.Size())
704 {
706 const FiniteElement *fe1, *fe2;
707
708 // Which boundary attributes need to be processed?
709 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
710 mesh->bdr_attributes.Max() : 0);
711 bdr_attr_marker = 0;
712 for (int k = 0; k < boundary_face_integs.Size(); k++)
713 {
714 if (boundary_face_integs_marker[k] == NULL)
715 {
716 bdr_attr_marker = 1;
717 break;
718 }
719 Array<int> &bdr_marker = *boundary_face_integs_marker[k];
720 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
721 "invalid boundary marker for boundary face integrator #"
722 << k << ", counting from zero");
723 for (int i = 0; i < bdr_attr_marker.Size(); i++)
724 {
725 bdr_attr_marker[i] |= bdr_marker[i];
726 }
727 }
728
729 for (int i = 0; i < fes -> GetNBE(); i++)
730 {
731 const int bdr_attr = mesh->GetBdrAttribute(i);
732 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
733
734 tr = mesh -> GetBdrFaceTransformations (i);
735 if (tr != NULL)
736 {
737 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
738 fe1 = fes -> GetFE (tr -> Elem1No);
739 // The fe2 object is really a dummy and not used on the boundaries,
740 // but we can't dereference a NULL pointer, and we don't want to
741 // actually make a fake element.
742 fe2 = fe1;
743 for (int k = 0; k < boundary_face_integs.Size(); k++)
744 {
746 (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
747 { continue; }
748
749 boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
750 elemmat);
751 mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
752 }
753 }
754 }
755 }
756
757#ifdef MFEM_USE_LEGACY_OPENMP
758 if (free_element_matrices)
759 {
761 }
762#endif
763}
764
766{
767 // Do not remove zero entries to preserve the symmetric structure of the
768 // matrix which in turn will give rise to symmetric structure in the new
769 // matrix. This ensures that subsequent calls to EliminateRowCol will work
770 // correctly.
771 Finalize(0);
772 MFEM_ASSERT(mat, "the BilinearForm is not assembled");
773
775 if (!P) { return; } // conforming mesh
776
777 SparseMatrix *R = Transpose(*P);
778 SparseMatrix *RA = mfem::Mult(*R, *mat);
779 delete mat;
780 if (mat_e)
781 {
782 SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
783 delete mat_e;
784 mat_e = RAe;
785 }
786 delete R;
787 mat = mfem::Mult(*RA, *P);
788 delete RA;
789 if (mat_e)
790 {
791 SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
792 delete mat_e;
793 mat_e = RAeP;
794 }
795
796 height = mat->Height();
797 width = mat->Width();
798}
799
801{
802 MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
803 "Vector for holding diagonal has wrong size!");
805 if (!ext)
806 {
807 MFEM_ASSERT(mat, "the BilinearForm is not assembled!");
808 MFEM_ASSERT(cP == nullptr || mat->Height() == cP->Width(),
809 "BilinearForm::ConformingAssemble() is not called!");
810 mat->GetDiag(diag);
811 return;
812 }
813 // Here, we have extension, ext.
814 if (!cP)
815 {
816 ext->AssembleDiagonal(diag);
817 return;
818 }
819 // Here, we have extension, ext, and conforming prolongation, cP.
820
821 // For an AMR mesh, a convergent diagonal is assembled with |P^T| d_l,
822 // where |P^T| has the entry-wise absolute values of the conforming
823 // prolongation transpose operator.
824 Vector local_diag(cP->Height());
825 ext->AssembleDiagonal(local_diag);
826 cP->AbsMultTranspose(local_diag, diag);
827}
828
830 Vector &b, OperatorHandle &A, Vector &X,
831 Vector &B, int copy_interior)
832{
833 if (ext)
834 {
835 if (hybridization)
836 {
837 FormSystemMatrix(ess_tdof_list, A);
838 ConstrainedOperator A_constrained(this, ess_tdof_list);
839 A_constrained.EliminateRHS(x, b);
840 hybridization->ReduceRHS(b, B);
841 X.SetSize(B.Size());
842 X = 0.0;
843 }
844 else
845 {
846 ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
847 }
848 return;
849 }
851 FormSystemMatrix(ess_tdof_list, A);
852
853 // Transform the system and perform the elimination in B, based on the
854 // essential BC values from x. Restrict the BC part of x in X, and set the
855 // non-BC part to zero. Since there is no good initial guess for the Lagrange
856 // multipliers, set X = 0.0 for hybridization.
857 if (static_cond)
858 {
859 // Schur complement reduction to the exposed dofs
860 static_cond->ReduceSystem(x, b, X, B, copy_interior);
861 }
862 else if (!P) // conforming space
863 {
864 if (hybridization)
865 {
866 // Reduction to the Lagrange multipliers system
867 EliminateVDofsInRHS(ess_tdof_list, x, b);
868 hybridization->ReduceRHS(b, B);
869 X.SetSize(B.Size());
870 X = 0.0;
871 }
872 else
873 {
874 // A, X and B point to the same data as mat, x and b
875 EliminateVDofsInRHS(ess_tdof_list, x, b);
876 X.MakeRef(x, 0, x.Size());
877 B.MakeRef(b, 0, b.Size());
878 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
879 }
880 }
881 else // non-conforming space
882 {
883 if (hybridization)
884 {
885 // Reduction to the Lagrange multipliers system
887 Vector conf_b(P->Width()), conf_x(P->Width());
888 P->MultTranspose(b, conf_b);
889 R->Mult(x, conf_x);
890 EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
891 R->MultTranspose(conf_b, b); // store eliminated rhs in b
892 hybridization->ReduceRHS(conf_b, B);
893 X.SetSize(B.Size());
894 X = 0.0;
895 }
896 else
897 {
898 // Variational restriction with P
900 B.SetSize(P->Width());
901 P->MultTranspose(b, B);
902 X.SetSize(R->Height());
903 R->Mult(x, X);
904 EliminateVDofsInRHS(ess_tdof_list, X, B);
905 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
906 }
907 }
908}
909
912{
913 if (ext)
914 {
915 if (hybridization)
916 {
917 const int remove_zeros = 0;
918 Finalize(remove_zeros);
919 A.Reset(&hybridization->GetMatrix(), false);
920 }
921 else
922 {
923 ext->FormSystemMatrix(ess_tdof_list, A);
924 }
925 return;
926 }
927
928 // Finish the matrix assembly and perform BC elimination, storing the
929 // eliminated part of the matrix.
930 if (static_cond)
931 {
932 if (!static_cond->HasEliminatedBC())
933 {
934 static_cond->SetEssentialTrueDofs(ess_tdof_list);
935 static_cond->Finalize(); // finalize Schur complement (to true dofs)
936 static_cond->EliminateReducedTrueDofs(diag_policy);
937 static_cond->Finalize(); // finalize eliminated part
938 }
939 A.Reset(&static_cond->GetMatrix(), false);
940 }
941 else
942 {
943 if (!mat_e)
944 {
946 if (P) { ConformingAssemble(); }
947 EliminateVDofs(ess_tdof_list, diag_policy);
948 const int remove_zeros = 0;
949 Finalize(remove_zeros);
950 }
951 if (hybridization)
952 {
953 A.Reset(&hybridization->GetMatrix(), false);
954 }
955 else
956 {
957 A.Reset(mat, false);
958 }
959 }
960}
961
963 const Vector &b, Vector &x)
964{
965 if (ext && !hybridization)
966 {
967 ext->RecoverFEMSolution(X, b, x);
968 return;
969 }
970
972 if (!P) // conforming space
973 {
974 if (static_cond)
975 {
976 // Private dofs back solve
977 static_cond->ComputeSolution(b, X, x);
978 }
979 else if (hybridization)
980 {
981 // Primal unknowns recovery
982 hybridization->ComputeSolution(b, X, x);
983 }
984 else
985 {
986 // X and x point to the same data
987
988 // If the validity flags of X's Memory were changed (e.g. if it was
989 // moved to device memory) then we need to tell x about that.
990 x.SyncMemory(X);
991 }
992 }
993 else // non-conforming space
994 {
995 if (static_cond)
996 {
997 // Private dofs back solve
998 static_cond->ComputeSolution(b, X, x);
999 }
1000 else if (hybridization)
1001 {
1002 // Primal unknowns recovery
1003 Vector conf_b(P->Width()), conf_x(P->Width());
1004 P->MultTranspose(b, conf_b);
1006 R->Mult(x, conf_x); // get essential b.c. from x
1007 hybridization->ComputeSolution(conf_b, X, conf_x);
1008 x.SetSize(P->Height());
1009 P->Mult(conf_x, x);
1010 }
1011 else
1012 {
1013 // Apply conforming prolongation
1014 x.SetSize(P->Height());
1015 P->Mult(X, x);
1016 }
1017 }
1018}
1019
1021{
1022 if (element_matrices) { return; }
1023
1024 if (auto *ea_ext = dynamic_cast<EABilinearFormExtension*>(ext.get()))
1025 {
1026 element_matrices.reset(new DenseTensor);
1027 ea_ext->GetElementMatrices(*element_matrices, ElementDofOrdering::NATIVE, true);
1028 return;
1029 }
1030
1031 if (domain_integs.Size() == 0 || fes->GetNE() == 0)
1032 {
1033 element_matrices.reset(new DenseTensor);
1034 return;
1035 }
1036
1037 int num_elements = fes->GetNE();
1038 int num_dofs_per_el = fes->GetTypicalFE()->GetDof() * fes->GetVDim();
1039
1040 element_matrices.reset(new DenseTensor(num_dofs_per_el, num_dofs_per_el,
1041 num_elements));
1042
1043 DenseMatrix tmp;
1045
1046#ifdef MFEM_USE_LEGACY_OPENMP
1047 #pragma omp parallel for private(tmp,eltrans)
1048#endif
1049 for (int i = 0; i < num_elements; i++)
1050 {
1051 DenseMatrix elmat(element_matrices->GetData(i),
1052 num_dofs_per_el, num_dofs_per_el);
1053 const FiniteElement &fe = *fes->GetFE(i);
1054#ifdef MFEM_DEBUG
1055 if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
1056 mfem_error("BilinearForm::ComputeElementMatrices:"
1057 " all elements must have same number of dofs");
1058#endif
1059 fes->GetElementTransformation(i, &eltrans);
1060
1061 domain_integs[0]->AssembleElementMatrix(fe, eltrans, elmat);
1062 for (int k = 1; k < domain_integs.Size(); k++)
1063 {
1064 // note: some integrators may not be thread-safe
1065 domain_integs[k]->AssembleElementMatrix(fe, eltrans, tmp);
1066 elmat += tmp;
1067 }
1068 elmat.ClearExternalData();
1069 }
1070}
1071
1073{
1074 ComputeElementMatrices(); // Won't recompute if element_matrices exists
1075 return *element_matrices;
1076}
1077
1079 const Vector &sol, Vector &rhs,
1080 DiagonalPolicy dpolicy)
1081{
1082 Array<int> ess_dofs, conf_ess_dofs;
1083 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
1084
1085 if (fes->GetVSize() == height)
1086 {
1087 EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
1088 }
1089 else
1090 {
1091 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
1092 EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
1093 }
1094}
1095
1097 DiagonalPolicy dpolicy)
1098{
1099 Array<int> ess_dofs, conf_ess_dofs;
1100 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
1101
1102 if (fes->GetVSize() == height)
1103 {
1104 EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
1105 }
1106 else
1107 {
1108 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
1109 EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
1110 }
1111}
1112
1114 real_t value)
1115{
1116 Array<int> ess_dofs, conf_ess_dofs;
1117 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
1118
1119 if (fes->GetVSize() == height)
1120 {
1121 EliminateEssentialBCFromDofsDiag(ess_dofs, value);
1122 }
1123 else
1124 {
1125 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
1126 EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
1127 }
1128}
1129
1131 const Vector &sol, Vector &rhs,
1132 DiagonalPolicy dpolicy)
1133{
1134 vdofs_.HostRead();
1135 for (int i = 0; i < vdofs_.Size(); i++)
1136 {
1137 int vdof = vdofs_[i];
1138 if ( vdof >= 0 )
1140 mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
1141 }
1142 else
1143 {
1144 mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
1145 }
1146 }
1147}
1148
1150 DiagonalPolicy dpolicy)
1151{
1152 if (mat_e == NULL)
1153 {
1154 mat_e = new SparseMatrix(height);
1155 }
1156
1157 vdofs_.HostRead();
1158 for (int i = 0; i < vdofs_.Size(); i++)
1159 {
1160 int vdof = vdofs_[i];
1161 if ( vdof >= 0 )
1162 {
1163 mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
1164 }
1165 else
1166 {
1167 mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
1168 }
1169 }
1170}
1171
1173 const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
1174 DiagonalPolicy dpolicy)
1175{
1176 MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
1177 MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
1178 MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
1179
1180 for (int i = 0; i < ess_dofs.Size(); i++)
1181 if (ess_dofs[i] < 0)
1182 {
1183 mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
1184 }
1185}
1186
1188 DiagonalPolicy dpolicy)
1189{
1190 MFEM_ASSERT(ess_dofs.Size() == height,
1191 "incorrect dof Array size: " << ess_dofs.Size() << ' ' << height);
1192
1193 for (int i = 0; i < ess_dofs.Size(); i++)
1194 if (ess_dofs[i] < 0)
1195 {
1196 mat -> EliminateRowCol (i, dpolicy);
1197 }
1198}
1199
1201 real_t value)
1202{
1203 MFEM_ASSERT(ess_dofs.Size() == height,
1204 "incorrect dof Array size: " << ess_dofs.Size() << ' ' << height);
1205
1206 for (int i = 0; i < ess_dofs.Size(); i++)
1207 if (ess_dofs[i] < 0)
1208 {
1209 mat -> EliminateRowColDiag (i, value);
1210 }
1211}
1212
1214 const Array<int> &vdofs_, const Vector &x, Vector &b)
1215{
1216 mat_e->AddMult(x, b, -1.);
1217 mat->PartMult(vdofs_, x, b);
1218}
1219
1220void BilinearForm::Mult(const Vector &x, Vector &y) const
1221{
1222 if (ext)
1223 {
1224 ext->Mult(x, y);
1225 }
1226 else
1227 {
1228 mat->Mult(x, y);
1229 }
1230}
1231
1233{
1234 if (ext)
1235 {
1236 ext->MultTranspose(x, y);
1237 }
1238 else
1239 {
1240 y = 0.0;
1241 AddMultTranspose (x, y);
1242 }
1243}
1244
1246{
1247 bool full_update;
1248
1249 if (nfes && nfes != fes)
1250 {
1251 full_update = true;
1252 fes = nfes;
1253 }
1254 else
1255 {
1256 // Check for different size (e.g. assembled form on non-conforming space)
1257 // or different sequence number.
1258 full_update = (fes->GetVSize() != Height() ||
1259 sequence < fes->GetSequence());
1260 }
1261
1262 delete mat_e;
1263 mat_e = NULL;
1265 static_cond.reset();
1266
1267 if (full_update)
1268 {
1269 delete mat;
1270 mat = NULL;
1271 hybridization.reset();
1273 }
1274 else
1275 {
1276 if (mat) { *mat = 0.0; }
1277 if (hybridization) { hybridization->Reset(); }
1278 }
1279
1280 height = width = fes->GetVSize();
1281
1282 if (ext) { ext->Update(); }
1283}
1284
1286{
1287 diag_policy = policy;
1288}
1289
1291{
1292 delete mat_e;
1293 delete mat;
1294
1295 if (!extern_bfs)
1296 {
1297 int k;
1298 for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
1299 for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
1300 for (k=0; k < interior_face_integs.Size(); k++)
1301 { delete interior_face_integs[k]; }
1302 for (k=0; k < boundary_face_integs.Size(); k++)
1303 { delete boundary_face_integs[k]; }
1304 }
1305}
1306
1307
1308MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1309 FiniteElementSpace *te_fes)
1310 : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1311{
1312 trial_fes = tr_fes;
1313 test_fes = te_fes;
1314 mat = NULL;
1315 mat_e = NULL;
1316 extern_bfs = 0;
1318 ext = NULL;
1319}
1320
1321MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1322 FiniteElementSpace *te_fes,
1323 MixedBilinearForm * mbf)
1324 : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1325{
1326 trial_fes = tr_fes;
1327 test_fes = te_fes;
1328 mat = NULL;
1329 mat_e = NULL;
1330 extern_bfs = 1;
1331
1332 // Copy the pointers to the integrators
1335
1338
1340
1343
1345 ext = NULL;
1346}
1347
1349{
1350 if (ext)
1351 {
1352 MFEM_ABORT("the assembly level has already been set!");
1353 }
1354 assembly = assembly_level;
1355 switch (assembly)
1356 {
1358 break;
1360 // ext.reset(new FAMixedBilinearFormExtension(this));
1361 // Use the original BilinearForm implementation for now
1362 break;
1364 MFEM_ABORT("Element assembly not supported yet... stay tuned!");
1365 // ext.reset(new EAMixedBilinearFormExtension(this));
1366 break;
1368 ext.reset(new PAMixedBilinearFormExtension(this));
1369 break;
1371 MFEM_ABORT("Matrix-free action not supported yet... stay tuned!");
1372 // ext.reset(new MFMixedBilinearFormExtension(this));
1373 break;
1374 default:
1375 MFEM_ABORT("Unknown assembly level");
1376 }
1377}
1378
1380{
1381 return (*mat)(i, j);
1382}
1383
1384const real_t & MixedBilinearForm::Elem (int i, int j) const
1385{
1386 return (*mat)(i, j);
1387}
1388
1389void MixedBilinearForm::Mult(const Vector & x, Vector & y) const
1390{
1391 y = 0.0;
1392 AddMult(x, y);
1393}
1394
1396 const real_t a) const
1397{
1398 if (ext)
1399 {
1400 ext->AddMult(x, y, a);
1401 }
1402 else
1403 {
1404 mat->AddMult(x, y, a);
1405 }
1406}
1407
1409{
1410 y = 0.0;
1411 AddMultTranspose(x, y);
1412}
1413
1415 const real_t a) const
1416{
1417 if (ext)
1418 {
1419 ext->AddMultTranspose(x, y, a);
1420 }
1421 else
1422 {
1423 mat->AddMultTranspose(x, y, a);
1424 }
1425}
1426
1428{
1430 {
1431 MFEM_WARNING("MixedBilinearForm::Inverse not possible with this "
1432 "assembly level!");
1433 return NULL;
1434 }
1435 else
1436 {
1437 return mat -> Inverse ();
1438 }
1439}
1440
1441void MixedBilinearForm::Finalize (int skip_zeros)
1442{
1444 {
1445 mat -> Finalize (skip_zeros);
1446 }
1447}
1448
1450{
1451 MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
1453 "MixedBilinearForm::GetBlocks: both trial and test spaces "
1454 "must use Ordering::byNODES!");
1455
1456 blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
1457
1458 mat->GetBlocks(blocks);
1459}
1460
1462{
1463 domain_integs.Append(bfi);
1464 domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
1465}
1466
1468 Array<int> &elem_marker)
1469{
1470 domain_integs.Append(bfi);
1471 domain_integs_marker.Append(&elem_marker);
1472}
1473
1475{
1476 boundary_integs.Append(bfi);
1477 boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
1478}
1479
1481 Array<int> &bdr_marker)
1482{
1483 boundary_integs.Append(bfi);
1484 boundary_integs_marker.Append(&bdr_marker);
1485}
1486
1491
1493{
1494 boundary_face_integs.Append(bfi);
1495 boundary_face_integs_marker.Append(NULL); // NULL marker means apply everywhere
1496}
1497
1499 Array<int> &bdr_marker)
1500{
1501 boundary_face_integs.Append(bfi);
1502 boundary_face_integs_marker.Append(&bdr_marker);
1503}
1504
1509
1511{
1512 boundary_trace_face_integs.Append(bfi);
1513 // NULL marker means apply everywhere
1515}
1516
1523
1525{
1526 if (ext)
1527 {
1528 ext->Assemble();
1529 return;
1530 }
1531
1532 ElementTransformation *eltrans;
1533 DofTransformation * dom_dof_trans;
1534 DofTransformation * ran_dof_trans;
1535 DenseMatrix elmat;
1536
1537 Mesh *mesh = test_fes -> GetMesh();
1538
1539 if (mat == NULL)
1540 {
1541 mat = new SparseMatrix(height, width);
1542 }
1543
1544 if (domain_integs.Size())
1545 {
1546 for (int k = 0; k < domain_integs.Size(); k++)
1547 {
1548 if (domain_integs_marker[k] != NULL)
1549 {
1550 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
1551 (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
1552 "invalid element marker for domain integrator #"
1553 << k << ", counting from zero");
1554 }
1555 }
1556
1557 for (int i = 0; i < test_fes -> GetNE(); i++)
1558 {
1559 const int elem_attr = mesh->GetAttribute(i);
1560 dom_dof_trans = trial_fes -> GetElementVDofs (i, trial_vdofs);
1561 ran_dof_trans = test_fes -> GetElementVDofs (i, test_vdofs);
1562 eltrans = test_fes -> GetElementTransformation (i);
1563
1565 elmat = 0.0;
1566 for (int k = 0; k < domain_integs.Size(); k++)
1567 {
1568 if (domain_integs_marker[k] == NULL ||
1569 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
1570 {
1571 domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1572 *test_fes -> GetFE(i),
1573 *eltrans, elemmat);
1574 elmat += elemmat;
1575 }
1576 }
1577 if (ran_dof_trans || dom_dof_trans)
1578 {
1579 TransformDual(ran_dof_trans, dom_dof_trans, elmat);
1580 }
1581 mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
1582 }
1583 }
1584
1585 if (boundary_integs.Size())
1586 {
1587 // Which boundary attributes need to be processed?
1588 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1589 mesh->bdr_attributes.Max() : 0);
1590 bdr_attr_marker = 0;
1591 for (int k = 0; k < boundary_integs.Size(); k++)
1592 {
1593 if (boundary_integs_marker[k] == NULL)
1594 {
1595 bdr_attr_marker = 1;
1596 break;
1597 }
1598 Array<int> &bdr_marker = *boundary_integs_marker[k];
1599 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1600 "invalid boundary marker for boundary integrator #"
1601 << k << ", counting from zero");
1602 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1603 {
1604 bdr_attr_marker[i] |= bdr_marker[i];
1605 }
1606 }
1607
1608 for (int i = 0; i < test_fes -> GetNBE(); i++)
1609 {
1610 const int bdr_attr = mesh->GetBdrAttribute(i);
1611 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1612
1613 dom_dof_trans = trial_fes -> GetBdrElementVDofs (i, trial_vdofs);
1614 ran_dof_trans = test_fes -> GetBdrElementVDofs (i, test_vdofs);
1615 eltrans = test_fes -> GetBdrElementTransformation (i);
1616
1618 elmat = 0.0;
1619 for (int k = 0; k < boundary_integs.Size(); k++)
1620 {
1621 if (boundary_integs_marker[k] &&
1622 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
1623
1624 boundary_integs[k]->AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1625 *test_fes -> GetBE(i),
1626 *eltrans, elemmat);
1627 elmat += elemmat;
1628 }
1629 if (ran_dof_trans || dom_dof_trans)
1630 {
1631 TransformDual(ran_dof_trans, dom_dof_trans, elmat);
1632 }
1633 mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
1634 }
1635 }
1636
1637 if (interior_face_integs.Size())
1638 {
1640 Array<int> trial_vdofs2, test_vdofs2;
1641 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
1642
1643 int nfaces = mesh->GetNumFaces();
1644 for (int i = 0; i < nfaces; i++)
1645 {
1646 ftr = mesh->GetInteriorFaceTransformations(i);
1647 if (ftr != NULL)
1648 {
1651 trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
1652 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1653 if (ftr->Elem2No >= 0)
1654 {
1655 trial_fes->GetElementVDofs(ftr->Elem2No, trial_vdofs2);
1656 test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
1657 trial_vdofs.Append(trial_vdofs2);
1658 test_vdofs.Append(test_vdofs2);
1659 trial_fe2 = trial_fes->GetFE(ftr->Elem2No);
1660 test_fe2 = test_fes->GetFE(ftr->Elem2No);
1661 }
1662 else
1663 {
1664 // The test_fe2 object is really a dummy and not used on the
1665 // boundaries, but we can't dereference a NULL pointer, and we don't
1666 // want to actually make a fake element.
1667 trial_fe2 = trial_fe1;
1668 test_fe2 = test_fe1;
1669 }
1670 for (int k = 0; k < interior_face_integs.Size(); k++)
1671 {
1672 interior_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
1673 *test_fe2,
1674 *ftr, elemmat);
1676 }
1677 }
1678 }
1679 }
1680
1681 if (boundary_face_integs.Size())
1682 {
1684 Array<int> tr_vdofs2, te_vdofs2;
1685 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
1686
1687 // Which boundary attributes need to be processed?
1688 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1689 mesh->bdr_attributes.Max() : 0);
1690 bdr_attr_marker = 0;
1691 for (int k = 0; k < boundary_face_integs.Size(); k++)
1692 {
1693 if (boundary_face_integs_marker[k] == NULL)
1694 {
1695 bdr_attr_marker = 1;
1696 break;
1697 }
1698 Array<int> &bdr_marker = *boundary_face_integs_marker[k];
1699 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1700 "invalid boundary marker for boundary face integrator #"
1701 << k << ", counting from zero");
1702 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1703 {
1704 bdr_attr_marker[i] |= bdr_marker[i];
1705 }
1706 }
1707
1708 for (int i = 0; i < trial_fes -> GetNBE(); i++)
1709 {
1710 const int bdr_attr = mesh->GetBdrAttribute(i);
1711 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1712
1713 ftr = mesh -> GetBdrFaceTransformations (i);
1714 if (ftr != NULL)
1715 {
1718 trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
1719 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1720 // The test_fe2 object is really a dummy and not used on the
1721 // boundaries, but we can't dereference a NULL pointer, and we don't
1722 // want to actually make a fake element.
1723 trial_fe2 = trial_fe1;
1724 test_fe2 = test_fe1;
1725 for (int k = 0; k < boundary_face_integs.Size(); k++)
1726 {
1728 (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
1729
1730 boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
1731 *test_fe2,
1732 *ftr, elemmat);
1734 }
1735 }
1736 }
1737 }
1738
1739 if (trace_face_integs.Size())
1740 {
1742 Array<int> test_vdofs2;
1743 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1744
1745 int nfaces = mesh->GetNumFaces();
1746 for (int i = 0; i < nfaces; i++)
1747 {
1748 ftr = mesh->GetFaceElementTransformations(i);
1751 trial_face_fe = trial_fes->GetFaceElement(i);
1752 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1753 if (ftr->Elem2No >= 0)
1754 {
1755 test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
1756 test_vdofs.Append(test_vdofs2);
1757 test_fe2 = test_fes->GetFE(ftr->Elem2No);
1758 }
1759 else
1760 {
1761 // The test_fe2 object is really a dummy and not used on the
1762 // boundaries, but we can't dereference a NULL pointer, and we don't
1763 // want to actually make a fake element.
1764 test_fe2 = test_fe1;
1765 }
1766 for (int k = 0; k < trace_face_integs.Size(); k++)
1767 {
1768 trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
1769 *test_fe2, *ftr, elemmat);
1771 }
1772 }
1773 }
1774
1775 if (boundary_trace_face_integs.Size())
1776 {
1778 Array<int> te_vdofs2;
1779 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1780
1781 // Which boundary attributes need to be processed?
1782 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1783 mesh->bdr_attributes.Max() : 0);
1784 bdr_attr_marker = 0;
1785 for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
1786 {
1787 if (boundary_trace_face_integs_marker[k] == NULL)
1788 {
1789 bdr_attr_marker = 1;
1790 break;
1791 }
1793 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1794 "invalid boundary marker for boundary trace face"
1795 "integrator #" << k << ", counting from zero");
1796 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1797 {
1798 bdr_attr_marker[i] |= bdr_marker[i];
1799 }
1800 }
1801
1802 for (int i = 0; i < trial_fes -> GetNBE(); i++)
1803 {
1804 const int bdr_attr = mesh->GetBdrAttribute(i);
1805 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1806
1807 ftr = mesh->GetBdrFaceTransformations(i);
1808 if (ftr)
1809 {
1810 const int iface = mesh->GetBdrElementFaceIndex(i);
1813 trial_face_fe = trial_fes->GetFaceElement(iface);
1814 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1815 // The test_fe2 object is really a dummy and not used on the
1816 // boundaries, but we can't dereference a NULL pointer, and we don't
1817 // want to actually make a fake element.
1818 test_fe2 = test_fe1;
1819 for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
1820 {
1822 (*boundary_trace_face_integs_marker[k])[bdr_attr-1] == 0)
1823 { continue; }
1824
1825 boundary_trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe,
1826 *test_fe1,
1827 *test_fe2,
1828 *ftr, elemmat);
1830 }
1831 }
1832 }
1833 }
1834}
1835
1837 Vector &diag) const
1838{
1839 if (ext)
1840 {
1841 MFEM_ASSERT(diag.Size() == test_fes->GetTrueVSize(),
1842 "Vector for holding diagonal has wrong size!");
1843 MFEM_ASSERT(D.Size() == trial_fes->GetTrueVSize(),
1844 "Vector for holding diagonal has wrong size!");
1845 const Operator *P_trial = trial_fes->GetProlongationMatrix();
1846 const Operator *P_test = test_fes->GetProlongationMatrix();
1847 if (!IsIdentityProlongation(P_trial))
1848 {
1849 Vector local_D(P_trial->Height());
1850 P_trial->Mult(D, local_D);
1851
1852 if (!IsIdentityProlongation(P_test))
1853 {
1854 Vector local_diag(P_test->Height());
1855 ext->AssembleDiagonal_ADAt(local_D, local_diag);
1856 P_test->MultTranspose(local_diag, diag);
1857 }
1858 else
1859 {
1860 ext->AssembleDiagonal_ADAt(local_D, diag);
1861 }
1862 }
1863 else
1864 {
1865 if (!IsIdentityProlongation(P_test))
1866 {
1867 Vector local_diag(P_test->Height());
1868 ext->AssembleDiagonal_ADAt(D, local_diag);
1869 P_test->MultTranspose(local_diag, diag);
1870 }
1871 else
1872 {
1873 ext->AssembleDiagonal_ADAt(D, diag);
1874 }
1875 }
1876 }
1877 else
1878 {
1879 MFEM_ABORT("Not implemented. Maybe assemble your bilinear form into a "
1880 "matrix and use SparseMatrix functions?");
1881 }
1882}
1883
1885{
1887 {
1888 MFEM_WARNING("Conforming assemble not supported for this assembly level!");
1889 return;
1890 }
1891
1892 Finalize();
1893
1895 if (P2)
1896 {
1897 SparseMatrix *R = Transpose(*P2);
1898 SparseMatrix *RA = mfem::Mult(*R, *mat);
1899 delete R;
1900 delete mat;
1901 mat = RA;
1902 }
1903
1905 if (P1)
1906 {
1907 SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1908 delete mat;
1909 mat = RAP;
1910 }
1911
1912 height = mat->Height();
1913 width = mat->Width();
1914}
1915
1916
1918{
1919 const FiniteElement &trial_fe = *trial_fes->GetFE(i);
1920 const FiniteElement &test_fe = *test_fes->GetFE(i);
1921
1922 if (domain_integs.Size())
1923 {
1925 domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1926 elmat);
1927 for (int k = 1; k < domain_integs.Size(); k++)
1928 {
1929 domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1930 elemmat);
1931 elmat += elemmat;
1932 }
1933 }
1934 else
1935 {
1936 const int tr_dofs = trial_fe.GetDof() * trial_fes->GetVDim();
1937 const int te_dofs = test_fe.GetDof() * test_fes->GetVDim();
1938
1939 elmat.SetSize(te_dofs, tr_dofs);
1940 elmat = 0.0;
1941 }
1942}
1943
1945{
1946 const FiniteElement &trial_be = *trial_fes->GetBE(i);
1947 const FiniteElement &test_be = *test_fes->GetBE(i);
1948
1949 if (boundary_integs.Size())
1950 {
1952 boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1953 elmat);
1954 for (int k = 1; k < boundary_integs.Size(); k++)
1955 {
1956 boundary_integs[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1957 elemmat);
1958 elmat += elemmat;
1959 }
1960 }
1961 else
1962 {
1963 const int tr_dofs = trial_be.GetDof() * trial_fes->GetVDim();
1964 const int te_dofs = test_be.GetDof() * test_fes->GetVDim();
1965
1966 elmat.SetSize(te_dofs, tr_dofs);
1967 elmat = 0.0;
1968 }
1969}
1970
1972{
1974 Mesh *mesh = test_fes -> GetMesh();
1975 ftr = mesh->GetFaceElementTransformations(i);
1976 MFEM_ASSERT(ftr, "No associated face transformations.");
1977
1978 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
1979
1980 trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
1981 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1982 if (ftr->Elem2No >= 0)
1983 {
1984 trial_fe2 = trial_fes->GetFE(ftr->Elem2No);
1985 test_fe2 = test_fes->GetFE(ftr->Elem2No);
1986 }
1987 else
1988 {
1989 // The test_fe2 object is really a dummy and not used on the
1990 // boundaries, but we can't dereference a NULL pointer, and we don't
1991 // want to actually make a fake element.
1992 trial_fe2 = trial_fe1;
1993 test_fe2 = test_fe1;
1994 }
1995
1996 if (interior_face_integs.Size())
1997 {
1998 interior_face_integs[0]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
1999 *test_fe2,
2000 *ftr, elmat);
2001 for (int k = 1; k < interior_face_integs.Size(); k++)
2002 {
2003 interior_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
2004 *test_fe2,
2005 *ftr, elemmat);
2006 elmat += elemmat;
2007 }
2008 }
2009 else
2010 {
2011 int tr_dofs = trial_fe1->GetDof() * trial_fes->GetVDim();
2012 int te_dofs = test_fe1->GetDof() * test_fes->GetVDim();
2013 if (ftr->Elem2No >= 0)
2014 {
2015 tr_dofs += trial_fe2->GetDof() * trial_fes->GetVDim();
2016 te_dofs += test_fe2->GetDof() * test_fes->GetVDim();
2017 }
2018
2019 elmat.SetSize(te_dofs, tr_dofs);
2020 elmat = 0.0;
2021 }
2022}
2023
2025{
2027 Mesh *mesh = test_fes -> GetMesh();
2028 ftr = mesh->GetBdrFaceTransformations(i);
2029 MFEM_ASSERT(ftr, "No associated boundary face.");
2030
2031 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
2032
2033 trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
2034 test_fe1 = test_fes->GetFE(ftr->Elem1No);
2035 // The test_fe2 object is really a dummy and not used on the
2036 // boundaries, but we can't dereference a NULL pointer, and we don't
2037 // want to actually make a fake element.
2038 trial_fe2 = trial_fe1;
2039 test_fe2 = test_fe1;
2040
2041 if (boundary_face_integs.Size())
2042 {
2043 boundary_face_integs[0]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
2044 *test_fe2,
2045 *ftr, elmat);
2046 for (int k = 1; k < boundary_face_integs.Size(); k++)
2047 {
2048 boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
2049 *test_fe2,
2050 *ftr, elemmat);
2051 elmat += elemmat;
2052 }
2053 }
2054 else
2055 {
2056 const int tr_dofs = trial_fe1->GetDof() * trial_fes->GetVDim();
2057 const int te_dofs = test_fe1->GetDof() * test_fes->GetVDim();
2058
2059 elmat.SetSize(te_dofs, tr_dofs);
2060 elmat = 0.0;
2061 }
2062}
2063
2065{
2067 Mesh *mesh = test_fes -> GetMesh();
2068 ftr = mesh->GetFaceElementTransformations(i);
2069 MFEM_ASSERT(ftr, "No associated face transformation.");
2070
2071 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
2072
2073 trial_face_fe = trial_fes->GetFaceElement(i);
2074 test_fe1 = test_fes->GetFE(ftr->Elem1No);
2075 if (ftr->Elem2No >= 0)
2076 {
2077 test_fe2 = test_fes->GetFE(ftr->Elem2No);
2078 }
2079 else
2080 {
2081 // The test_fe2 object is really a dummy and not used on the
2082 // boundaries, but we can't dereference a NULL pointer, and we don't
2083 // want to actually make a fake element.
2084 test_fe2 = test_fe1;
2085 }
2086
2087 if (trace_face_integs.Size())
2088 {
2089 trace_face_integs[0]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
2090 *ftr, elmat);
2091 for (int k = 1; k < trace_face_integs.Size(); k++)
2092 {
2093 trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
2094 *ftr, elemmat);
2095 elmat += elemmat;
2096 }
2097 }
2098 else
2099 {
2100 const int tr_face_dofs = trial_face_fe->GetDof() * trial_fes->GetVDim();
2101 int te_dofs = test_fe1->GetDof() * test_fes->GetVDim();
2102 if (ftr->Elem2No >= 0)
2103 {
2104 te_dofs += test_fe2->GetDof() * test_fes->GetVDim();
2105 }
2106
2107 elmat.SetSize(te_dofs, tr_face_dofs);
2108 elmat = 0.0;
2109 }
2110}
2111
2113 DenseMatrix &elmat) const
2114{
2116 Mesh *mesh = test_fes -> GetMesh();
2117 ftr = mesh->GetBdrFaceTransformations(i);
2118 MFEM_ASSERT(ftr, "No associated boundary face.");
2119
2120 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
2121 int iface = mesh->GetBdrElementFaceIndex(i);
2122 trial_face_fe = trial_fes->GetFaceElement(iface);
2123 test_fe1 = test_fes->GetFE(ftr->Elem1No);
2124 // The test_fe2 object is really a dummy and not used on the
2125 // boundaries, but we can't dereference a NULL pointer, and we don't
2126 // want to actually make a fake element.
2127 test_fe2 = test_fe1;
2128
2129 if (boundary_trace_face_integs.Size())
2130 {
2131 boundary_trace_face_integs[0]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
2132 *test_fe2,
2133 *ftr, elmat);
2134 for (int k = 1; k < boundary_trace_face_integs.Size(); k++)
2135 {
2136 boundary_trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
2137 *test_fe2,
2138 *ftr, elemmat);
2139 elmat += elemmat;
2140 }
2141 }
2142 else
2143 {
2144 const int tr_face_dofs = trial_face_fe->GetDof() * trial_fes->GetVDim();
2145 int te_dofs = test_fe1->GetDof() * test_fes->GetVDim();
2146
2147 elmat.SetSize(te_dofs, tr_face_dofs);
2148 elmat = 0.0;
2149 }
2150}
2151
2153 int i, const DenseMatrix &elmat, int skip_zeros)
2154{
2155 AssembleElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
2156}
2157
2159 int i, const DenseMatrix &elmat, Array<int> &trial_vdofs_,
2160 Array<int> &test_vdofs_, int skip_zeros)
2161{
2162 trial_fes->GetElementVDofs(i, trial_vdofs_);
2163 test_fes->GetElementVDofs(i, test_vdofs_);
2164 if (mat == NULL)
2165 {
2166 mat = new SparseMatrix(height, width);
2167 }
2168 mat->AddSubMatrix(test_vdofs_, trial_vdofs_, elmat, skip_zeros);
2169}
2170
2172 int i, const DenseMatrix &elmat, int skip_zeros)
2173{
2174 AssembleBdrElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
2175}
2176
2178 int i, const DenseMatrix &elmat, Array<int> &trial_vdofs_,
2179 Array<int> &test_vdofs_, int skip_zeros)
2180{
2181 trial_fes->GetBdrElementVDofs(i, trial_vdofs_);
2182 test_fes->GetBdrElementVDofs(i, test_vdofs_);
2183 if (mat == NULL)
2184 {
2185 mat = new SparseMatrix(height, width);
2186 }
2187 mat->AddSubMatrix(test_vdofs_, trial_vdofs_, elmat, skip_zeros);
2188}
2189
2191 const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
2192{
2193 Array<int> trial_ess_dofs;
2194 trial_fes->GetEssentialVDofs(bdr_attr_is_ess, trial_ess_dofs);
2195 mat->EliminateCols(trial_ess_dofs, &sol, &rhs);
2196}
2197
2199 &bdr_attr_is_ess)
2200{
2201 Array<int> trial_ess_dofs;
2202 trial_fes->GetEssentialVDofs(bdr_attr_is_ess, trial_ess_dofs);
2203 mat->EliminateCols(trial_ess_dofs);
2204}
2205
2207 const Vector &sol, Vector &rhs)
2208{
2209 Array<int> trial_vdofs_marker;
2211 trial_vdofs_marker);
2212 mat->EliminateCols(trial_vdofs_marker, &sol, &rhs);
2213}
2214
2216{
2217 if (mat_e == NULL)
2218 {
2219 mat_e = new SparseMatrix(mat->Height(), mat->Width());
2220 }
2221
2222 Array<int> trial_vdofs_marker;
2224 trial_vdofs_marker);
2225 mat->EliminateCols(trial_vdofs_marker, *mat_e);
2226 mat_e->Finalize();
2227}
2228
2230 const Vector &x, Vector &b)
2231{
2232 mat_e->AddMult(x, b, -1.);
2233}
2234
2236 const Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
2237{
2238 mat->EliminateCols(marked_vdofs, &sol, &rhs);
2239}
2240
2242 &bdr_attr_is_ess)
2243{
2244 int i, j, k;
2245 Array<int> te_vdofs;
2246
2247 for (i = 0; i < test_fes -> GetNBE(); i++)
2248 if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
2249 {
2250 test_fes -> GetBdrElementVDofs (i, te_vdofs);
2251 for (j = 0; j < te_vdofs.Size(); j++)
2252 {
2253 if ( (k = te_vdofs[j]) < 0 )
2254 {
2255 k = -1-k;
2256 }
2257 mat -> EliminateRow (k);
2258 }
2259 }
2260}
2261
2263{
2264 for (int i=0; i<test_vdofs_.Size(); ++i)
2265 {
2266 mat->EliminateRow(test_vdofs_[i]);
2267 }
2268}
2269
2271 const Array<int> &trial_tdof_list,
2272 const Array<int> &test_tdof_list,
2273 OperatorHandle &A)
2274
2275{
2276 if (ext)
2277 {
2278 ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
2279 return;
2280 }
2281
2284
2285 mat->Finalize();
2286
2287 if (test_P && trial_P)
2288 {
2289 SparseMatrix *m = RAP(*test_P, *mat, *trial_P);
2290 delete mat;
2291 mat = m;
2292 }
2293 else if (test_P)
2294 {
2295 SparseMatrix *m = TransposeMult(*test_P, *mat);
2296 delete mat;
2297 mat = m;
2298 }
2299 else if (trial_P)
2300 {
2301 SparseMatrix *m = mfem::Mult(*mat, *trial_P);
2302 delete mat;
2303 mat = m;
2304 }
2305
2306 EliminateTrialVDofs(trial_tdof_list);
2307 EliminateTestVDofs(test_tdof_list);
2308
2309 A.Reset(mat, false);
2310}
2311
2313 const Array<int> &trial_tdof_list,
2314 const Array<int> &test_tdof_list,
2315 Vector &x, Vector &b,
2316 OperatorHandle &A,
2317 Vector &X, Vector &B)
2318{
2319 if (ext)
2320 {
2321 ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
2322 x, b, A, X, B);
2323 return;
2324 }
2325
2326 const Operator *Pi = this->GetProlongation();
2327 const Operator *Po = this->GetOutputProlongation();
2328 const Operator *Ri = this->GetRestriction();
2329 InitTVectors(Po, Ri, Pi, x, b, X, B);
2330
2331 if (!mat_e)
2332 {
2333 FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list,
2334 A); // Set A = mat_e
2335 }
2336 // Eliminate essential BCs with B -= Ab xb
2337 EliminateTrialVDofsInRHS(trial_tdof_list, X, B);
2338
2339 B.SetSubVector(test_tdof_list, 0.0);
2340}
2341
2343{
2344 delete mat;
2345 mat = NULL;
2346 delete mat_e;
2347 mat_e = NULL;
2350 if (ext) { ext->Update(); }
2351}
2352
2354{
2355 if (mat) { delete mat; }
2356 if (mat_e) { delete mat_e; }
2357 if (!extern_bfs)
2358 {
2359 int i;
2360 for (i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
2361 for (i = 0; i < boundary_integs.Size(); i++)
2362 { delete boundary_integs[i]; }
2363 for (i = 0; i < interior_face_integs.Size(); i++)
2364 { delete interior_face_integs[i]; }
2365 for (i = 0; i < boundary_face_integs.Size(); i++)
2366 { delete boundary_face_integs[i]; }
2367 for (i = 0; i < trace_face_integs.Size(); i++)
2368 { delete trace_face_integs[i]; }
2369 for (i = 0; i < boundary_trace_face_integs.Size(); i++)
2370 { delete boundary_trace_face_integs[i]; }
2371 }
2372}
2373
2375{
2376 if (ext)
2377 {
2378 MFEM_ABORT("the assembly level has already been set!");
2379 }
2380 assembly = assembly_level;
2381 switch (assembly)
2382 {
2385 // Use the original implementation for now
2386 break;
2388 MFEM_ABORT("Element assembly not supported yet... stay tuned!");
2389 break;
2391 ext.reset(new PADiscreteLinearOperatorExtension(this));
2392 break;
2394 MFEM_ABORT("Matrix-free action not supported yet... stay tuned!");
2395 break;
2396 default:
2397 MFEM_ABORT("Unknown assembly level");
2398 }
2399}
2400
2402{
2403 if (ext)
2404 {
2405 ext->Assemble();
2406 return;
2407 }
2408
2409 ElementTransformation *eltrans;
2410 DofTransformation * dom_dof_trans;
2411 DofTransformation * ran_dof_trans;
2412 DenseMatrix elmat;
2413
2414 Mesh *mesh = test_fes->GetMesh();
2415
2416 if (mat == NULL)
2417 {
2418 mat = new SparseMatrix(height, width);
2419 }
2420
2421 if (domain_integs.Size())
2422 {
2423 for (int k = 0; k < domain_integs.Size(); k++)
2424 {
2425 if (domain_integs_marker[k] != NULL)
2426 {
2427 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
2428 (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
2429 "invalid element marker for domain integrator #"
2430 << k << ", counting from zero");
2431 }
2432 }
2433
2434 for (int i = 0; i < test_fes->GetNE(); i++)
2435 {
2436 const int elem_attr = mesh->GetAttribute(i);
2437 dom_dof_trans = trial_fes->GetElementVDofs(i, trial_vdofs);
2438 ran_dof_trans = test_fes->GetElementVDofs(i, test_vdofs);
2439 eltrans = test_fes->GetElementTransformation(i);
2440
2442 elmat = 0.0;
2443 for (int k = 0; k < domain_integs.Size(); k++)
2444 {
2445 if (domain_integs_marker[k] == NULL ||
2446 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
2447 {
2448 domain_integs[k]->AssembleElementMatrix2(*trial_fes->GetFE(i),
2449 *test_fes->GetFE(i),
2450 *eltrans, elemmat);
2451 elmat += elemmat;
2452 }
2453 }
2454 if (ran_dof_trans || dom_dof_trans)
2455 {
2456 TransformPrimal(ran_dof_trans, dom_dof_trans, elemmat);
2457 }
2459 }
2460 }
2461
2462 if (trace_face_integs.Size())
2463 {
2464 const int nfaces = test_fes->GetMesh()->GetNumFaces();
2465 for (int i = 0; i < nfaces; i++)
2466 {
2469 eltrans = test_fes->GetMesh()->GetFaceTransformation(i);
2470
2472 elmat = 0.0;
2473 for (int k = 0; k < trace_face_integs.Size(); k++)
2474 {
2475 trace_face_integs[k]->AssembleElementMatrix2(*trial_fes->GetFaceElement(i),
2477 *eltrans, elemmat);
2478 elmat += elemmat;
2479 }
2480 mat->SetSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
2481 }
2482 }
2483}
2484
2485}
Dynamic 2D array using row-major layout.
Definition array.hpp:392
void SetSize(int m, int n)
Definition array.hpp:405
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:341
int Size() const
Return the logical size of the array.
Definition array.hpp:147
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
Abstract base class BilinearFormIntegrator.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void AssembleDiagonal(Vector &diag) const override
Assemble the diagonal of the bilinear form into diag. Note that diag is a tdof Vector.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
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.
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
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> &,...
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 AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const override
Add the matrix transpose vector multiplication: .
void ComputeFaceMatrix(int i, DenseMatrix &elmat) const
Compute the face matrix of the given face element.
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.
void Finalize(int skip_zeros=1) override
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
virtual ~BilinearForm()
Deletes internal matrices, bilinear integrators, and the BilinearFormExtension.
BilinearForm()
may be used in the construction of derived classes
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 ComputeBdrFaceMatrix(int i, DenseMatrix &elmat) const
Compute the boundary face matrix of the given boundary element.
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.
void EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
Enable hybridization.
real_t & Elem(int i, int j) override
Returns a reference to: .
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.
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.
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,...
void MultTranspose(const Vector &x, Vector &y) const override
Matrix transpose vector multiplication: .
std::unique_ptr< StaticCondensation > static_cond
const DenseTensor & GetElementMatrices()
Return a DenseTensor containing the assembled element matrices.
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x) override
Recover the solution of a linear system formed with FormLinearSystem().
std::unique_ptr< DenseTensor > element_matrices
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(....
void ComputeElementMatrix(int i, DenseMatrix &elmat) const
Compute the element matrix of the given element.
int Size() const
Get the size of the BilinearForm as a square matrix.
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication: .
std::unique_ptr< BilinearFormExtension > ext
Extension for supporting Full Assembly (FA), Element Assembly (EA),Partial Assembly (PA),...
MatrixInverse * Inverse() const override
Returns a pointer to (approximation) of the matrix inverse: (currently returns NULL)
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.
std::unique_ptr< Hybridization > hybridization
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, real_t value)
Perform elimination and set the diagonal entry to the given value.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat) const
Compute the boundary element matrix of the given boundary element.
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.
Square Operator for imposing essential boundary conditions using only the action, Mult(),...
Definition operator.hpp:992
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate "essential boundary condition" values specified in x from the given right-hand side b.
Definition operator.cpp:559
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:118
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
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)
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:750
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1463
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:1287
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition fespace.cpp:3881
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:851
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:924
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:750
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:809
const NURBSExtension * GetNURBSext() const
Definition fespace.hpp:681
virtual const Operator * GetProlongationMatrix() const
Definition fespace.hpp:727
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:332
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:3835
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:876
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition fespace.hpp:933
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:891
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:361
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
Definition fespace.cpp:3914
const SparseMatrix * GetConformingProlongation() const
Definition fespace.cpp:1456
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:848
int GetVDim() const
Returns the vector dimension of the finite element space.
Definition fespace.hpp:841
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3871
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:348
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:584
void GetPatchVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom in vdofs for NURBS patch i.
Definition fespace.cpp:355
Abstract class for all finite elements.
Definition fe_base.hpp:244
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
Auxiliary class Hybridization, used to implement BilinearForm hybridization.
A standard isoparametric element transformation.
Definition eltrans.hpp:629
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:64
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition mesh.cpp:1004
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:290
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:298
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6531
int GetAttribute(int i) const
Return the attribute of element i.
Definition mesh.hpp:1389
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition mesh.hpp:1395
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
Definition mesh.hpp:1588
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.
Definition mesh.cpp:607
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition mesh.cpp:1123
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Definition mesh.cpp:1103
Table * GetFaceToElementTable() const
Definition mesh.cpp:7480
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:288
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds a boundary face integrator. Assumes ownership of bfi.
void ConformingAssemble()
For partially conforming trial and/or test FE spaces, complete the assembly process by performing wh...
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const override
Add the matrix vector multiple to a vector: .
void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const override
Add the matrix transpose vector multiplication: .
Array< BilinearFormIntegrator * > boundary_face_integs
Boundary face integrators.
void EliminateTrialVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
Use the stored eliminated part of the matrix (see EliminateTrialVDofs(const Array<int> &)) to modify ...
void Assemble(int skip_zeros=1)
real_t & Elem(int i, int j) override
Returns a reference to: .
void Finalize(int skip_zeros=1) override
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds a boundary integrator. Assumes ownership of bfi.
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds an interior face integrator. Assumes ownership of bfi.
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 ComputeFaceMatrix(int i, DenseMatrix &elmat) const
Compute the face matrix of the given face element.
void ComputeBdrTraceFaceMatrix(int i, DenseMatrix &elmat) const
Compute the boundary trace face matrix of the given boundary element.
void MultTranspose(const Vector &x, Vector &y) const override
Matrix transpose vector multiplication: .
void ComputeElementMatrix(int i, DenseMatrix &elmat) const
Compute the element matrix of the given element.
Array< Array< int > * > boundary_face_integs_marker
Entries are not owned.
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::...
MatrixInverse * Inverse() const override
Returns a pointer to (approximation) of the matrix inverse: (currently unimplemented and returns NUL...
void EliminateTrialVDofs(const Array< int > &vdofs, const Vector &sol, Vector &rhs)
Eliminate the given trial vdofs. NOTE: here, vdofs is a list of DOFs.
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.
const Operator * GetOutputProlongation() const override
Get the test finite element space prolongation matrix.
Array< BilinearFormIntegrator * > interior_face_integs
Interior face integrators.
void ComputeTraceFaceMatrix(int i, DenseMatrix &elmat) const
Compute the trace face matrix of the given face element.
void EliminateTrialEssentialBC(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
Eliminate essential boundary trial DOFs from the system.
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.
Array< BilinearFormIntegrator * > boundary_trace_face_integs
Boundary trace face (skeleton) integrators.
const Operator * GetRestriction() const override
Get the input finite element space restriction matrix.
Array< BilinearFormIntegrator * > domain_integs
Domain integrators.
SparseMatrix * mat
Owned.
void Update()
Must be called after making changes to trial_fes or test_fes.
Array< Array< int > * > domain_integs_marker
Entries are not owned.
const Operator * GetProlongation() const override
Get the input finite element space prolongation matrix.
void ComputeBdrFaceMatrix(int i, DenseMatrix &elmat) const
Compute the boundary face matrix of the given boundary element.
Array< BilinearFormIntegrator * > boundary_integs
Boundary integrators.
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Similar to EliminateTrialVDofs(const Array<int> &, const Vector &, Vector &) but here ess_dofs is a m...
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(....
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...
std::unique_ptr< MixedBilinearFormExtension > ext
void EliminateTestEssentialBC(const Array< int > &bdr_attr_is_ess)
Eliminate essential boundary test DOFs from the system matrix.
void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi)
Add a trace face integrator. Assumes ownership of bfi.
void EliminateTestVDofs(const Array< int > &vdofs)
Eliminate the given test vdofs. NOTE: here, vdofs is a list of DOFs.
void Mult(const Vector &x, Vector &y) const override
Matrix multiplication: .
virtual ~MixedBilinearForm()
Deletes internal matrices, bilinear integrators, and the BilinearFormExtension.
FiniteElementSpace * test_fes
Not owned.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
AssemblyLevel assembly
The form assembly level (full, partial, etc.)
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat) const
Compute the boundary element matrix of the given boundary element.
int GetNP() const
Return the number of patches.
Definition nurbs.hpp:745
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
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 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
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
int GetRow(const int row, Array< int > &cols, Vector &srow) const override
Extract all column indices and values from a given row.
void GetDiag(Vector &d) const
Returns the Diagonal of A.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
void MultTranspose(const Vector &x, Vector &y) const override
Multiply a vector with the transposed matrix. y = At * x.
bool Finalized() const
Returns whether or not CSR format has been finalized.
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).
void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const override
y += At * x (default) or y += a * At * x
void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const override
y += A * x (default) or y += a * A * x
void Mult(const Vector &x, Vector &y) const override
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
void Finalize(int skip_zeros=1) override
Finalize the matrix initialization, switching the storage format from LIL to CSR.
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 LoseData()
Call this if data has been stolen.
Definition table.hpp:188
int * GetJ()
Definition table.hpp:122
int * GetI()
Definition table.hpp:121
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition table.cpp:271
Vector data type.
Definition vector.hpp:82
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:679
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition vector.hpp:264
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
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:814
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:622
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:548
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:486
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:817
float real_t
Definition config.hpp:43
@ NATIVE
Native ordering as defined by the FiniteElement.
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)