MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
weakform.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#include "weakform.hpp"
13
14namespace mfem
15{
16
18{
19 trial_integs.SetSize(trial_fes.Size(), test_fecols.Size());
20 for (int i = 0; i < trial_integs.NumRows(); i++)
21 {
22 for (int j = 0; j < trial_integs.NumCols(); j++)
23 {
25 }
26 }
27
28 test_integs.SetSize(test_fecols.Size(), test_fecols.Size());
29 for (int i = 0; i < test_integs.NumRows(); i++)
30 {
31 for (int j = 0; j < test_integs.NumCols(); j++)
32 {
34 }
35 }
36
37 lfis.SetSize(test_fecols.Size());
38 for (int j = 0; j < lfis.Size(); j++)
39 {
41 }
42
43
45
46 mat = mat_e = NULL;
49 width = height;
50
51 initialized = true;
52 static_cond = nullptr;
53
55 {
56 Bmat.SetSize(mesh->GetNE());
57 fvec.SetSize(mesh->GetNE());
58 }
59}
60
62{
65 dof_offsets[0] = 0;
66 tdof_offsets[0] = 0;
67 for (int i =0; i<nblocks; i++)
68 {
69 dof_offsets[i+1] = trial_fes[i]->GetVSize();
70 tdof_offsets[i+1] = trial_fes[i]->GetTrueVSize();
71 }
74}
75
76// Allocate SparseMatrix and RHS
78{
79 if (static_cond) { return; }
80
82 mat->owns_blocks = 1;
83
84 for (int i = 0; i<mat->NumRowBlocks(); i++)
85 {
86 int h = dof_offsets[i+1] - dof_offsets[i];
87 for (int j = 0; j<mat->NumColBlocks(); j++)
88 {
89 int w = dof_offsets[j+1] - dof_offsets[j];
90 mat->SetBlock(i,j,new SparseMatrix(h, w));
91 }
92 }
94 *y = 0.;
95}
96
97void DPGWeakForm::Finalize(int skip_zeros)
98{
99 if (mat) { mat->Finalize(skip_zeros); }
100 if (mat_e) { mat_e->Finalize(skip_zeros); }
101 if (static_cond) { static_cond->Finalize(); }
102}
103
104/// Adds new Domain BF Integrator. Assumes ownership of @a bfi.
106 BilinearFormIntegrator *bfi, int n, int m)
107{
108 MFEM_VERIFY(n>=0 && n<trial_fes.Size(),
109 "DPGWeakFrom::AddTrialIntegrator: trial fespace index out of bounds");
110 MFEM_VERIFY(m>=0 && m<test_fecols.Size(),
111 "DPGWeakFrom::AddTrialIntegrator: test fecol index out of bounds");
112 trial_integs(n,m)->Append(bfi);
113}
114
115/// Adds new Domain BF Integrator. Assumes ownership of @a bfi.
117(BilinearFormIntegrator *bfi, int n, int m)
118{
119 MFEM_VERIFY(n>=0 && n<test_fecols.Size() && m>=0 && m<test_fecols.Size(),
120 "DPGWeakFrom::AdTestIntegrator: test fecol index out of bounds");
121 test_integs(n,m)->Append(bfi);
122}
123
124/// Adds new Domain LF Integrator. Assumes ownership of @a bfi.
126 LinearFormIntegrator *lfi, int n)
127{
128 MFEM_VERIFY(n>=0 && n<test_fecols.Size(),
129 "DPGWeakFrom::AddDomainLFIntegrator: test fecol index out of bounds");
130 lfis[n]->Append(lfi);
131}
132
134{
137 P->owns_blocks = 0;
138 R->owns_blocks = 0;
139 for (int i = 0; i<nblocks; i++)
140 {
141 const SparseMatrix *P_ = trial_fes[i]->GetConformingProlongation();
142 if (P_)
143 {
144 const SparseMatrix *R_ = trial_fes[i]->GetRestrictionMatrix();
145 P->SetBlock(i,i,const_cast<SparseMatrix*>(P_));
146 R->SetBlock(i,i,const_cast<SparseMatrix*>(R_));
147 }
148 }
149}
150
152{
153 Finalize(0);
154 if (!P) { BuildProlongation(); }
155
156 BlockMatrix * Pt = Transpose(*P);
157 BlockMatrix * PtA = mfem::Mult(*Pt, *mat);
158 mat->owns_blocks = 0;
159 for (int i = 0; i<nblocks; i++)
160 {
161 for (int j = 0; j<nblocks; j++)
162 {
163 SparseMatrix * tmp = &mat->GetBlock(i,j);
164 if (Pt->IsZeroBlock(i,i))
165 {
166 PtA->SetBlock(i,j,tmp);
167 }
168 else
169 {
170 delete tmp;
171 }
172 }
173 }
174 delete mat;
175 if (mat_e)
176 {
177 BlockMatrix *PtAe = mfem::Mult(*Pt, *mat_e);
178 mat_e->owns_blocks = 0;
179 for (int i = 0; i<nblocks; i++)
180 {
181 for (int j = 0; j<nblocks; j++)
182 {
183 SparseMatrix * tmp = &mat_e->GetBlock(i,j);
184 if (Pt->IsZeroBlock(i,i))
185 {
186 PtAe->SetBlock(i,j,tmp);
187 }
188 else
189 {
190 delete tmp;
191 }
192 }
193 }
194 delete mat_e;
195 mat_e = PtAe;
196 }
197 delete Pt;
198
199 mat = mfem::Mult(*PtA, *P);
200
201 PtA->owns_blocks = 0;
202 for (int i = 0; i<nblocks; i++)
203 {
204 for (int j = 0; j<nblocks; j++)
205 {
206 SparseMatrix * tmp = &PtA->GetBlock(j,i);
207 if (P->IsZeroBlock(i,i))
208 {
209 mat->SetBlock(j,i,tmp);
210 }
211 else
212 {
213 delete tmp;
214 }
215 }
216 }
217 delete PtA;
218
219 if (mat_e)
220 {
221 BlockMatrix *PtAeP = mfem::Mult(*mat_e, *P);
222 mat_e->owns_blocks = 0;
223 for (int i = 0; i<nblocks; i++)
224 {
225 for (int j = 0; j<nblocks; j++)
226 {
227 SparseMatrix * tmp = &mat_e->GetBlock(j,i);
228 if (P->IsZeroBlock(i,i))
229 {
230 PtAeP->SetBlock(j,i,tmp);
231 }
232 else
233 {
234 delete tmp;
235 }
236 }
237 }
238
239 delete mat_e;
240 mat_e = PtAeP;
241 }
242 height = mat->Height();
243 width = mat->Width();
244}
245
246/// Assembles the form i.e. sums over all domain integrators.
247void DPGWeakForm::Assemble(int skip_zeros)
248{
249 ElementTransformation *eltrans;
250 Array<int> faces, ori;
251
252 DofTransformation * doftrans_i, *doftrans_j;
253 if (mat == NULL)
254 {
255 AllocMat();
256 }
257
258 // loop through the elements
259 int dim = mesh->Dimension();
260 DenseMatrix B, Be, G, Ge, A;
261 Vector vec_e, vec, Gvec, b;
262 Array<int> vdofs;
263
264 // loop through elements
265 for (int iel = 0; iel < mesh -> GetNE(); iel++)
266 {
267 if (dim == 1)
268 {
269 mesh->GetElementVertices(iel, faces);
270 }
271 else if (dim == 2)
272 {
273 mesh->GetElementEdges(iel, faces, ori);
274 }
275 else if (dim == 3)
276 {
277 mesh->GetElementFaces(iel,faces,ori);
278 }
279 else
280 {
281 MFEM_ABORT("DPGWeakForm::Assemble: dim > 3 not supported");
282 }
283 int numfaces = faces.Size();
284
285 Array<int> test_offs(test_fecols.Size()+1); test_offs[0] = 0;
286 Array<int> trial_offs(trial_fes.Size()+1); trial_offs = 0;
287
288 eltrans = mesh->GetElementTransformation(iel);
289 for (int j = 0; j < test_fecols.Size(); j++)
290 {
291 int order = test_fecols[j]->GetOrder(); // assuming uniform order
292 test_offs[j+1] = test_fecols_vdims[j]*test_fecols[j]->GetFE(
293 eltrans->GetGeometryType(),
294 order)->GetDof();
295 }
296 for (int j = 0; j < trial_fes.Size(); j++)
297 {
298 if (IsTraceFes[j])
299 {
300 for (int ie = 0; ie<faces.Size(); ie++)
301 {
302 trial_offs[j+1] += trial_fes[j]->GetVDim()*trial_fes[j]->GetFaceElement(
303 faces[ie])->GetDof();
304 }
305 }
306 else
307 {
308 trial_offs[j+1] = trial_fes[j]->GetVDim() * trial_fes[j]->GetFE(
309 iel)->GetDof();
310 }
311 }
312 test_offs.PartialSum();
313 trial_offs.PartialSum();
314
315 G.SetSize(test_offs.Last()); G = 0.0;
316 vec.SetSize(test_offs.Last()); vec = 0.0;
317 B.SetSize(test_offs.Last(),trial_offs.Last()); B = 0.0;
318
319
320 for (int j = 0; j < test_fecols.Size(); j++)
321 {
322 int order_j = test_fecols[j]->GetOrder();
323
324 eltrans = mesh->GetElementTransformation(iel);
325 const FiniteElement & test_fe =
326 *test_fecols[j]->GetFE(eltrans->GetGeometryType(), order_j);
327
328 for (int k = 0; k < lfis[j]->Size(); k++)
329 {
330 (*lfis[j])[k]->AssembleRHSElementVect(test_fe,*eltrans,vec_e);
331 vec.AddSubVector(vec_e,test_offs[j]);
332 }
333
334 for (int i = 0; i < test_fecols.Size(); i++)
335 {
336 int order_i = test_fecols[i]->GetOrder();
337 eltrans = mesh->GetElementTransformation(iel);
338 const FiniteElement & test_fe_i =
339 *test_fecols[i]->GetFE(eltrans->GetGeometryType(), order_i);
340
341 for (int k = 0; k < test_integs(i,j)->Size(); k++)
342 {
343 if (i==j)
344 {
345 (*test_integs(i,j))[k]->AssembleElementMatrix(test_fe,*eltrans,Ge);
346 }
347 else
348 {
349 (*test_integs(i,j))[k]->AssembleElementMatrix2(test_fe_i,test_fe,*eltrans,
350 Ge);
351 }
352 G.AddSubMatrix(test_offs[j], test_offs[i], Ge);
353 }
354 }
355
356 for (int i = 0; i < trial_fes.Size(); i++)
357 {
358 if (IsTraceFes[i])
359 {
360 for (int k = 0; k < trial_integs(i,j)->Size(); k++)
361 {
362 int face_dof_offs = 0;
363 for (int ie = 0; ie < numfaces; ie++)
364 {
365 int iface = faces[ie];
367 const FiniteElement & tfe = *trial_fes[i]->GetFaceElement(iface);
368 (*trial_integs(i,j))[k]->AssembleTraceFaceMatrix(iel,tfe,test_fe,*ftr,Be);
369 B.AddSubMatrix(test_offs[j], trial_offs[i]+face_dof_offs, Be);
370 face_dof_offs+=Be.Width();
371 }
372 }
373 }
374 else
375 {
376 const FiniteElement & fe = *trial_fes[i]->GetFE(iel);
377 eltrans = mesh->GetElementTransformation(iel);
378 for (int k = 0; k < trial_integs(i,j)->Size(); k++)
379 {
380 (*trial_integs(i,j))[k]->AssembleElementMatrix2(fe,test_fe,*eltrans,Be);
381 B.AddSubMatrix(test_offs[j], trial_offs[i], Be);
382 }
383 }
384 }
385 }
386
387 // Form Normal Equations B^T G^-1 B = B^T G^-1 l
388 Gvec.SetSize(G.Height());
389 b.SetSize(B.Width());
390 A.SetSize(B.Width());
391
392 CholeskyFactors chol(G.GetData());
393 chol.Factor(G.Height());
394
395 chol.LSolve(B.Height(), B.Width(), B.GetData());
396 chol.LSolve(vec.Size(), 1, vec.GetData());
397 if (store_matrices)
398 {
399 Bmat[iel] = new DenseMatrix(B);
400 fvec[iel] = new Vector(vec);
401 }
402 mfem::MultAtB(B,B,A);
403 B.MultTranspose(vec,b);
404
405 if (static_cond)
406 {
408 }
409 else
410 {
411 // Assembly
412 for (int i = 0; i<trial_fes.Size(); i++)
413 {
414 Array<int> vdofs_i;
415 doftrans_i = nullptr;
416 if (IsTraceFes[i])
417 {
418 Array<int> face_vdofs;
419 for (int k = 0; k < numfaces; k++)
420 {
421 int iface = faces[k];
422 trial_fes[i]->GetFaceVDofs(iface, face_vdofs);
423 vdofs_i.Append(face_vdofs);
424 }
425 }
426 else
427 {
428 doftrans_i = trial_fes[i]->GetElementVDofs(iel, vdofs_i);
429 }
430 for (int j = 0; j<trial_fes.Size(); j++)
431 {
432 Array<int> vdofs_j;
433 doftrans_j = nullptr;
434
435 if (IsTraceFes[j])
436 {
437 Array<int> face_vdofs;
438 for (int k = 0; k < numfaces; k++)
439 {
440 int iface = faces[k];
441 trial_fes[j]->GetFaceVDofs(iface, face_vdofs);
442 vdofs_j.Append(face_vdofs);
443 }
444 }
445 else
446 {
447 doftrans_j = trial_fes[j]->GetElementVDofs(iel, vdofs_j);
448 }
449
450 DenseMatrix Ae;
451 A.GetSubMatrix(trial_offs[i],trial_offs[i+1],
452 trial_offs[j],trial_offs[j+1], Ae);
453 if (doftrans_i || doftrans_j)
454 {
455 TransformDual(doftrans_i, doftrans_j, Ae);
456 }
457 mat->GetBlock(i,j).AddSubMatrix(vdofs_i,vdofs_j, Ae);
458 }
459
460 // assemble rhs
461 real_t * data = b.GetData();
462 Vector vec1;
463 // ref subvector
464 vec1.SetDataAndSize(&data[trial_offs[i]],
465 trial_offs[i+1]-trial_offs[i]);
466 if (doftrans_i)
467 {
468 doftrans_i->TransformDual(vec1);
469 }
470 y->GetBlock(i).AddElementVector(vdofs_i,vec1);
471 }
472 }
473 }
474}
475
477 &ess_tdof_list,
478 Vector &x,
479 OperatorHandle &A, Vector &X,
480 Vector &B, int copy_interior)
481{
482 FormSystemMatrix(ess_tdof_list, A);
483 if (static_cond)
484 {
485 // Schur complement reduction to the exposed dofs
486 static_cond->ReduceSystem(x, X, B, copy_interior);
487 }
488 else if (!P)
489 {
490 EliminateVDofsInRHS(ess_tdof_list, x, *y);
491 X.MakeRef(x, 0, x.Size());
492 B.MakeRef(*y, 0, y->Size());
493 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
494 }
495 else // non conforming space
496 {
497 B.SetSize(P->Width());
498
499 P->MultTranspose(*y, B);
500 real_t *data = y->GetData();
501 Vector tmp;
502 for (int i = 0; i<nblocks; i++)
503 {
504 if (P->IsZeroBlock(i,i))
505 {
506 int offset = tdof_offsets[i];
507 tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
508 B.SetVector(tmp,offset);
509 }
510 }
511
512 X.SetSize(R->Height());
513
514 R->Mult(x, X);
515 data = x.GetData();
516 for (int i = 0; i<nblocks; i++)
517 {
518 if (R->IsZeroBlock(i,i))
519 {
520 int offset = tdof_offsets[i];
521 tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
522 X.SetVector(tmp,offset);
523 }
524 }
525
526 EliminateVDofsInRHS(ess_tdof_list, X, B);
527 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
528 }
529}
530
532 &ess_tdof_list,
534{
535 if (static_cond)
536 {
538 {
539 static_cond->SetEssentialTrueDofs(ess_tdof_list);
541 }
542 A.Reset(&static_cond->GetSchurMatrix(), false);
543 }
544 else
545 {
546 if (!mat_e)
547 {
548 bool conforming = true;
549 for (int i = 0; i<nblocks; i++)
550 {
551 const SparseMatrix *P_ = trial_fes[i]->GetConformingProlongation();
552 if (P_)
553 {
554 conforming = false;
555 break;
556 }
557 }
558 if (!conforming) { ConformingAssemble(); }
559 const int remove_zeros = 0;
560 EliminateVDofs(ess_tdof_list, diag_policy);
561 Finalize(remove_zeros);
562 }
563 A.Reset(mat, false);
564 }
565}
566
568 const Array<int> &vdofs, const Vector &x, Vector &b)
569{
570 mat_e->AddMult(x,b,-1.);
571 mat->PartMult(vdofs,x,b);
572}
573
576{
577 if (mat_e == NULL)
578 {
579 Array<int> offsets;
580
581 offsets.MakeRef( (P) ? tdof_offsets : dof_offsets);
582
583 mat_e = new BlockMatrix(offsets);
584 mat_e->owns_blocks = 1;
585 for (int i = 0; i<mat_e->NumRowBlocks(); i++)
586 {
587 int h = offsets[i+1] - offsets[i];
588 for (int j = 0; j<mat_e->NumColBlocks(); j++)
589 {
590 int w = offsets[j+1] - offsets[j];
591 mat_e->SetBlock(i,j,new SparseMatrix(h, w));
592 }
593 }
594 }
596}
597
599 Vector &x)
600{
601
602 if (static_cond)
603 {
604 // Private dofs back solve
606 }
607 else if (!P)
608 {
609 x.SyncMemory(X);
610 }
611 else
612 {
613 x.SetSize(P->Height());
614 P->Mult(X, x);
615 real_t *data = X.GetData();
616 Vector tmp;
617 for (int i = 0; i<nblocks; i++)
618 {
619 if (P->IsZeroBlock(i,i))
620 {
621 int offset = tdof_offsets[i];
622 tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
623 x.SetVector(tmp,offset);
624 }
625 }
626 }
627}
628
630{
631 if (initialized)
632 {
633 for (int k = 0; k< trial_integs.NumRows(); k++)
634 {
635 for (int l = 0; l<trial_integs.NumCols(); l++)
636 {
637 for (int i = 0; i<trial_integs(k,l)->Size(); i++)
638 {
639 delete (*trial_integs(k,l))[i];
640 }
641 delete trial_integs(k,l);
642 }
643 }
644 trial_integs.DeleteAll();
645
646 for (int k = 0; k < test_integs.NumRows(); k++)
647 {
648 for (int l = 0; l < test_integs.NumCols(); l++)
649 {
650 for (int i = 0; i < test_integs(k,l)->Size(); i++)
651 {
652 delete (*test_integs(k,l))[i];
653 }
654 delete test_integs(k,l);
655 }
656 }
657 test_integs.DeleteAll();
658
659 for (int k = 0; k < lfis.Size(); k++)
660 {
661 for (int i = 0; i < lfis[k]->Size(); i++)
662 {
663 delete (*lfis[k])[i];
664 }
665 delete lfis[k];
666 }
667 lfis.DeleteAll();
668 }
669}
670
672{
673 delete mat_e; mat_e = nullptr;
674 delete mat; mat = nullptr;
675 delete y; y = nullptr;
676
677 if (P)
678 {
679 delete P; P = nullptr;
680 delete R; R = nullptr;
681 }
682
683 if (static_cond)
684 {
686 }
687 else
688 {
689 delete static_cond; static_cond = nullptr;
690 }
691
693
696 width = height;
697
698 initialized = true;
699
700 if (store_matrices)
701 {
702 for (int i = 0; i<Bmat.Size(); i++)
703 {
704 delete Bmat[i]; Bmat[i] = nullptr;
705 delete fvec[i]; fvec[i] = nullptr;
706 }
708 fvec.SetSize(mesh->GetNE());
709 for (int i = 0; i<Bmat.Size(); i++)
710 {
711 Bmat[i] = nullptr;
712 fvec[i] = nullptr;
713 }
714 }
715}
716
722
724{
725 // Element vector of trial space size
726 Vector u;
727 Array<int> vdofs;
728 Array<int> faces, ori;
729 int dim = mesh->Dimension();
731 // loop through elements
732 for (int iel = 0; iel < mesh -> GetNE(); iel++)
733 {
734 if (dim == 1)
735 {
736 mesh->GetElementVertices(iel, faces);
737 }
738 else if (dim == 2)
739 {
740 mesh->GetElementEdges(iel, faces, ori);
741 }
742 else if (dim == 3)
743 {
744 mesh->GetElementFaces(iel,faces,ori);
745 }
746 else
747 {
748 MFEM_ABORT("DPGWeakForm::ComputeResidual: "
749 "dim > 3 not supported");
750 }
751 int numfaces = faces.Size();
752
753 Array<int> trial_offs(trial_fes.Size()+1); trial_offs = 0;
754
755 for (int j = 0; j < trial_fes.Size(); j++)
756 {
757 if (IsTraceFes[j])
758 {
759 for (int ie = 0; ie<faces.Size(); ie++)
760 {
761 trial_offs[j+1] += trial_fes[j]->GetFaceElement(faces[ie])->GetDof();
762 }
763 }
764 else
765 {
766 trial_offs[j+1] = trial_fes[j]->GetVDim() * trial_fes[j]->GetFE(
767 iel)->GetDof();
768 }
769 }
770 trial_offs.PartialSum();
771
772 u.SetSize(trial_offs.Last());
773 real_t * data = u.GetData();
774 DofTransformation * doftrans = nullptr;
775 for (int i = 0; i<trial_fes.Size(); i++)
776 {
777 vdofs.SetSize(0);
778 doftrans = nullptr;
779 if (IsTraceFes[i])
780 {
781 Array<int> face_vdofs;
782 for (int k = 0; k < numfaces; k++)
783 {
784 int iface = faces[k];
785 trial_fes[i]->GetFaceVDofs(iface, face_vdofs);
786 vdofs.Append(face_vdofs);
787 }
788 }
789 else
790 {
791 doftrans = trial_fes[i]->GetElementVDofs(iel, vdofs);
792 }
793 Vector vec1;
794 vec1.SetDataAndSize(&data[trial_offs[i]],
795 trial_offs[i+1]-trial_offs[i]);
796 x.GetBlock(i).GetSubVector(vdofs,vec1);
797 if (doftrans)
798 {
799 doftrans->InvTransformPrimal(vec1);
800 }
801 } // end of loop through trial spaces
802
803 Vector v(Bmat[iel]->Height());
804 Bmat[iel]->Mult(u,v);
805 v -= *fvec[iel];
806 residuals[iel] = v.Norml2();
807 } // end of loop through elements
808 return residuals;
809}
810
812{
813 delete mat_e; mat_e = nullptr;
814 delete mat; mat = nullptr;
815 delete y; y = nullptr;
816
818
819 if (P)
820 {
821 delete P;
822 delete R;
823 }
824
825 delete static_cond;
826
827 if (store_matrices)
828 {
829 for (int i = 0; i<mesh->GetNE(); i++)
830 {
831 delete Bmat[i]; Bmat[i] = nullptr;
832 delete fvec[i]; fvec[i] = nullptr;
833 }
834 }
835}
836
837} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
Definition array.hpp:882
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
T & Last()
Return the last element in the array.
Definition array.hpp:802
Abstract base class BilinearFormIntegrator.
virtual void AddMult(const Vector &x, Vector &y, const real_t val=1.) const
Matrix-Vector Multiplication y = y + val*A*x.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Partial matrix vector multiplication of (*this) with x involving only the rows given by rows....
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
int NumColBlocks() const
Return the number of column blocks.
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
void EliminateRowCols(const Array< int > &vdofs, BlockMatrix *Ae, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the rows and columns corresponding to the entries in vdofs + save the eliminated entries in...
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
int owns_blocks
If owns_blocks the SparseMatrix objects Aij will be deallocated.
virtual void Finalize(int skip_zeros=1)
Finalize all the submatrices.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A'*x.
int NumRowBlocks() const
Return the number of row blocks.
Class that performs static condensation of interior dofs for multiple FE spaces. This class is used i...
BlockMatrix & GetSchurMatrix()
Return the serial Schur complement matrix.
void FormSystemMatrix(Operator::DiagonalPolicy diag_policy)
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
void ReduceSystem(Vector &x, 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 &sc_sol, Vector &sol) const
void Finalize(int skip_zeros=0)
Finalize the construction of the Schur complement matrix.
void AssembleReducedSystem(int el, DenseMatrix &elmat, Vector &elvect)
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
virtual bool Factor(int m, real_t TOL=0.0)
Compute the Cholesky factorization of the current matrix.
void LSolve(int m, int n, real_t *X) const
Array2D< Array< BilinearFormIntegrator * > * > test_integs
Set of Test Space (broken) Integrators to be applied for matrix G.
Definition weakform.hpp:73
BlockStaticCondensation * static_cond
Owned.
Definition weakform.hpp:38
BlockMatrix * mat
Block matrix to be associated with the Block bilinear form. Owned.
Definition weakform.hpp:49
Array< int > tdof_offsets
Definition weakform.hpp:46
Array< FiniteElementSpace * > trial_fes
Trial FE spaces.
Definition weakform.hpp:60
virtual void Update()
Update the DPGWeakForm after mesh modifications (AMR)
Definition weakform.cpp:671
Array< Array< LinearFormIntegrator * > * > lfis
Set of Linear Form Integrators to be applied.
Definition weakform.hpp:76
void AddTrialIntegrator(BilinearFormIntegrator *bfi, int n, int m)
Adds new Trial Integrator. Assumes ownership of bfi.
Definition weakform.cpp:105
Array< FiniteElementCollection * > test_fecols
Test FE Collections (Broken)
Definition weakform.hpp:66
mfem::Operator::DiagonalPolicy diag_policy
Definition weakform.hpp:83
Array< DenseMatrix * > Bmat
Definition weakform.hpp:101
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Definition weakform.cpp:531
Vector & ComputeResidual(const BlockVector &x)
Compute DPG residual based error estimator.
Definition weakform.cpp:723
Array2D< Array< BilinearFormIntegrator * > * > trial_integs
Set of Trial Integrators to be applied for matrix B.
Definition weakform.hpp:70
virtual void BuildProlongation()
Definition weakform.cpp:133
Array< int > dof_offsets
Definition weakform.hpp:45
void AddDomainLFIntegrator(LinearFormIntegrator *lfi, int n)
Adds new Domain LF Integrator. Assumes ownership of bfi.
Definition weakform.cpp:125
void AllocMat()
Allocate appropriate BlockMatrix and assign it to mat.
Definition weakform.cpp:77
void ConformingAssemble()
Definition weakform.cpp:151
Array< Vector * > fvec
Definition weakform.hpp:102
void EnableStaticCondensation()
Definition weakform.cpp:717
BlockMatrix * R
Block Restriction.
Definition weakform.hpp:81
virtual void RecoverFEMSolution(const Vector &X, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
Definition weakform.cpp:598
Array< int > test_fecols_vdims
Definition weakform.hpp:67
virtual ~DPGWeakForm()
Definition weakform.cpp:811
void AddTestIntegrator(BilinearFormIntegrator *bfi, int n, int m)
Adds new Test Integrator. Assumes ownership of bfi.
Definition weakform.cpp:117
void ComputeOffsets()
Definition weakform.cpp:61
BlockVector * y
Block vector to be associated with the Block linear form.
Definition weakform.hpp:52
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this DPG weak form.
Definition weakform.cpp:476
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
Use the stored eliminated part of the matrix (see EliminateVDofs(const Array<int> &,...
Definition weakform.cpp:567
void ReleaseInitMemory()
Definition weakform.cpp:629
void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Definition weakform.cpp:97
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all integrators.
Definition weakform.cpp:247
BlockMatrix * mat_e
Block Matrix used to store the eliminations from the b.c. Owned. .
Definition weakform.hpp:57
void EliminateVDofs(const Array< int > &vdofs, Operator::DiagonalPolicy dpolicy=Operator::DIAG_ONE)
Eliminate the given vdofs, storing the eliminated part internally in .
Definition weakform.cpp:574
Array< int > IsTraceFes
Flags to determine if a FiniteElementSpace is Trace.
Definition weakform.hpp:63
BlockMatrix * P
Block Prolongation.
Definition weakform.hpp:79
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:262
void GetSubMatrix(const Array< int > &idx, DenseMatrix &A) const
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:115
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void AddSubMatrix(const Array< int > &idx, const DenseMatrix &A)
(*this)(idx[i],idx[j]) += A(i,j)
void TransformDual(real_t *v) const
Definition doftrans.cpp:81
void InvTransformPrimal(real_t *v) const
Definition doftrans.cpp:49
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition eltrans.hpp:162
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:484
Abstract class for all finite elements.
Definition fe_base.hpp:239
Abstract base class LinearFormIntegrator.
Definition lininteg.hpp:25
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition mesh.cpp:980
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1438
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition mesh.cpp:7201
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition mesh.cpp:6985
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
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
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
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Data type sparse matrix.
Definition sparsemat.hpp:51
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Vector data type.
Definition vector.hpp:80
void SetVector(const Vector &v, int offset)
Definition vector.cpp:274
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:175
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition vector.cpp:670
void AddSubVector(const Vector &v, int offset)
Definition vector.cpp:287
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
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
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
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 GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
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
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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 Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
float real_t
Definition config.hpp:43