1 // Copyright (c) 2010-2021, 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 <cmath>
17
18 namespace mfem
19 {
20
AllocMat()21 void BilinearForm::AllocMat()
22 {
23 if (static_cond) { return; }
24
25 if (precompute_sparsity == 0 || fes->GetVDim() > 1)
26 {
27 mat = new SparseMatrix(height);
28 return;
29 }
30
31 const Table &elem_dof = fes->GetElementToDofTable();
32 Table dof_dof;
33
34 if (interior_face_integs.Size() > 0)
35 {
36 // the sparsity pattern is defined from the map: face->element->dof
37 Table face_dof, dof_face;
38 {
39 Table *face_elem = fes->GetMesh()->GetFaceToElementTable();
40 mfem::Mult(*face_elem, elem_dof, face_dof);
41 delete face_elem;
42 }
43 Transpose(face_dof, dof_face, height);
44 mfem::Mult(dof_face, face_dof, dof_dof);
45 }
46 else
47 {
48 // the sparsity pattern is defined from the map: element->dof
49 Table dof_elem;
50 Transpose(elem_dof, dof_elem, height);
51 mfem::Mult(dof_elem, elem_dof, dof_dof);
52 }
53
54 dof_dof.SortRows();
55
56 int *I = dof_dof.GetI();
57 int *J = dof_dof.GetJ();
58 double *data = Memory<double>(I[height]);
59
60 mat = new SparseMatrix(I, J, data, height, height, true, true, true);
61 *mat = 0.0;
62
63 dof_dof.LoseData();
64 }
65
BilinearForm(FiniteElementSpace * f)66 BilinearForm::BilinearForm(FiniteElementSpace * f)
67 : Matrix (f->GetVSize())
68 {
69 fes = f;
70 sequence = f->GetSequence();
71 mat = mat_e = NULL;
72 extern_bfs = 0;
73 element_matrices = NULL;
74 static_cond = NULL;
75 hybridization = NULL;
76 precompute_sparsity = 0;
77 diag_policy = DIAG_KEEP;
78
79 assembly = AssemblyLevel::LEGACY;
80 batch = 1;
81 ext = NULL;
82 }
83
BilinearForm(FiniteElementSpace * f,BilinearForm * bf,int ps)84 BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
85 : Matrix (f->GetVSize())
86 {
87 fes = f;
88 sequence = f->GetSequence();
89 mat_e = NULL;
90 extern_bfs = 1;
91 element_matrices = NULL;
92 static_cond = NULL;
93 hybridization = NULL;
94 precompute_sparsity = ps;
95 diag_policy = DIAG_KEEP;
96
97 assembly = AssemblyLevel::LEGACY;
98 batch = 1;
99 ext = NULL;
100
101 // Copy the pointers to the integrators
102 domain_integs = bf->domain_integs;
103
104 boundary_integs = bf->boundary_integs;
105 boundary_integs_marker = bf->boundary_integs_marker;
106
107 interior_face_integs = bf->interior_face_integs;
108
109 boundary_face_integs = bf->boundary_face_integs;
110 boundary_face_integs_marker = bf->boundary_face_integs_marker;
111
112 AllocMat();
113 }
114
SetAssemblyLevel(AssemblyLevel assembly_level)115 void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
116 {
117 if (ext)
118 {
119 MFEM_ABORT("the assembly level has already been set!");
120 }
121 assembly = assembly_level;
122 switch (assembly)
123 {
124 case AssemblyLevel::LEGACY:
125 break;
126 case AssemblyLevel::FULL:
127 ext = new FABilinearFormExtension(this);
128 break;
129 case AssemblyLevel::ELEMENT:
130 ext = new EABilinearFormExtension(this);
131 break;
132 case AssemblyLevel::PARTIAL:
133 ext = new PABilinearFormExtension(this);
134 break;
135 case AssemblyLevel::NONE:
136 ext = new MFBilinearFormExtension(this);
137 break;
138 default:
139 mfem_error("Unknown assembly level");
140 }
141 }
142
EnableStaticCondensation()143 void BilinearForm::EnableStaticCondensation()
144 {
145 delete static_cond;
146 if (assembly != AssemblyLevel::LEGACY)
147 {
148 static_cond = NULL;
149 MFEM_WARNING("Static condensation not supported for this assembly level");
150 return;
151 }
152 static_cond = new StaticCondensation(fes);
153 if (static_cond->ReducesTrueVSize())
154 {
155 bool symmetric = false; // TODO
156 bool block_diagonal = false; // TODO
157 static_cond->Init(symmetric, block_diagonal);
158 }
159 else
160 {
161 delete static_cond;
162 static_cond = NULL;
163 }
164 }
165
EnableHybridization(FiniteElementSpace * constr_space,BilinearFormIntegrator * constr_integ,const Array<int> & ess_tdof_list)166 void BilinearForm::EnableHybridization(FiniteElementSpace *constr_space,
167 BilinearFormIntegrator *constr_integ,
168 const Array<int> &ess_tdof_list)
169 {
170 delete hybridization;
171 if (assembly != AssemblyLevel::LEGACY)
172 {
173 delete constr_integ;
174 hybridization = NULL;
175 MFEM_WARNING("Hybridization not supported for this assembly level");
176 return;
177 }
178 hybridization = new Hybridization(fes, constr_space);
179 hybridization->SetConstraintIntegrator(constr_integ);
180 hybridization->Init(ess_tdof_list);
181 }
182
UseSparsity(int * I,int * J,bool isSorted)183 void BilinearForm::UseSparsity(int *I, int *J, bool isSorted)
184 {
185 if (static_cond) { return; }
186
187 if (mat)
188 {
189 if (mat->Finalized() && mat->GetI() == I && mat->GetJ() == J)
190 {
191 return; // mat is already using the given sparsity
192 }
193 delete mat;
194 }
195 height = width = fes->GetVSize();
196 mat = new SparseMatrix(I, J, NULL, height, width, false, true, isSorted);
197 }
198
UseSparsity(SparseMatrix & A)199 void BilinearForm::UseSparsity(SparseMatrix &A)
200 {
201 MFEM_ASSERT(A.Height() == fes->GetVSize() && A.Width() == fes->GetVSize(),
202 "invalid matrix A dimensions: "
203 << A.Height() << " x " << A.Width());
204 MFEM_ASSERT(A.Finalized(), "matrix A must be Finalized");
205
206 UseSparsity(A.GetI(), A.GetJ(), A.ColumnsAreSorted());
207 }
208
Elem(int i,int j)209 double& BilinearForm::Elem (int i, int j)
210 {
211 return mat -> Elem(i,j);
212 }
213
Elem(int i,int j) const214 const double& BilinearForm::Elem (int i, int j) const
215 {
216 return mat -> Elem(i,j);
217 }
218
Inverse() const219 MatrixInverse * BilinearForm::Inverse() const
220 {
221 return mat -> Inverse();
222 }
223
Finalize(int skip_zeros)224 void BilinearForm::Finalize (int skip_zeros)
225 {
226 if (assembly == AssemblyLevel::LEGACY)
227 {
228 if (!static_cond) { mat->Finalize(skip_zeros); }
229 if (mat_e) { mat_e->Finalize(skip_zeros); }
230 if (static_cond) { static_cond->Finalize(); }
231 if (hybridization) { hybridization->Finalize(); }
232 }
233 }
234
AddDomainIntegrator(BilinearFormIntegrator * bfi)235 void BilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi)
236 {
237 domain_integs.Append(bfi);
238 domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
239 }
240
AddDomainIntegrator(BilinearFormIntegrator * bfi,Array<int> & elem_marker)241 void BilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi,
242 Array<int> &elem_marker)
243 {
244 domain_integs.Append(bfi);
245 domain_integs_marker.Append(&elem_marker);
246 }
247
AddBoundaryIntegrator(BilinearFormIntegrator * bfi)248 void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
249 {
250 boundary_integs.Append (bfi);
251 boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
252 }
253
AddBoundaryIntegrator(BilinearFormIntegrator * bfi,Array<int> & bdr_marker)254 void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
255 Array<int> &bdr_marker)
256 {
257 boundary_integs.Append (bfi);
258 boundary_integs_marker.Append(&bdr_marker);
259 }
260
AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi)261 void BilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi)
262 {
263 interior_face_integs.Append (bfi);
264 }
265
AddBdrFaceIntegrator(BilinearFormIntegrator * bfi)266 void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
267 {
268 boundary_face_integs.Append(bfi);
269 // NULL marker means apply everywhere
270 boundary_face_integs_marker.Append(NULL);
271 }
272
AddBdrFaceIntegrator(BilinearFormIntegrator * bfi,Array<int> & bdr_marker)273 void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
274 Array<int> &bdr_marker)
275 {
276 boundary_face_integs.Append(bfi);
277 boundary_face_integs_marker.Append(&bdr_marker);
278 }
279
ComputeElementMatrix(int i,DenseMatrix & elmat)280 void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
281 {
282 if (element_matrices)
283 {
284 elmat.SetSize(element_matrices->SizeI(), element_matrices->SizeJ());
285 elmat = element_matrices->GetData(i);
286 return;
287 }
288
289 if (domain_integs.Size())
290 {
291 const FiniteElement &fe = *fes->GetFE(i);
292 ElementTransformation *eltrans = fes->GetElementTransformation(i);
293 domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
294 for (int k = 1; k < domain_integs.Size(); k++)
295 {
296 domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
297 elmat += elemmat;
298 }
299 }
300 else
301 {
302 fes->GetElementVDofs(i, vdofs);
303 elmat.SetSize(vdofs.Size());
304 elmat = 0.0;
305 }
306 }
307
ComputeBdrElementMatrix(int i,DenseMatrix & elmat)308 void BilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
309 {
310 if (boundary_integs.Size())
311 {
312 const FiniteElement &be = *fes->GetBE(i);
313 ElementTransformation *eltrans = fes->GetBdrElementTransformation(i);
314 boundary_integs[0]->AssembleElementMatrix(be, *eltrans, elmat);
315 for (int k = 1; k < boundary_integs.Size(); k++)
316 {
317 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
318 elmat += elemmat;
319 }
320 }
321 else
322 {
323 fes->GetBdrElementVDofs(i, vdofs);
324 elmat.SetSize(vdofs.Size());
325 elmat = 0.0;
326 }
327 }
328
AssembleElementMatrix(int i,const DenseMatrix & elmat,int skip_zeros)329 void BilinearForm::AssembleElementMatrix(
330 int i, const DenseMatrix &elmat, int skip_zeros)
331 {
332 AssembleElementMatrix(i, elmat, vdofs, skip_zeros);
333 }
334
AssembleElementMatrix(int i,const DenseMatrix & elmat,Array<int> & vdofs,int skip_zeros)335 void BilinearForm::AssembleElementMatrix(
336 int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
337 {
338 fes->GetElementVDofs(i, vdofs);
339 if (static_cond)
340 {
341 static_cond->AssembleMatrix(i, elmat);
342 }
343 else
344 {
345 if (mat == NULL)
346 {
347 AllocMat();
348 }
349 mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
350 if (hybridization)
351 {
352 hybridization->AssembleMatrix(i, elmat);
353 }
354 }
355 }
356
AssembleBdrElementMatrix(int i,const DenseMatrix & elmat,int skip_zeros)357 void BilinearForm::AssembleBdrElementMatrix(
358 int i, const DenseMatrix &elmat, int skip_zeros)
359 {
360 AssembleBdrElementMatrix(i, elmat, vdofs, skip_zeros);
361 }
362
AssembleBdrElementMatrix(int i,const DenseMatrix & elmat,Array<int> & vdofs,int skip_zeros)363 void BilinearForm::AssembleBdrElementMatrix(
364 int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
365 {
366 fes->GetBdrElementVDofs(i, vdofs);
367 if (static_cond)
368 {
369 static_cond->AssembleBdrMatrix(i, elmat);
370 }
371 else
372 {
373 if (mat == NULL)
374 {
375 AllocMat();
376 }
377 mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
378 if (hybridization)
379 {
380 hybridization->AssembleBdrMatrix(i, elmat);
381 }
382 }
383 }
384
Assemble(int skip_zeros)385 void BilinearForm::Assemble(int skip_zeros)
386 {
387 if (ext)
388 {
389 ext->Assemble();
390 return;
391 }
392
393 ElementTransformation *eltrans;
394 Mesh *mesh = fes -> GetMesh();
395 DenseMatrix elmat, *elmat_p;
396
397 if (mat == NULL)
398 {
399 AllocMat();
400 }
401
402 #ifdef MFEM_USE_LEGACY_OPENMP
403 int free_element_matrices = 0;
404 if (!element_matrices)
405 {
406 ComputeElementMatrices();
407 free_element_matrices = 1;
408 }
409 #endif
410
411 if (domain_integs.Size())
412 {
413 for (int k = 0; k < domain_integs.Size(); k++)
414 {
415 if (domain_integs_marker[k] != NULL)
416 {
417 MFEM_VERIFY(mesh->attributes.Size() ==
418 domain_integs_marker[k]->Size(),
419 "invalid element marker for domain integrator #"
420 << k << ", counting from zero");
421 }
422 }
423
424 for (int i = 0; i < fes -> GetNE(); i++)
425 {
426 int elem_attr = fes->GetMesh()->GetAttribute(i);
427 fes->GetElementVDofs(i, vdofs);
428 if (element_matrices)
429 {
430 elmat_p = &(*element_matrices)(i);
431 }
432 else
433 {
434 elmat.SetSize(0);
435 for (int k = 0; k < domain_integs.Size(); k++)
436 {
437 if ( domain_integs_marker[k] == NULL ||
438 (*(domain_integs_marker[k]))[elem_attr-1] == 1)
439 {
440 const FiniteElement &fe = *fes->GetFE(i);
441 eltrans = fes->GetElementTransformation(i);
442 domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
443 if (elmat.Size() == 0)
444 {
445 elmat = elemmat;
446 }
447 else
448 {
449 elmat += elemmat;
450 }
451 }
452 }
453 if (elmat.Size() == 0)
454 {
455 continue;
456 }
457 else
458 {
459 elmat_p = &elmat;
460 }
461 }
462 if (static_cond)
463 {
464 static_cond->AssembleMatrix(i, *elmat_p);
465 }
466 else
467 {
468 mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
469 if (hybridization)
470 {
471 hybridization->AssembleMatrix(i, *elmat_p);
472 }
473 }
474 }
475 }
476
477 if (boundary_integs.Size())
478 {
479 // Which boundary attributes need to be processed?
480 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
481 mesh->bdr_attributes.Max() : 0);
482 bdr_attr_marker = 0;
483 for (int k = 0; k < boundary_integs.Size(); k++)
484 {
485 if (boundary_integs_marker[k] == NULL)
486 {
487 bdr_attr_marker = 1;
488 break;
489 }
490 Array<int> &bdr_marker = *boundary_integs_marker[k];
491 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
492 "invalid boundary marker for boundary integrator #"
493 << k << ", counting from zero");
494 for (int i = 0; i < bdr_attr_marker.Size(); i++)
495 {
496 bdr_attr_marker[i] |= bdr_marker[i];
497 }
498 }
499
500 for (int i = 0; i < fes -> GetNBE(); i++)
501 {
502 const int bdr_attr = mesh->GetBdrAttribute(i);
503 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
504
505 const FiniteElement &be = *fes->GetBE(i);
506 fes -> GetBdrElementVDofs (i, vdofs);
507 eltrans = fes -> GetBdrElementTransformation (i);
508 int k = 0;
509 for (; k < boundary_integs.Size(); k++)
510 {
511 if (boundary_integs_marker[k] &&
512 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
513
514 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elmat);
515 k++;
516 break;
517 }
518 for (; k < boundary_integs.Size(); k++)
519 {
520 if (boundary_integs_marker[k] &&
521 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
522
523 boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
524 elmat += elemmat;
525 }
526 if (!static_cond)
527 {
528 mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
529 if (hybridization)
530 {
531 hybridization->AssembleBdrMatrix(i, elmat);
532 }
533 }
534 else
535 {
536 static_cond->AssembleBdrMatrix(i, elmat);
537 }
538 }
539 }
540
541 if (interior_face_integs.Size())
542 {
543 FaceElementTransformations *tr;
544 Array<int> vdofs2;
545
546 int nfaces = mesh->GetNumFaces();
547 for (int i = 0; i < nfaces; i++)
548 {
549 tr = mesh -> GetInteriorFaceTransformations (i);
550 if (tr != NULL)
551 {
552 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
553 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
554 vdofs.Append (vdofs2);
555 for (int k = 0; k < interior_face_integs.Size(); k++)
556 {
557 interior_face_integs[k]->
558 AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
559 *fes->GetFE(tr->Elem2No),
560 *tr, elemmat);
561 mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
562 }
563 }
564 }
565 }
566
567 if (boundary_face_integs.Size())
568 {
569 FaceElementTransformations *tr;
570 const FiniteElement *fe1, *fe2;
571
572 // Which boundary attributes need to be processed?
573 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
574 mesh->bdr_attributes.Max() : 0);
575 bdr_attr_marker = 0;
576 for (int k = 0; k < boundary_face_integs.Size(); k++)
577 {
578 if (boundary_face_integs_marker[k] == NULL)
579 {
580 bdr_attr_marker = 1;
581 break;
582 }
583 Array<int> &bdr_marker = *boundary_face_integs_marker[k];
584 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
585 "invalid boundary marker for boundary face integrator #"
586 << k << ", counting from zero");
587 for (int i = 0; i < bdr_attr_marker.Size(); i++)
588 {
589 bdr_attr_marker[i] |= bdr_marker[i];
590 }
591 }
592
593 for (int i = 0; i < fes -> GetNBE(); i++)
594 {
595 const int bdr_attr = mesh->GetBdrAttribute(i);
596 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
597
598 tr = mesh -> GetBdrFaceTransformations (i);
599 if (tr != NULL)
600 {
601 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
602 fe1 = fes -> GetFE (tr -> Elem1No);
603 // The fe2 object is really a dummy and not used on the boundaries,
604 // but we can't dereference a NULL pointer, and we don't want to
605 // actually make a fake element.
606 fe2 = fe1;
607 for (int k = 0; k < boundary_face_integs.Size(); k++)
608 {
609 if (boundary_face_integs_marker[k] &&
610 (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
611 { continue; }
612
613 boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
614 elemmat);
615 mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
616 }
617 }
618 }
619 }
620
621 #ifdef MFEM_USE_LEGACY_OPENMP
622 if (free_element_matrices)
623 {
624 FreeElementMatrices();
625 }
626 #endif
627 }
628
ConformingAssemble()629 void BilinearForm::ConformingAssemble()
630 {
631 // Do not remove zero entries to preserve the symmetric structure of the
632 // matrix which in turn will give rise to symmetric structure in the new
633 // matrix. This ensures that subsequent calls to EliminateRowCol will work
634 // correctly.
635 Finalize(0);
636 MFEM_ASSERT(mat, "the BilinearForm is not assembled");
637
638 const SparseMatrix *P = fes->GetConformingProlongation();
639 if (!P) { return; } // conforming mesh
640
641 SparseMatrix *R = Transpose(*P);
642 SparseMatrix *RA = mfem::Mult(*R, *mat);
643 delete mat;
644 if (mat_e)
645 {
646 SparseMatrix *RAe = mfem::Mult(*R, *mat_e);
647 delete mat_e;
648 mat_e = RAe;
649 }
650 delete R;
651 mat = mfem::Mult(*RA, *P);
652 delete RA;
653 if (mat_e)
654 {
655 SparseMatrix *RAeP = mfem::Mult(*mat_e, *P);
656 delete mat_e;
657 mat_e = RAeP;
658 }
659
660 height = mat->Height();
661 width = mat->Width();
662 }
663
AssembleDiagonal(Vector & diag) const664 void BilinearForm::AssembleDiagonal(Vector &diag) const
665 {
666 MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
667 "Vector for holding diagonal has wrong size!");
668 const SparseMatrix *cP = fes->GetConformingProlongation();
669 if (!ext)
670 {
671 MFEM_ASSERT(mat, "the BilinearForm is not assembled!");
672 MFEM_ASSERT(cP == nullptr || mat->Height() == cP->Width(),
673 "BilinearForm::ConformingAssemble() is not called!");
674 mat->GetDiag(diag);
675 return;
676 }
677 // Here, we have extension, ext.
678 if (!cP)
679 {
680 ext->AssembleDiagonal(diag);
681 return;
682 }
683 // Here, we have extension, ext, and conforming prolongation, cP.
684
685 // For an AMR mesh, a convergent diagonal is assembled with |P^T| d_l,
686 // where |P^T| has the entry-wise absolute values of the conforming
687 // prolongation transpose operator.
688 Vector local_diag(cP->Height());
689 ext->AssembleDiagonal(local_diag);
690 cP->AbsMultTranspose(local_diag, diag);
691 }
692
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B,int copy_interior)693 void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
694 Vector &b, OperatorHandle &A, Vector &X,
695 Vector &B, int copy_interior)
696 {
697 if (ext)
698 {
699 ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
700 return;
701 }
702 const SparseMatrix *P = fes->GetConformingProlongation();
703 FormSystemMatrix(ess_tdof_list, A);
704
705 // Transform the system and perform the elimination in B, based on the
706 // essential BC values from x. Restrict the BC part of x in X, and set the
707 // non-BC part to zero. Since there is no good initial guess for the Lagrange
708 // multipliers, set X = 0.0 for hybridization.
709 if (static_cond)
710 {
711 // Schur complement reduction to the exposed dofs
712 static_cond->ReduceSystem(x, b, X, B, copy_interior);
713 }
714 else if (!P) // conforming space
715 {
716 if (hybridization)
717 {
718 // Reduction to the Lagrange multipliers system
719 EliminateVDofsInRHS(ess_tdof_list, x, b);
720 hybridization->ReduceRHS(b, B);
721 X.SetSize(B.Size());
722 X = 0.0;
723 }
724 else
725 {
726 // A, X and B point to the same data as mat, x and b
727 EliminateVDofsInRHS(ess_tdof_list, x, b);
728 X.MakeRef(x, 0, x.Size());
729 B.MakeRef(b, 0, b.Size());
730 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
731 }
732 }
733 else // non-conforming space
734 {
735 if (hybridization)
736 {
737 // Reduction to the Lagrange multipliers system
738 const SparseMatrix *R = fes->GetConformingRestriction();
739 Vector conf_b(P->Width()), conf_x(P->Width());
740 P->MultTranspose(b, conf_b);
741 R->Mult(x, conf_x);
742 EliminateVDofsInRHS(ess_tdof_list, conf_x, conf_b);
743 R->MultTranspose(conf_b, b); // store eliminated rhs in b
744 hybridization->ReduceRHS(conf_b, B);
745 X.SetSize(B.Size());
746 X = 0.0;
747 }
748 else
749 {
750 // Variational restriction with P
751 const SparseMatrix *R = fes->GetConformingRestriction();
752 B.SetSize(P->Width());
753 P->MultTranspose(b, B);
754 X.SetSize(R->Height());
755 R->Mult(x, X);
756 EliminateVDofsInRHS(ess_tdof_list, X, B);
757 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
758 }
759 }
760 }
761
FormSystemMatrix(const Array<int> & ess_tdof_list,OperatorHandle & A)762 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
763 OperatorHandle &A)
764 {
765 if (ext)
766 {
767 ext->FormSystemMatrix(ess_tdof_list, A);
768 return;
769 }
770
771 // Finish the matrix assembly and perform BC elimination, storing the
772 // eliminated part of the matrix.
773 if (static_cond)
774 {
775 if (!static_cond->HasEliminatedBC())
776 {
777 static_cond->SetEssentialTrueDofs(ess_tdof_list);
778 static_cond->Finalize(); // finalize Schur complement (to true dofs)
779 static_cond->EliminateReducedTrueDofs(diag_policy);
780 static_cond->Finalize(); // finalize eliminated part
781 }
782 A.Reset(&static_cond->GetMatrix(), false);
783 }
784 else
785 {
786 if (!mat_e)
787 {
788 const SparseMatrix *P = fes->GetConformingProlongation();
789 if (P) { ConformingAssemble(); }
790 EliminateVDofs(ess_tdof_list, diag_policy);
791 const int remove_zeros = 0;
792 Finalize(remove_zeros);
793 }
794 if (hybridization)
795 {
796 A.Reset(&hybridization->GetMatrix(), false);
797 }
798 else
799 {
800 A.Reset(mat, false);
801 }
802 }
803 }
804
RecoverFEMSolution(const Vector & X,const Vector & b,Vector & x)805 void BilinearForm::RecoverFEMSolution(const Vector &X,
806 const Vector &b, Vector &x)
807 {
808 if (ext)
809 {
810 ext->RecoverFEMSolution(X, b, x);
811 return;
812 }
813
814 const SparseMatrix *P = fes->GetConformingProlongation();
815 if (!P) // conforming space
816 {
817 if (static_cond)
818 {
819 // Private dofs back solve
820 static_cond->ComputeSolution(b, X, x);
821 }
822 else if (hybridization)
823 {
824 // Primal unknowns recovery
825 hybridization->ComputeSolution(b, X, x);
826 }
827 else
828 {
829 // X and x point to the same data
830
831 // If the validity flags of X's Memory were changed (e.g. if it was
832 // moved to device memory) then we need to tell x about that.
833 x.SyncMemory(X);
834 }
835 }
836 else // non-conforming space
837 {
838 if (static_cond)
839 {
840 // Private dofs back solve
841 static_cond->ComputeSolution(b, X, x);
842 }
843 else if (hybridization)
844 {
845 // Primal unknowns recovery
846 Vector conf_b(P->Width()), conf_x(P->Width());
847 P->MultTranspose(b, conf_b);
848 const SparseMatrix *R = fes->GetConformingRestriction();
849 R->Mult(x, conf_x); // get essential b.c. from x
850 hybridization->ComputeSolution(conf_b, X, conf_x);
851 x.SetSize(P->Height());
852 P->Mult(conf_x, x);
853 }
854 else
855 {
856 // Apply conforming prolongation
857 x.SetSize(P->Height());
858 P->Mult(X, x);
859 }
860 }
861 }
862
ComputeElementMatrices()863 void BilinearForm::ComputeElementMatrices()
864 {
865 if (element_matrices || domain_integs.Size() == 0 || fes->GetNE() == 0)
866 {
867 return;
868 }
869
870 int num_elements = fes->GetNE();
871 int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
872
873 element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
874 num_elements);
875
876 DenseMatrix tmp;
877 IsoparametricTransformation eltrans;
878
879 #ifdef MFEM_USE_LEGACY_OPENMP
880 #pragma omp parallel for private(tmp,eltrans)
881 #endif
882 for (int i = 0; i < num_elements; i++)
883 {
884 DenseMatrix elmat(element_matrices->GetData(i),
885 num_dofs_per_el, num_dofs_per_el);
886 const FiniteElement &fe = *fes->GetFE(i);
887 #ifdef MFEM_DEBUG
888 if (num_dofs_per_el != fe.GetDof()*fes->GetVDim())
889 mfem_error("BilinearForm::ComputeElementMatrices:"
890 " all elements must have same number of dofs");
891 #endif
892 fes->GetElementTransformation(i, &eltrans);
893
894 domain_integs[0]->AssembleElementMatrix(fe, eltrans, elmat);
895 for (int k = 1; k < domain_integs.Size(); k++)
896 {
897 // note: some integrators may not be thread-safe
898 domain_integs[k]->AssembleElementMatrix(fe, eltrans, tmp);
899 elmat += tmp;
900 }
901 elmat.ClearExternalData();
902 }
903 }
904
EliminateEssentialBC(const Array<int> & bdr_attr_is_ess,const Vector & sol,Vector & rhs,DiagonalPolicy dpolicy)905 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
906 const Vector &sol, Vector &rhs,
907 DiagonalPolicy dpolicy)
908 {
909 Array<int> ess_dofs, conf_ess_dofs;
910 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
911
912 if (fes->GetVSize() == height)
913 {
914 EliminateEssentialBCFromDofs(ess_dofs, sol, rhs, dpolicy);
915 }
916 else
917 {
918 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
919 EliminateEssentialBCFromDofs(conf_ess_dofs, sol, rhs, dpolicy);
920 }
921 }
922
EliminateEssentialBC(const Array<int> & bdr_attr_is_ess,DiagonalPolicy dpolicy)923 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
924 DiagonalPolicy dpolicy)
925 {
926 Array<int> ess_dofs, conf_ess_dofs;
927 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
928
929 if (fes->GetVSize() == height)
930 {
931 EliminateEssentialBCFromDofs(ess_dofs, dpolicy);
932 }
933 else
934 {
935 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
936 EliminateEssentialBCFromDofs(conf_ess_dofs, dpolicy);
937 }
938 }
939
EliminateEssentialBCDiag(const Array<int> & bdr_attr_is_ess,double value)940 void BilinearForm::EliminateEssentialBCDiag (const Array<int> &bdr_attr_is_ess,
941 double value)
942 {
943 Array<int> ess_dofs, conf_ess_dofs;
944 fes->GetEssentialVDofs(bdr_attr_is_ess, ess_dofs);
945
946 if (fes->GetVSize() == height)
947 {
948 EliminateEssentialBCFromDofsDiag(ess_dofs, value);
949 }
950 else
951 {
952 fes->GetRestrictionMatrix()->BooleanMult(ess_dofs, conf_ess_dofs);
953 EliminateEssentialBCFromDofsDiag(conf_ess_dofs, value);
954 }
955 }
956
EliminateVDofs(const Array<int> & vdofs,const Vector & sol,Vector & rhs,DiagonalPolicy dpolicy)957 void BilinearForm::EliminateVDofs(const Array<int> &vdofs,
958 const Vector &sol, Vector &rhs,
959 DiagonalPolicy dpolicy)
960 {
961 for (int i = 0; i < vdofs.Size(); i++)
962 {
963 int vdof = vdofs[i];
964 if ( vdof >= 0 )
965 {
966 mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
967 }
968 else
969 {
970 mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
971 }
972 }
973 }
974
EliminateVDofs(const Array<int> & vdofs,DiagonalPolicy dpolicy)975 void BilinearForm::EliminateVDofs(const Array<int> &vdofs,
976 DiagonalPolicy dpolicy)
977 {
978 if (mat_e == NULL)
979 {
980 mat_e = new SparseMatrix(height);
981 }
982
983 for (int i = 0; i < vdofs.Size(); i++)
984 {
985 int vdof = vdofs[i];
986 if ( vdof >= 0 )
987 {
988 mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
989 }
990 else
991 {
992 mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
993 }
994 }
995 }
996
EliminateEssentialBCFromDofs(const Array<int> & ess_dofs,const Vector & sol,Vector & rhs,DiagonalPolicy dpolicy)997 void BilinearForm::EliminateEssentialBCFromDofs(
998 const Array<int> &ess_dofs, const Vector &sol, Vector &rhs,
999 DiagonalPolicy dpolicy)
1000 {
1001 MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
1002 MFEM_ASSERT(sol.Size() == height, "incorrect sol Vector size");
1003 MFEM_ASSERT(rhs.Size() == height, "incorrect rhs Vector size");
1004
1005 for (int i = 0; i < ess_dofs.Size(); i++)
1006 if (ess_dofs[i] < 0)
1007 {
1008 mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
1009 }
1010 }
1011
EliminateEssentialBCFromDofs(const Array<int> & ess_dofs,DiagonalPolicy dpolicy)1012 void BilinearForm::EliminateEssentialBCFromDofs (const Array<int> &ess_dofs,
1013 DiagonalPolicy dpolicy)
1014 {
1015 MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
1016
1017 for (int i = 0; i < ess_dofs.Size(); i++)
1018 if (ess_dofs[i] < 0)
1019 {
1020 mat -> EliminateRowCol (i, dpolicy);
1021 }
1022 }
1023
EliminateEssentialBCFromDofsDiag(const Array<int> & ess_dofs,double value)1024 void BilinearForm::EliminateEssentialBCFromDofsDiag (const Array<int> &ess_dofs,
1025 double value)
1026 {
1027 MFEM_ASSERT(ess_dofs.Size() == height, "incorrect dof Array size");
1028
1029 for (int i = 0; i < ess_dofs.Size(); i++)
1030 if (ess_dofs[i] < 0)
1031 {
1032 mat -> EliminateRowColDiag (i, value);
1033 }
1034 }
1035
EliminateVDofsInRHS(const Array<int> & vdofs,const Vector & x,Vector & b)1036 void BilinearForm::EliminateVDofsInRHS(
1037 const Array<int> &vdofs, const Vector &x, Vector &b)
1038 {
1039 mat_e->AddMult(x, b, -1.);
1040 mat->PartMult(vdofs, x, b);
1041 }
1042
Mult(const Vector & x,Vector & y) const1043 void BilinearForm::Mult(const Vector &x, Vector &y) const
1044 {
1045 if (ext)
1046 {
1047 ext->Mult(x, y);
1048 }
1049 else
1050 {
1051 mat->Mult(x, y);
1052 }
1053 }
1054
Update(FiniteElementSpace * nfes)1055 void BilinearForm::Update(FiniteElementSpace *nfes)
1056 {
1057 bool full_update;
1058
1059 if (nfes && nfes != fes)
1060 {
1061 full_update = true;
1062 fes = nfes;
1063 }
1064 else
1065 {
1066 // Check for different size (e.g. assembled form on non-conforming space)
1067 // or different sequence number.
1068 full_update = (fes->GetVSize() != Height() ||
1069 sequence < fes->GetSequence());
1070 }
1071
1072 delete mat_e;
1073 mat_e = NULL;
1074 FreeElementMatrices();
1075 delete static_cond;
1076 static_cond = NULL;
1077
1078 if (full_update)
1079 {
1080 delete mat;
1081 mat = NULL;
1082 delete hybridization;
1083 hybridization = NULL;
1084 sequence = fes->GetSequence();
1085 }
1086 else
1087 {
1088 if (mat) { *mat = 0.0; }
1089 if (hybridization) { hybridization->Reset(); }
1090 }
1091
1092 height = width = fes->GetVSize();
1093
1094 if (ext) { ext->Update(); }
1095 }
1096
SetDiagonalPolicy(DiagonalPolicy policy)1097 void BilinearForm::SetDiagonalPolicy(DiagonalPolicy policy)
1098 {
1099 diag_policy = policy;
1100 }
1101
~BilinearForm()1102 BilinearForm::~BilinearForm()
1103 {
1104 delete mat_e;
1105 delete mat;
1106 delete element_matrices;
1107 delete static_cond;
1108 delete hybridization;
1109
1110 if (!extern_bfs)
1111 {
1112 int k;
1113 for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
1114 for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
1115 for (k=0; k < interior_face_integs.Size(); k++)
1116 { delete interior_face_integs[k]; }
1117 for (k=0; k < boundary_face_integs.Size(); k++)
1118 { delete boundary_face_integs[k]; }
1119 }
1120
1121 delete ext;
1122 }
1123
1124
MixedBilinearForm(FiniteElementSpace * tr_fes,FiniteElementSpace * te_fes)1125 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1126 FiniteElementSpace *te_fes)
1127 : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1128 {
1129 trial_fes = tr_fes;
1130 test_fes = te_fes;
1131 mat = NULL;
1132 mat_e = NULL;
1133 extern_bfs = 0;
1134 assembly = AssemblyLevel::LEGACY;
1135 ext = NULL;
1136 }
1137
MixedBilinearForm(FiniteElementSpace * tr_fes,FiniteElementSpace * te_fes,MixedBilinearForm * mbf)1138 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
1139 FiniteElementSpace *te_fes,
1140 MixedBilinearForm * mbf)
1141 : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1142 {
1143 trial_fes = tr_fes;
1144 test_fes = te_fes;
1145 mat = NULL;
1146 mat_e = NULL;
1147 extern_bfs = 1;
1148 ext = NULL;
1149
1150 // Copy the pointers to the integrators
1151 domain_integs = mbf->domain_integs;
1152 boundary_integs = mbf->boundary_integs;
1153 trace_face_integs = mbf->trace_face_integs;
1154 boundary_trace_face_integs = mbf->boundary_trace_face_integs;
1155
1156 boundary_integs_marker = mbf->boundary_integs_marker;
1157 boundary_trace_face_integs_marker = mbf->boundary_trace_face_integs_marker;
1158
1159 assembly = AssemblyLevel::LEGACY;
1160 ext = NULL;
1161 }
1162
SetAssemblyLevel(AssemblyLevel assembly_level)1163 void MixedBilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
1164 {
1165 if (ext)
1166 {
1167 MFEM_ABORT("the assembly level has already been set!");
1168 }
1169 assembly = assembly_level;
1170 switch (assembly)
1171 {
1172 case AssemblyLevel::LEGACY:
1173 break;
1174 case AssemblyLevel::FULL:
1175 // ext = new FAMixedBilinearFormExtension(this);
1176 // Use the original BilinearForm implementation for now
1177 break;
1178 case AssemblyLevel::ELEMENT:
1179 mfem_error("Element assembly not supported yet... stay tuned!");
1180 // ext = new EAMixedBilinearFormExtension(this);
1181 break;
1182 case AssemblyLevel::PARTIAL:
1183 ext = new PAMixedBilinearFormExtension(this);
1184 break;
1185 case AssemblyLevel::NONE:
1186 mfem_error("Matrix-free action not supported yet... stay tuned!");
1187 // ext = new MFMixedBilinearFormExtension(this);
1188 break;
1189 default:
1190 mfem_error("Unknown assembly level");
1191 }
1192 }
1193
Elem(int i,int j)1194 double & MixedBilinearForm::Elem (int i, int j)
1195 {
1196 return (*mat)(i, j);
1197 }
1198
Elem(int i,int j) const1199 const double & MixedBilinearForm::Elem (int i, int j) const
1200 {
1201 return (*mat)(i, j);
1202 }
1203
Mult(const Vector & x,Vector & y) const1204 void MixedBilinearForm::Mult(const Vector & x, Vector & y) const
1205 {
1206 y = 0.0;
1207 AddMult(x, y);
1208 }
1209
AddMult(const Vector & x,Vector & y,const double a) const1210 void MixedBilinearForm::AddMult(const Vector & x, Vector & y,
1211 const double a) const
1212 {
1213 if (ext)
1214 {
1215 ext->AddMult(x, y, a);
1216 }
1217 else
1218 {
1219 mat->AddMult(x, y, a);
1220 }
1221 }
1222
MultTranspose(const Vector & x,Vector & y) const1223 void MixedBilinearForm::MultTranspose(const Vector & x, Vector & y) const
1224 {
1225 y = 0.0;
1226 AddMultTranspose(x, y);
1227 }
1228
AddMultTranspose(const Vector & x,Vector & y,const double a) const1229 void MixedBilinearForm::AddMultTranspose(const Vector & x, Vector & y,
1230 const double a) const
1231 {
1232 if (ext)
1233 {
1234 ext->AddMultTranspose(x, y, a);
1235 }
1236 else
1237 {
1238 mat->AddMultTranspose(x, y, a);
1239 }
1240 }
1241
Inverse() const1242 MatrixInverse * MixedBilinearForm::Inverse() const
1243 {
1244 if (assembly != AssemblyLevel::LEGACY)
1245 {
1246 MFEM_WARNING("MixedBilinearForm::Inverse not possible with this "
1247 "assembly level!");
1248 return NULL;
1249 }
1250 else
1251 {
1252 return mat -> Inverse ();
1253 }
1254 }
1255
Finalize(int skip_zeros)1256 void MixedBilinearForm::Finalize (int skip_zeros)
1257 {
1258 if (assembly == AssemblyLevel::LEGACY)
1259 {
1260 mat -> Finalize (skip_zeros);
1261 }
1262 }
1263
GetBlocks(Array2D<SparseMatrix * > & blocks) const1264 void MixedBilinearForm::GetBlocks(Array2D<SparseMatrix *> &blocks) const
1265 {
1266 MFEM_VERIFY(trial_fes->GetOrdering() == Ordering::byNODES &&
1267 test_fes->GetOrdering() == Ordering::byNODES,
1268 "MixedBilinearForm::GetBlocks: both trial and test spaces "
1269 "must use Ordering::byNODES!");
1270
1271 blocks.SetSize(test_fes->GetVDim(), trial_fes->GetVDim());
1272
1273 mat->GetBlocks(blocks);
1274 }
1275
AddDomainIntegrator(BilinearFormIntegrator * bfi)1276 void MixedBilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi)
1277 {
1278 domain_integs.Append (bfi);
1279 }
1280
AddBoundaryIntegrator(BilinearFormIntegrator * bfi)1281 void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
1282 {
1283 boundary_integs.Append (bfi);
1284 boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
1285 }
1286
AddBoundaryIntegrator(BilinearFormIntegrator * bfi,Array<int> & bdr_marker)1287 void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
1288 Array<int> &bdr_marker)
1289 {
1290 boundary_integs.Append (bfi);
1291 boundary_integs_marker.Append(&bdr_marker);
1292 }
1293
AddTraceFaceIntegrator(BilinearFormIntegrator * bfi)1294 void MixedBilinearForm::AddTraceFaceIntegrator (BilinearFormIntegrator * bfi)
1295 {
1296 trace_face_integs.Append (bfi);
1297 }
1298
AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi)1299 void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi)
1300 {
1301 boundary_trace_face_integs.Append(bfi);
1302 // NULL marker means apply everywhere
1303 boundary_trace_face_integs_marker.Append(NULL);
1304 }
1305
AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi,Array<int> & bdr_marker)1306 void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi,
1307 Array<int> &bdr_marker)
1308 {
1309 boundary_trace_face_integs.Append(bfi);
1310 boundary_trace_face_integs_marker.Append(&bdr_marker);
1311 }
1312
Assemble(int skip_zeros)1313 void MixedBilinearForm::Assemble (int skip_zeros)
1314 {
1315 if (ext)
1316 {
1317 ext->Assemble();
1318 return;
1319 }
1320
1321 Array<int> tr_vdofs, te_vdofs;
1322 ElementTransformation *eltrans;
1323 DenseMatrix elemmat;
1324
1325 Mesh *mesh = test_fes -> GetMesh();
1326
1327 if (mat == NULL)
1328 {
1329 mat = new SparseMatrix(height, width);
1330 }
1331
1332 if (domain_integs.Size())
1333 {
1334 for (int i = 0; i < test_fes -> GetNE(); i++)
1335 {
1336 trial_fes -> GetElementVDofs (i, tr_vdofs);
1337 test_fes -> GetElementVDofs (i, te_vdofs);
1338 eltrans = test_fes -> GetElementTransformation (i);
1339 for (int k = 0; k < domain_integs.Size(); k++)
1340 {
1341 domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
1342 *test_fes -> GetFE(i),
1343 *eltrans, elemmat);
1344 mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1345 }
1346 }
1347 }
1348
1349 if (boundary_integs.Size())
1350 {
1351 // Which boundary attributes need to be processed?
1352 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1353 mesh->bdr_attributes.Max() : 0);
1354 bdr_attr_marker = 0;
1355 for (int k = 0; k < boundary_integs.Size(); k++)
1356 {
1357 if (boundary_integs_marker[k] == NULL)
1358 {
1359 bdr_attr_marker = 1;
1360 break;
1361 }
1362 Array<int> &bdr_marker = *boundary_integs_marker[k];
1363 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1364 "invalid boundary marker for boundary integrator #"
1365 << k << ", counting from zero");
1366 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1367 {
1368 bdr_attr_marker[i] |= bdr_marker[i];
1369 }
1370 }
1371
1372 for (int i = 0; i < test_fes -> GetNBE(); i++)
1373 {
1374 const int bdr_attr = mesh->GetBdrAttribute(i);
1375 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1376
1377 trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1378 test_fes -> GetBdrElementVDofs (i, te_vdofs);
1379 eltrans = test_fes -> GetBdrElementTransformation (i);
1380 for (int k = 0; k < boundary_integs.Size(); k++)
1381 {
1382 if (boundary_integs_marker[k] &&
1383 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
1384
1385 boundary_integs[k]->AssembleElementMatrix2 (*trial_fes -> GetBE(i),
1386 *test_fes -> GetBE(i),
1387 *eltrans, elemmat);
1388 mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1389 }
1390 }
1391 }
1392
1393 if (trace_face_integs.Size())
1394 {
1395 FaceElementTransformations *ftr;
1396 Array<int> te_vdofs2;
1397 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1398
1399 int nfaces = mesh->GetNumFaces();
1400 for (int i = 0; i < nfaces; i++)
1401 {
1402 ftr = mesh->GetFaceElementTransformations(i);
1403 trial_fes->GetFaceVDofs(i, tr_vdofs);
1404 test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1405 trial_face_fe = trial_fes->GetFaceElement(i);
1406 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1407 if (ftr->Elem2No >= 0)
1408 {
1409 test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
1410 te_vdofs.Append(te_vdofs2);
1411 test_fe2 = test_fes->GetFE(ftr->Elem2No);
1412 }
1413 else
1414 {
1415 // The test_fe2 object is really a dummy and not used on the
1416 // boundaries, but we can't dereference a NULL pointer, and we don't
1417 // want to actually make a fake element.
1418 test_fe2 = test_fe1;
1419 }
1420 for (int k = 0; k < trace_face_integs.Size(); k++)
1421 {
1422 trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
1423 *test_fe2, *ftr, elemmat);
1424 mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
1425 }
1426 }
1427 }
1428
1429 if (boundary_trace_face_integs.Size())
1430 {
1431 FaceElementTransformations *ftr;
1432 Array<int> te_vdofs2;
1433 const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
1434
1435 // Which boundary attributes need to be processed?
1436 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1437 mesh->bdr_attributes.Max() : 0);
1438 bdr_attr_marker = 0;
1439 for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
1440 {
1441 if (boundary_trace_face_integs_marker[k] == NULL)
1442 {
1443 bdr_attr_marker = 1;
1444 break;
1445 }
1446 Array<int> &bdr_marker = *boundary_trace_face_integs_marker[k];
1447 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1448 "invalid boundary marker for boundary trace face"
1449 "integrator #" << k << ", counting from zero");
1450 for (int i = 0; i < bdr_attr_marker.Size(); i++)
1451 {
1452 bdr_attr_marker[i] |= bdr_marker[i];
1453 }
1454 }
1455
1456 for (int i = 0; i < trial_fes -> GetNBE(); i++)
1457 {
1458 const int bdr_attr = mesh->GetBdrAttribute(i);
1459 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1460
1461 ftr = mesh->GetBdrFaceTransformations(i);
1462 if (ftr)
1463 {
1464 trial_fes->GetFaceVDofs(ftr->ElementNo, tr_vdofs);
1465 test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
1466 trial_face_fe = trial_fes->GetFaceElement(ftr->ElementNo);
1467 test_fe1 = test_fes->GetFE(ftr->Elem1No);
1468 // The test_fe2 object is really a dummy and not used on the
1469 // boundaries, but we can't dereference a NULL pointer, and we don't
1470 // want to actually make a fake element.
1471 test_fe2 = test_fe1;
1472 for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
1473 {
1474 if (boundary_trace_face_integs_marker[k] &&
1475 (*boundary_trace_face_integs_marker[k])[bdr_attr-1] == 0)
1476 { continue; }
1477
1478 boundary_trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe,
1479 *test_fe1,
1480 *test_fe2,
1481 *ftr, elemmat);
1482 mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
1483 }
1484 }
1485 }
1486 }
1487 }
1488
AssembleDiagonal_ADAt(const Vector & D,Vector & diag) const1489 void MixedBilinearForm::AssembleDiagonal_ADAt(const Vector &D,
1490 Vector &diag) const
1491 {
1492 if (ext)
1493 {
1494 MFEM_ASSERT(diag.Size() == test_fes->GetTrueVSize(),
1495 "Vector for holding diagonal has wrong size!");
1496 MFEM_ASSERT(D.Size() == trial_fes->GetTrueVSize(),
1497 "Vector for holding diagonal has wrong size!");
1498 const Operator *P_trial = trial_fes->GetProlongationMatrix();
1499 const Operator *P_test = test_fes->GetProlongationMatrix();
1500 if (!IsIdentityProlongation(P_trial))
1501 {
1502 Vector local_D(P_trial->Height());
1503 P_trial->Mult(D, local_D);
1504
1505 if (!IsIdentityProlongation(P_test))
1506 {
1507 Vector local_diag(P_test->Height());
1508 ext->AssembleDiagonal_ADAt(local_D, local_diag);
1509 P_test->MultTranspose(local_diag, diag);
1510 }
1511 else
1512 {
1513 ext->AssembleDiagonal_ADAt(local_D, diag);
1514 }
1515 }
1516 else
1517 {
1518 if (!IsIdentityProlongation(P_test))
1519 {
1520 Vector local_diag(P_test->Height());
1521 ext->AssembleDiagonal_ADAt(D, local_diag);
1522 P_test->MultTranspose(local_diag, diag);
1523 }
1524 else
1525 {
1526 ext->AssembleDiagonal_ADAt(D, diag);
1527 }
1528 }
1529 }
1530 else
1531 {
1532 MFEM_ABORT("Not implemented. Maybe assemble your bilinear form into a "
1533 "matrix and use SparseMatrix functions?");
1534 }
1535 }
1536
ConformingAssemble()1537 void MixedBilinearForm::ConformingAssemble()
1538 {
1539 if (assembly != AssemblyLevel::LEGACY)
1540 {
1541 MFEM_WARNING("Conforming assemble not supported for this assembly level!");
1542 return;
1543 }
1544
1545 Finalize();
1546
1547 const SparseMatrix *P2 = test_fes->GetConformingProlongation();
1548 if (P2)
1549 {
1550 SparseMatrix *R = Transpose(*P2);
1551 SparseMatrix *RA = mfem::Mult(*R, *mat);
1552 delete R;
1553 delete mat;
1554 mat = RA;
1555 }
1556
1557 const SparseMatrix *P1 = trial_fes->GetConformingProlongation();
1558 if (P1)
1559 {
1560 SparseMatrix *RAP = mfem::Mult(*mat, *P1);
1561 delete mat;
1562 mat = RAP;
1563 }
1564
1565 height = mat->Height();
1566 width = mat->Width();
1567 }
1568
1569
ComputeElementMatrix(int i,DenseMatrix & elmat)1570 void MixedBilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
1571 {
1572 if (domain_integs.Size())
1573 {
1574 const FiniteElement &trial_fe = *trial_fes->GetFE(i);
1575 const FiniteElement &test_fe = *test_fes->GetFE(i);
1576 ElementTransformation *eltrans = test_fes->GetElementTransformation(i);
1577 domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1578 elmat);
1579 for (int k = 1; k < domain_integs.Size(); k++)
1580 {
1581 domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1582 elemmat);
1583 elmat += elemmat;
1584 }
1585 }
1586 else
1587 {
1588 trial_fes->GetElementVDofs(i, trial_vdofs);
1589 test_fes->GetElementVDofs(i, test_vdofs);
1590 elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
1591 elmat = 0.0;
1592 }
1593 }
1594
ComputeBdrElementMatrix(int i,DenseMatrix & elmat)1595 void MixedBilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
1596 {
1597 if (boundary_integs.Size())
1598 {
1599 const FiniteElement &trial_be = *trial_fes->GetBE(i);
1600 const FiniteElement &test_be = *test_fes->GetBE(i);
1601 ElementTransformation *eltrans = test_fes->GetBdrElementTransformation(i);
1602 boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1603 elmat);
1604 for (int k = 1; k < boundary_integs.Size(); k++)
1605 {
1606 boundary_integs[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1607 elemmat);
1608 elmat += elemmat;
1609 }
1610 }
1611 else
1612 {
1613 trial_fes->GetBdrElementVDofs(i, trial_vdofs);
1614 test_fes->GetBdrElementVDofs(i, test_vdofs);
1615 elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
1616 elmat = 0.0;
1617 }
1618 }
1619
AssembleElementMatrix(int i,const DenseMatrix & elmat,int skip_zeros)1620 void MixedBilinearForm::AssembleElementMatrix(
1621 int i, const DenseMatrix &elmat, int skip_zeros)
1622 {
1623 AssembleElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1624 }
1625
AssembleElementMatrix(int i,const DenseMatrix & elmat,Array<int> & trial_vdofs,Array<int> & test_vdofs,int skip_zeros)1626 void MixedBilinearForm::AssembleElementMatrix(
1627 int i, const DenseMatrix &elmat, Array<int> &trial_vdofs,
1628 Array<int> &test_vdofs, int skip_zeros)
1629 {
1630 trial_fes->GetElementVDofs(i, trial_vdofs);
1631 test_fes->GetElementVDofs(i, test_vdofs);
1632 if (mat == NULL)
1633 {
1634 mat = new SparseMatrix(height, width);
1635 }
1636 mat->AddSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
1637 }
1638
AssembleBdrElementMatrix(int i,const DenseMatrix & elmat,int skip_zeros)1639 void MixedBilinearForm::AssembleBdrElementMatrix(
1640 int i, const DenseMatrix &elmat, int skip_zeros)
1641 {
1642 AssembleBdrElementMatrix(i, elmat, trial_vdofs, test_vdofs, skip_zeros);
1643 }
1644
AssembleBdrElementMatrix(int i,const DenseMatrix & elmat,Array<int> & trial_vdofs,Array<int> & test_vdofs,int skip_zeros)1645 void MixedBilinearForm::AssembleBdrElementMatrix(
1646 int i, const DenseMatrix &elmat, Array<int> &trial_vdofs,
1647 Array<int> &test_vdofs, int skip_zeros)
1648 {
1649 trial_fes->GetBdrElementVDofs(i, trial_vdofs);
1650 test_fes->GetBdrElementVDofs(i, test_vdofs);
1651 if (mat == NULL)
1652 {
1653 mat = new SparseMatrix(height, width);
1654 }
1655 mat->AddSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
1656 }
1657
EliminateTrialDofs(const Array<int> & bdr_attr_is_ess,const Vector & sol,Vector & rhs)1658 void MixedBilinearForm::EliminateTrialDofs (
1659 const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
1660 {
1661 int i, j, k;
1662 Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
1663
1664 cols_marker = 0;
1665 for (i = 0; i < trial_fes -> GetNBE(); i++)
1666 if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
1667 {
1668 trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1669 for (j = 0; j < tr_vdofs.Size(); j++)
1670 {
1671 if ( (k = tr_vdofs[j]) < 0 )
1672 {
1673 k = -1-k;
1674 }
1675 cols_marker[k] = 1;
1676 }
1677 }
1678 mat -> EliminateCols (cols_marker, &sol, &rhs);
1679 }
1680
EliminateEssentialBCFromTrialDofs(const Array<int> & marked_vdofs,const Vector & sol,Vector & rhs)1681 void MixedBilinearForm::EliminateEssentialBCFromTrialDofs (
1682 const Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
1683 {
1684 mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1685 }
1686
EliminateTestDofs(const Array<int> & bdr_attr_is_ess)1687 void MixedBilinearForm::EliminateTestDofs (const Array<int> &bdr_attr_is_ess)
1688 {
1689 int i, j, k;
1690 Array<int> te_vdofs;
1691
1692 for (i = 0; i < test_fes -> GetNBE(); i++)
1693 if (bdr_attr_is_ess[test_fes -> GetBdrAttribute (i)-1])
1694 {
1695 test_fes -> GetBdrElementVDofs (i, te_vdofs);
1696 for (j = 0; j < te_vdofs.Size(); j++)
1697 {
1698 if ( (k = te_vdofs[j]) < 0 )
1699 {
1700 k = -1-k;
1701 }
1702 mat -> EliminateRow (k);
1703 }
1704 }
1705 }
1706
FormRectangularSystemMatrix(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,OperatorHandle & A)1707 void MixedBilinearForm::FormRectangularSystemMatrix(
1708 const Array<int> &trial_tdof_list,
1709 const Array<int> &test_tdof_list,
1710 OperatorHandle &A)
1711
1712 {
1713 if (ext)
1714 {
1715 ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
1716 return;
1717 }
1718
1719 const SparseMatrix *test_P = test_fes->GetConformingProlongation();
1720 const SparseMatrix *trial_P = trial_fes->GetConformingProlongation();
1721
1722 mat->Finalize();
1723
1724 if (test_P) // TODO: Must actually check for trial_P too
1725 {
1726 SparseMatrix *m = RAP(*test_P, *mat, *trial_P);
1727 delete mat;
1728 mat = m;
1729 }
1730
1731 Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
1732 FiniteElementSpace::ListToMarker(trial_tdof_list, trial_fes->GetTrueVSize(),
1733 ess_trial_tdof_marker);
1734 FiniteElementSpace::ListToMarker(test_tdof_list, test_fes->GetTrueVSize(),
1735 ess_test_tdof_marker);
1736
1737 mat_e = new SparseMatrix(mat->Height(), mat->Width());
1738 mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
1739
1740 for (int i=0; i<test_tdof_list.Size(); ++i)
1741 {
1742 mat->EliminateRow(test_tdof_list[i]);
1743 }
1744 mat_e->Finalize();
1745 A.Reset(mat, false);
1746 }
1747
FormRectangularLinearSystem(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B)1748 void MixedBilinearForm::FormRectangularLinearSystem(
1749 const Array<int> &trial_tdof_list,
1750 const Array<int> &test_tdof_list,
1751 Vector &x, Vector &b,
1752 OperatorHandle &A,
1753 Vector &X, Vector &B)
1754 {
1755 if (ext)
1756 {
1757 ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
1758 x, b, A, X, B);
1759 return;
1760 }
1761
1762 const Operator *Pi = this->GetProlongation();
1763 const Operator *Po = this->GetOutputProlongation();
1764 const Operator *Ri = this->GetRestriction();
1765 InitTVectors(Po, Ri, Pi, x, b, X, B);
1766
1767 if (!mat_e)
1768 {
1769 FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list,
1770 A); // Set A = mat_e
1771 }
1772 // Eliminate essential BCs with B -= Ab xb
1773 mat_e->AddMult(X, B, -1.0);
1774
1775 B.SetSubVector(test_tdof_list, 0.0);
1776 }
1777
Update()1778 void MixedBilinearForm::Update()
1779 {
1780 delete mat;
1781 mat = NULL;
1782 delete mat_e;
1783 mat_e = NULL;
1784 height = test_fes->GetVSize();
1785 width = trial_fes->GetVSize();
1786 if (ext) { ext->Update(); }
1787 }
1788
~MixedBilinearForm()1789 MixedBilinearForm::~MixedBilinearForm()
1790 {
1791 if (mat) { delete mat; }
1792 if (mat_e) { delete mat_e; }
1793 if (!extern_bfs)
1794 {
1795 int i;
1796 for (i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
1797 for (i = 0; i < boundary_integs.Size(); i++)
1798 { delete boundary_integs[i]; }
1799 for (i = 0; i < trace_face_integs.Size(); i++)
1800 { delete trace_face_integs[i]; }
1801 for (i = 0; i < boundary_trace_face_integs.Size(); i++)
1802 { delete boundary_trace_face_integs[i]; }
1803 }
1804 delete ext;
1805 }
1806
SetAssemblyLevel(AssemblyLevel assembly_level)1807 void DiscreteLinearOperator::SetAssemblyLevel(AssemblyLevel assembly_level)
1808 {
1809 if (ext)
1810 {
1811 MFEM_ABORT("the assembly level has already been set!");
1812 }
1813 assembly = assembly_level;
1814 switch (assembly)
1815 {
1816 case AssemblyLevel::LEGACY:
1817 case AssemblyLevel::FULL:
1818 // Use the original implementation for now
1819 break;
1820 case AssemblyLevel::ELEMENT:
1821 mfem_error("Element assembly not supported yet... stay tuned!");
1822 break;
1823 case AssemblyLevel::PARTIAL:
1824 ext = new PADiscreteLinearOperatorExtension(this);
1825 break;
1826 case AssemblyLevel::NONE:
1827 mfem_error("Matrix-free action not supported yet... stay tuned!");
1828 break;
1829 default:
1830 mfem_error("Unknown assembly level");
1831 }
1832 }
1833
Assemble(int skip_zeros)1834 void DiscreteLinearOperator::Assemble(int skip_zeros)
1835 {
1836 if (ext)
1837 {
1838 ext->Assemble();
1839 return;
1840 }
1841
1842 Array<int> dom_vdofs, ran_vdofs;
1843 ElementTransformation *T;
1844 const FiniteElement *dom_fe, *ran_fe;
1845 DenseMatrix totelmat, elmat;
1846
1847 if (mat == NULL)
1848 {
1849 mat = new SparseMatrix(height, width);
1850 }
1851
1852 if (domain_integs.Size() > 0)
1853 {
1854 for (int i = 0; i < test_fes->GetNE(); i++)
1855 {
1856 trial_fes->GetElementVDofs(i, dom_vdofs);
1857 test_fes->GetElementVDofs(i, ran_vdofs);
1858 T = test_fes->GetElementTransformation(i);
1859 dom_fe = trial_fes->GetFE(i);
1860 ran_fe = test_fes->GetFE(i);
1861
1862 domain_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
1863 totelmat);
1864 for (int j = 1; j < domain_integs.Size(); j++)
1865 {
1866 domain_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
1867 elmat);
1868 totelmat += elmat;
1869 }
1870 mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1871 }
1872 }
1873
1874 if (trace_face_integs.Size())
1875 {
1876 const int nfaces = test_fes->GetMesh()->GetNumFaces();
1877 for (int i = 0; i < nfaces; i++)
1878 {
1879 trial_fes->GetFaceVDofs(i, dom_vdofs);
1880 test_fes->GetFaceVDofs(i, ran_vdofs);
1881 T = test_fes->GetMesh()->GetFaceTransformation(i);
1882 dom_fe = trial_fes->GetFaceElement(i);
1883 ran_fe = test_fes->GetFaceElement(i);
1884
1885 trace_face_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
1886 totelmat);
1887 for (int j = 1; j < trace_face_integs.Size(); j++)
1888 {
1889 trace_face_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
1890 elmat);
1891 totelmat += elmat;
1892 }
1893 mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
1894 }
1895 }
1896 }
1897
1898 }
1899