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 // Implementations of classes FABilinearFormExtension, EABilinearFormExtension,
13 // PABilinearFormExtension and MFBilinearFormExtension.
14
15 #include "../general/forall.hpp"
16 #include "bilinearform.hpp"
17 #include "pbilinearform.hpp"
18 #include "pgridfunc.hpp"
19 #include "ceed/util.hpp"
20
21 namespace mfem
22 {
23
BilinearFormExtension(BilinearForm * form)24 BilinearFormExtension::BilinearFormExtension(BilinearForm *form)
25 : Operator(form->Size()), a(form)
26 {
27 // empty
28 }
29
GetProlongation() const30 const Operator *BilinearFormExtension::GetProlongation() const
31 {
32 return a->GetProlongation();
33 }
34
GetRestriction() const35 const Operator *BilinearFormExtension::GetRestriction() const
36 {
37 return a->GetRestriction();
38 }
39
40 // Data and methods for partially-assembled bilinear forms
MFBilinearFormExtension(BilinearForm * form)41 MFBilinearFormExtension::MFBilinearFormExtension(BilinearForm *form)
42 : BilinearFormExtension(form),
43 trialFes(a->FESpace()),
44 testFes(a->FESpace())
45 {
46 elem_restrict = NULL;
47 int_face_restrict_lex = NULL;
48 bdr_face_restrict_lex = NULL;
49 }
50
Assemble()51 void MFBilinearFormExtension::Assemble()
52 {
53 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
54 const int integratorCount = integrators.Size();
55 for (int i = 0; i < integratorCount; ++i)
56 {
57 integrators[i]->AssembleMF(*a->FESpace());
58 }
59 }
60
AssembleDiagonal(Vector & y) const61 void MFBilinearFormExtension::AssembleDiagonal(Vector &y) const
62 {
63 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
64
65 const int iSz = integrators.Size();
66 if (elem_restrict && !DeviceCanUseCeed())
67 {
68 localY = 0.0;
69 for (int i = 0; i < iSz; ++i)
70 {
71 integrators[i]->AssembleDiagonalMF(localY);
72 }
73 const ElementRestriction* H1elem_restrict =
74 dynamic_cast<const ElementRestriction*>(elem_restrict);
75 if (H1elem_restrict)
76 {
77 H1elem_restrict->MultTransposeUnsigned(localY, y);
78 }
79 else
80 {
81 elem_restrict->MultTranspose(localY, y);
82 }
83 }
84 else
85 {
86 y.UseDevice(true); // typically this is a large vector, so store on device
87 y = 0.0;
88 for (int i = 0; i < iSz; ++i)
89 {
90 integrators[i]->AssembleDiagonalMF(y);
91 }
92 }
93 }
94
Update()95 void MFBilinearFormExtension::Update()
96 {
97 FiniteElementSpace *fes = a->FESpace();
98 height = width = fes->GetVSize();
99 trialFes = fes;
100 testFes = fes;
101
102 elem_restrict = nullptr;
103 int_face_restrict_lex = nullptr;
104 bdr_face_restrict_lex = nullptr;
105 }
106
FormSystemMatrix(const Array<int> & ess_tdof_list,OperatorHandle & A)107 void MFBilinearFormExtension::FormSystemMatrix(const Array<int> &ess_tdof_list,
108 OperatorHandle &A)
109 {
110 Operator *oper;
111 Operator::FormSystemOperator(ess_tdof_list, oper);
112 A.Reset(oper); // A will own oper
113 }
114
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B,int copy_interior)115 void MFBilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
116 Vector &x, Vector &b,
117 OperatorHandle &A,
118 Vector &X, Vector &B,
119 int copy_interior)
120 {
121 Operator *oper;
122 Operator::FormLinearSystem(ess_tdof_list, x, b, oper, X, B, copy_interior);
123 A.Reset(oper); // A will own oper
124 }
125
Mult(const Vector & x,Vector & y) const126 void MFBilinearFormExtension::Mult(const Vector &x, Vector &y) const
127 {
128 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
129
130 const int iSz = integrators.Size();
131 if (DeviceCanUseCeed() || !elem_restrict)
132 {
133 y.UseDevice(true); // typically this is a large vector, so store on device
134 y = 0.0;
135 for (int i = 0; i < iSz; ++i)
136 {
137 integrators[i]->AddMultMF(x, y);
138 }
139 }
140 else
141 {
142 elem_restrict->Mult(x, localX);
143 localY = 0.0;
144 for (int i = 0; i < iSz; ++i)
145 {
146 integrators[i]->AddMultMF(localX, localY);
147 }
148 elem_restrict->MultTranspose(localY, y);
149 }
150
151 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
152 const int iFISz = intFaceIntegrators.Size();
153 if (int_face_restrict_lex && iFISz>0)
154 {
155 int_face_restrict_lex->Mult(x, faceIntX);
156 if (faceIntX.Size()>0)
157 {
158 faceIntY = 0.0;
159 for (int i = 0; i < iFISz; ++i)
160 {
161 intFaceIntegrators[i]->AddMultMF(faceIntX, faceIntY);
162 }
163 int_face_restrict_lex->AddMultTranspose(faceIntY, y);
164 }
165 }
166
167 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
168 const int bFISz = bdrFaceIntegrators.Size();
169 if (bdr_face_restrict_lex && bFISz>0)
170 {
171 bdr_face_restrict_lex->Mult(x, faceBdrX);
172 if (faceBdrX.Size()>0)
173 {
174 faceBdrY = 0.0;
175 for (int i = 0; i < bFISz; ++i)
176 {
177 bdrFaceIntegrators[i]->AddMultMF(faceBdrX, faceBdrY);
178 }
179 bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
180 }
181 }
182 }
183
MultTranspose(const Vector & x,Vector & y) const184 void MFBilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
185 {
186 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
187 const int iSz = integrators.Size();
188 if (elem_restrict)
189 {
190 elem_restrict->Mult(x, localX);
191 localY = 0.0;
192 for (int i = 0; i < iSz; ++i)
193 {
194 integrators[i]->AddMultTransposeMF(localX, localY);
195 }
196 elem_restrict->MultTranspose(localY, y);
197 }
198 else
199 {
200 y.UseDevice(true);
201 y = 0.0;
202 for (int i = 0; i < iSz; ++i)
203 {
204 integrators[i]->AddMultTransposeMF(x, y);
205 }
206 }
207
208 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
209 const int iFISz = intFaceIntegrators.Size();
210 if (int_face_restrict_lex && iFISz>0)
211 {
212 int_face_restrict_lex->Mult(x, faceIntX);
213 if (faceIntX.Size()>0)
214 {
215 faceIntY = 0.0;
216 for (int i = 0; i < iFISz; ++i)
217 {
218 intFaceIntegrators[i]->AddMultTransposeMF(faceIntX, faceIntY);
219 }
220 int_face_restrict_lex->AddMultTranspose(faceIntY, y);
221 }
222 }
223
224 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
225 const int bFISz = bdrFaceIntegrators.Size();
226 if (bdr_face_restrict_lex && bFISz>0)
227 {
228 bdr_face_restrict_lex->Mult(x, faceBdrX);
229 if (faceBdrX.Size()>0)
230 {
231 faceBdrY = 0.0;
232 for (int i = 0; i < bFISz; ++i)
233 {
234 bdrFaceIntegrators[i]->AddMultTransposeMF(faceBdrX, faceBdrY);
235 }
236 bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
237 }
238 }
239 }
240
241 // Data and methods for partially-assembled bilinear forms
PABilinearFormExtension(BilinearForm * form)242 PABilinearFormExtension::PABilinearFormExtension(BilinearForm *form)
243 : BilinearFormExtension(form),
244 trialFes(a->FESpace()),
245 testFes(a->FESpace())
246 {
247 elem_restrict = NULL;
248 int_face_restrict_lex = NULL;
249 bdr_face_restrict_lex = NULL;
250 }
251
SetupRestrictionOperators(const L2FaceValues m)252 void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
253 {
254 ElementDofOrdering ordering = UsesTensorBasis(*a->FESpace())?
255 ElementDofOrdering::LEXICOGRAPHIC:
256 ElementDofOrdering::NATIVE;
257 elem_restrict = trialFes->GetElementRestriction(ordering);
258 if (elem_restrict)
259 {
260 localX.SetSize(elem_restrict->Height(), Device::GetDeviceMemoryType());
261 localY.SetSize(elem_restrict->Height(), Device::GetDeviceMemoryType());
262 localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
263 }
264
265 // Construct face restriction operators only if the bilinear form has
266 // interior or boundary face integrators
267 if (int_face_restrict_lex == NULL && a->GetFBFI()->Size() > 0)
268 {
269 int_face_restrict_lex = trialFes->GetFaceRestriction(
270 ElementDofOrdering::LEXICOGRAPHIC,
271 FaceType::Interior);
272 faceIntX.SetSize(int_face_restrict_lex->Height(), Device::GetMemoryType());
273 faceIntY.SetSize(int_face_restrict_lex->Height(), Device::GetMemoryType());
274 faceIntY.UseDevice(true); // ensure 'faceIntY = 0.0' is done on device
275 }
276
277 if (bdr_face_restrict_lex == NULL && a->GetBFBFI()->Size() > 0)
278 {
279 bdr_face_restrict_lex = trialFes->GetFaceRestriction(
280 ElementDofOrdering::LEXICOGRAPHIC,
281 FaceType::Boundary,
282 m);
283 faceBdrX.SetSize(bdr_face_restrict_lex->Height(), Device::GetMemoryType());
284 faceBdrY.SetSize(bdr_face_restrict_lex->Height(), Device::GetMemoryType());
285 faceBdrY.UseDevice(true); // ensure 'faceBoundY = 0.0' is done on device
286 }
287 }
288
Assemble()289 void PABilinearFormExtension::Assemble()
290 {
291 SetupRestrictionOperators(L2FaceValues::DoubleValued);
292
293 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
294 const int integratorCount = integrators.Size();
295 for (int i = 0; i < integratorCount; ++i)
296 {
297 integrators[i]->AssemblePA(*a->FESpace());
298 }
299
300 MFEM_VERIFY(a->GetBBFI()->Size() == 0,
301 "Partial assembly does not support AddBoundaryIntegrator yet.");
302
303 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
304 const int intFaceIntegratorCount = intFaceIntegrators.Size();
305 for (int i = 0; i < intFaceIntegratorCount; ++i)
306 {
307 intFaceIntegrators[i]->AssemblePAInteriorFaces(*a->FESpace());
308 }
309
310 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
311 const int boundFaceIntegratorCount = bdrFaceIntegrators.Size();
312 for (int i = 0; i < boundFaceIntegratorCount; ++i)
313 {
314 bdrFaceIntegrators[i]->AssemblePABoundaryFaces(*a->FESpace());
315 }
316 }
317
AssembleDiagonal(Vector & y) const318 void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
319 {
320 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
321
322 const int iSz = integrators.Size();
323 if (elem_restrict && !DeviceCanUseCeed())
324 {
325 localY = 0.0;
326 for (int i = 0; i < iSz; ++i)
327 {
328 integrators[i]->AssembleDiagonalPA(localY);
329 }
330 const ElementRestriction* H1elem_restrict =
331 dynamic_cast<const ElementRestriction*>(elem_restrict);
332 if (H1elem_restrict)
333 {
334 H1elem_restrict->MultTransposeUnsigned(localY, y);
335 }
336 else
337 {
338 elem_restrict->MultTranspose(localY, y);
339 }
340 }
341 else
342 {
343 y.UseDevice(true); // typically this is a large vector, so store on device
344 y = 0.0;
345 for (int i = 0; i < iSz; ++i)
346 {
347 integrators[i]->AssembleDiagonalPA(y);
348 }
349 }
350 }
351
Update()352 void PABilinearFormExtension::Update()
353 {
354 FiniteElementSpace *fes = a->FESpace();
355 height = width = fes->GetVSize();
356 trialFes = fes;
357 testFes = fes;
358
359 elem_restrict = nullptr;
360 int_face_restrict_lex = nullptr;
361 bdr_face_restrict_lex = nullptr;
362 }
363
FormSystemMatrix(const Array<int> & ess_tdof_list,OperatorHandle & A)364 void PABilinearFormExtension::FormSystemMatrix(const Array<int> &ess_tdof_list,
365 OperatorHandle &A)
366 {
367 Operator *oper;
368 Operator::FormSystemOperator(ess_tdof_list, oper);
369 A.Reset(oper); // A will own oper
370 }
371
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B,int copy_interior)372 void PABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
373 Vector &x, Vector &b,
374 OperatorHandle &A,
375 Vector &X, Vector &B,
376 int copy_interior)
377 {
378 Operator *oper;
379 Operator::FormLinearSystem(ess_tdof_list, x, b, oper, X, B, copy_interior);
380 A.Reset(oper); // A will own oper
381 }
382
Mult(const Vector & x,Vector & y) const383 void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
384 {
385 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
386
387 const int iSz = integrators.Size();
388 if (DeviceCanUseCeed() || !elem_restrict)
389 {
390 y.UseDevice(true); // typically this is a large vector, so store on device
391 y = 0.0;
392 for (int i = 0; i < iSz; ++i)
393 {
394 integrators[i]->AddMultPA(x, y);
395 }
396 }
397 else
398 {
399 elem_restrict->Mult(x, localX);
400 localY = 0.0;
401 for (int i = 0; i < iSz; ++i)
402 {
403 integrators[i]->AddMultPA(localX, localY);
404 }
405 elem_restrict->MultTranspose(localY, y);
406 }
407
408 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
409 const int iFISz = intFaceIntegrators.Size();
410 if (int_face_restrict_lex && iFISz>0)
411 {
412 int_face_restrict_lex->Mult(x, faceIntX);
413 if (faceIntX.Size()>0)
414 {
415 faceIntY = 0.0;
416 for (int i = 0; i < iFISz; ++i)
417 {
418 intFaceIntegrators[i]->AddMultPA(faceIntX, faceIntY);
419 }
420 int_face_restrict_lex->AddMultTranspose(faceIntY, y);
421 }
422 }
423
424 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
425 const int bFISz = bdrFaceIntegrators.Size();
426 if (bdr_face_restrict_lex && bFISz>0)
427 {
428 bdr_face_restrict_lex->Mult(x, faceBdrX);
429 if (faceBdrX.Size()>0)
430 {
431 faceBdrY = 0.0;
432 for (int i = 0; i < bFISz; ++i)
433 {
434 bdrFaceIntegrators[i]->AddMultPA(faceBdrX, faceBdrY);
435 }
436 bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
437 }
438 }
439 }
440
MultTranspose(const Vector & x,Vector & y) const441 void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
442 {
443 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
444 const int iSz = integrators.Size();
445 if (elem_restrict)
446 {
447 elem_restrict->Mult(x, localX);
448 localY = 0.0;
449 for (int i = 0; i < iSz; ++i)
450 {
451 integrators[i]->AddMultTransposePA(localX, localY);
452 }
453 elem_restrict->MultTranspose(localY, y);
454 }
455 else
456 {
457 y.UseDevice(true);
458 y = 0.0;
459 for (int i = 0; i < iSz; ++i)
460 {
461 integrators[i]->AddMultTransposePA(x, y);
462 }
463 }
464
465 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
466 const int iFISz = intFaceIntegrators.Size();
467 if (int_face_restrict_lex && iFISz>0)
468 {
469 int_face_restrict_lex->Mult(x, faceIntX);
470 if (faceIntX.Size()>0)
471 {
472 faceIntY = 0.0;
473 for (int i = 0; i < iFISz; ++i)
474 {
475 intFaceIntegrators[i]->AddMultTransposePA(faceIntX, faceIntY);
476 }
477 int_face_restrict_lex->AddMultTranspose(faceIntY, y);
478 }
479 }
480
481 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
482 const int bFISz = bdrFaceIntegrators.Size();
483 if (bdr_face_restrict_lex && bFISz>0)
484 {
485 bdr_face_restrict_lex->Mult(x, faceBdrX);
486 if (faceBdrX.Size()>0)
487 {
488 faceBdrY = 0.0;
489 for (int i = 0; i < bFISz; ++i)
490 {
491 bdrFaceIntegrators[i]->AddMultTransposePA(faceBdrX, faceBdrY);
492 }
493 bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
494 }
495 }
496 }
497
498 // Data and methods for element-assembled bilinear forms
EABilinearFormExtension(BilinearForm * form)499 EABilinearFormExtension::EABilinearFormExtension(BilinearForm *form)
500 : PABilinearFormExtension(form),
501 factorize_face_terms(form->FESpace()->IsDGSpace())
502 {
503 }
504
Assemble()505 void EABilinearFormExtension::Assemble()
506 {
507 SetupRestrictionOperators(L2FaceValues::SingleValued);
508
509 ne = trialFes->GetMesh()->GetNE();
510 elemDofs = trialFes->GetFE(0)->GetDof();
511
512 ea_data.SetSize(ne*elemDofs*elemDofs, Device::GetMemoryType());
513 ea_data.UseDevice(true);
514
515 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
516 const int integratorCount = integrators.Size();
517 for (int i = 0; i < integratorCount; ++i)
518 {
519 integrators[i]->AssembleEA(*a->FESpace(), ea_data, i);
520 }
521
522 faceDofs = trialFes ->
523 GetTraceElement(0, trialFes->GetMesh()->GetFaceBaseGeometry(0)) ->
524 GetDof();
525
526 MFEM_VERIFY(a->GetBBFI()->Size() == 0,
527 "Element assembly does not support AddBoundaryIntegrator yet.");
528
529 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
530 const int intFaceIntegratorCount = intFaceIntegrators.Size();
531 if (intFaceIntegratorCount>0)
532 {
533 nf_int = trialFes->GetNFbyType(FaceType::Interior);
534 ea_data_int.SetSize(2*nf_int*faceDofs*faceDofs, Device::GetMemoryType());
535 ea_data_ext.SetSize(2*nf_int*faceDofs*faceDofs, Device::GetMemoryType());
536 }
537 for (int i = 0; i < intFaceIntegratorCount; ++i)
538 {
539 intFaceIntegrators[i]->AssembleEAInteriorFaces(*a->FESpace(),
540 ea_data_int,
541 ea_data_ext,
542 i);
543 }
544
545 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
546 const int boundFaceIntegratorCount = bdrFaceIntegrators.Size();
547 if (boundFaceIntegratorCount>0)
548 {
549 nf_bdr = trialFes->GetNFbyType(FaceType::Boundary);
550 ea_data_bdr.SetSize(nf_bdr*faceDofs*faceDofs, Device::GetMemoryType());
551 ea_data_bdr = 0.0;
552 }
553 for (int i = 0; i < boundFaceIntegratorCount; ++i)
554 {
555 bdrFaceIntegrators[i]->AssembleEABoundaryFaces(*a->FESpace(),ea_data_bdr,i);
556 }
557
558 if (factorize_face_terms && int_face_restrict_lex)
559 {
560 auto restFint = dynamic_cast<const L2FaceRestriction&>(*int_face_restrict_lex);
561 restFint.AddFaceMatricesToElementMatrices(ea_data_int, ea_data);
562 }
563 if (factorize_face_terms && bdr_face_restrict_lex)
564 {
565 auto restFbdr = dynamic_cast<const L2FaceRestriction&>(*bdr_face_restrict_lex);
566 restFbdr.AddFaceMatricesToElementMatrices(ea_data_bdr, ea_data);
567 }
568 }
569
Mult(const Vector & x,Vector & y) const570 void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
571 {
572 // Apply the Element Restriction
573 const bool useRestrict = !DeviceCanUseCeed() && elem_restrict;
574 if (!useRestrict)
575 {
576 y.UseDevice(true); // typically this is a large vector, so store on device
577 y = 0.0;
578 }
579 else
580 {
581 elem_restrict->Mult(x, localX);
582 localY = 0.0;
583 }
584 // Apply the Element Matrices
585 const int NDOFS = elemDofs;
586 auto X = Reshape(useRestrict?localX.Read():x.Read(), NDOFS, ne);
587 auto Y = Reshape(useRestrict?localY.ReadWrite():y.ReadWrite(), NDOFS, ne);
588 auto A = Reshape(ea_data.Read(), NDOFS, NDOFS, ne);
589 MFEM_FORALL(glob_j, ne*NDOFS,
590 {
591 const int e = glob_j/NDOFS;
592 const int j = glob_j%NDOFS;
593 double res = 0.0;
594 for (int i = 0; i < NDOFS; i++)
595 {
596 res += A(i, j, e)*X(i, e);
597 }
598 Y(j, e) += res;
599 });
600 // Apply the Element Restriction transposed
601 if (useRestrict)
602 {
603 elem_restrict->MultTranspose(localY, y);
604 }
605
606 // Treatment of interior faces
607 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
608 const int iFISz = intFaceIntegrators.Size();
609 if (int_face_restrict_lex && iFISz>0)
610 {
611 // Apply the Interior Face Restriction
612 int_face_restrict_lex->Mult(x, faceIntX);
613 if (faceIntX.Size()>0)
614 {
615 faceIntY = 0.0;
616 // Apply the interior face matrices
617 const int NDOFS = faceDofs;
618 auto X = Reshape(faceIntX.Read(), NDOFS, 2, nf_int);
619 auto Y = Reshape(faceIntY.ReadWrite(), NDOFS, 2, nf_int);
620 if (!factorize_face_terms)
621 {
622 auto A_int = Reshape(ea_data_int.Read(), NDOFS, NDOFS, 2, nf_int);
623 MFEM_FORALL(glob_j, nf_int*NDOFS,
624 {
625 const int f = glob_j/NDOFS;
626 const int j = glob_j%NDOFS;
627 double res = 0.0;
628 for (int i = 0; i < NDOFS; i++)
629 {
630 res += A_int(i, j, 0, f)*X(i, 0, f);
631 }
632 Y(j, 0, f) += res;
633 res = 0.0;
634 for (int i = 0; i < NDOFS; i++)
635 {
636 res += A_int(i, j, 1, f)*X(i, 1, f);
637 }
638 Y(j, 1, f) += res;
639 });
640 }
641 auto A_ext = Reshape(ea_data_ext.Read(), NDOFS, NDOFS, 2, nf_int);
642 MFEM_FORALL(glob_j, nf_int*NDOFS,
643 {
644 const int f = glob_j/NDOFS;
645 const int j = glob_j%NDOFS;
646 double res = 0.0;
647 for (int i = 0; i < NDOFS; i++)
648 {
649 res += A_ext(i, j, 0, f)*X(i, 0, f);
650 }
651 Y(j, 1, f) += res;
652 res = 0.0;
653 for (int i = 0; i < NDOFS; i++)
654 {
655 res += A_ext(i, j, 1, f)*X(i, 1, f);
656 }
657 Y(j, 0, f) += res;
658 });
659 // Apply the Interior Face Restriction transposed
660 int_face_restrict_lex->AddMultTranspose(faceIntY, y);
661 }
662 }
663
664 // Treatment of boundary faces
665 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
666 const int bFISz = bdrFaceIntegrators.Size();
667 if (!factorize_face_terms && bdr_face_restrict_lex && bFISz>0)
668 {
669 // Apply the Boundary Face Restriction
670 bdr_face_restrict_lex->Mult(x, faceBdrX);
671 if (faceBdrX.Size()>0)
672 {
673 faceBdrY = 0.0;
674 // Apply the boundary face matrices
675 const int NDOFS = faceDofs;
676 auto X = Reshape(faceBdrX.Read(), NDOFS, nf_bdr);
677 auto Y = Reshape(faceBdrY.ReadWrite(), NDOFS, nf_bdr);
678 auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
679 MFEM_FORALL(glob_j, nf_bdr*NDOFS,
680 {
681 const int f = glob_j/NDOFS;
682 const int j = glob_j%NDOFS;
683 double res = 0.0;
684 for (int i = 0; i < NDOFS; i++)
685 {
686 res += A(i, j, f)*X(i, f);
687 }
688 Y(j, f) += res;
689 });
690 // Apply the Boundary Face Restriction transposed
691 bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
692 }
693 }
694 }
695
MultTranspose(const Vector & x,Vector & y) const696 void EABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
697 {
698 // Apply the Element Restriction
699 const bool useRestrict = DeviceCanUseCeed() || !elem_restrict;
700 if (!useRestrict)
701 {
702 y.UseDevice(true); // typically this is a large vector, so store on device
703 y = 0.0;
704 }
705 else
706 {
707 elem_restrict->Mult(x, localX);
708 localY = 0.0;
709 }
710 // Apply the Element Matrices transposed
711 const int NDOFS = elemDofs;
712 auto X = Reshape(useRestrict?localX.Read():x.Read(), NDOFS, ne);
713 auto Y = Reshape(useRestrict?localY.ReadWrite():y.ReadWrite(), NDOFS, ne);
714 auto A = Reshape(ea_data.Read(), NDOFS, NDOFS, ne);
715 MFEM_FORALL(glob_j, ne*NDOFS,
716 {
717 const int e = glob_j/NDOFS;
718 const int j = glob_j%NDOFS;
719 double res = 0.0;
720 for (int i = 0; i < NDOFS; i++)
721 {
722 res += A(j, i, e)*X(i, e);
723 }
724 Y(j, e) += res;
725 });
726 // Apply the Element Restriction transposed
727 if (useRestrict)
728 {
729 elem_restrict->MultTranspose(localY, y);
730 }
731
732 // Treatment of interior faces
733 Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
734 const int iFISz = intFaceIntegrators.Size();
735 if (int_face_restrict_lex && iFISz>0)
736 {
737 // Apply the Interior Face Restriction
738 int_face_restrict_lex->Mult(x, faceIntX);
739 if (faceIntX.Size()>0)
740 {
741 faceIntY = 0.0;
742 // Apply the interior face matrices transposed
743 const int NDOFS = faceDofs;
744 auto X = Reshape(faceIntX.Read(), NDOFS, 2, nf_int);
745 auto Y = Reshape(faceIntY.ReadWrite(), NDOFS, 2, nf_int);
746 if (!factorize_face_terms)
747 {
748 auto A_int = Reshape(ea_data_int.Read(), NDOFS, NDOFS, 2, nf_int);
749 MFEM_FORALL(glob_j, nf_int*NDOFS,
750 {
751 const int f = glob_j/NDOFS;
752 const int j = glob_j%NDOFS;
753 double res = 0.0;
754 for (int i = 0; i < NDOFS; i++)
755 {
756 res += A_int(j, i, 0, f)*X(i, 0, f);
757 }
758 Y(j, 0, f) += res;
759 res = 0.0;
760 for (int i = 0; i < NDOFS; i++)
761 {
762 res += A_int(j, i, 1, f)*X(i, 1, f);
763 }
764 Y(j, 1, f) += res;
765 });
766 }
767 auto A_ext = Reshape(ea_data_ext.Read(), NDOFS, NDOFS, 2, nf_int);
768 MFEM_FORALL(glob_j, nf_int*NDOFS,
769 {
770 const int f = glob_j/NDOFS;
771 const int j = glob_j%NDOFS;
772 double res = 0.0;
773 for (int i = 0; i < NDOFS; i++)
774 {
775 res += A_ext(j, i, 0, f)*X(i, 0, f);
776 }
777 Y(j, 1, f) += res;
778 res = 0.0;
779 for (int i = 0; i < NDOFS; i++)
780 {
781 res += A_ext(j, i, 1, f)*X(i, 1, f);
782 }
783 Y(j, 0, f) += res;
784 });
785 // Apply the Interior Face Restriction transposed
786 int_face_restrict_lex->AddMultTranspose(faceIntY, y);
787 }
788 }
789
790 // Treatment of boundary faces
791 Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
792 const int bFISz = bdrFaceIntegrators.Size();
793 if (!factorize_face_terms && bdr_face_restrict_lex && bFISz>0)
794 {
795 // Apply the Boundary Face Restriction
796 bdr_face_restrict_lex->Mult(x, faceBdrX);
797 if (faceBdrX.Size()>0)
798 {
799 faceBdrY = 0.0;
800 // Apply the boundary face matrices transposed
801 const int NDOFS = faceDofs;
802 auto X = Reshape(faceBdrX.Read(), NDOFS, nf_bdr);
803 auto Y = Reshape(faceBdrY.ReadWrite(), NDOFS, nf_bdr);
804 auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
805 MFEM_FORALL(glob_j, nf_bdr*NDOFS,
806 {
807 const int f = glob_j/NDOFS;
808 const int j = glob_j%NDOFS;
809 double res = 0.0;
810 for (int i = 0; i < NDOFS; i++)
811 {
812 res += A(j, i, f)*X(i, f);
813 }
814 Y(j, f) += res;
815 });
816 // Apply the Boundary Face Restriction transposed
817 bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
818 }
819 }
820 }
821
822 // Data and methods for fully-assembled bilinear forms
FABilinearFormExtension(BilinearForm * form)823 FABilinearFormExtension::FABilinearFormExtension(BilinearForm *form)
824 : EABilinearFormExtension(form),
825 mat(a->mat)
826 {
827 #ifdef MFEM_USE_MPI
828 ParFiniteElementSpace *pfes = nullptr;
829 if ( a->GetFBFI()->Size()>0 &&
830 (pfes = dynamic_cast<ParFiniteElementSpace*>(form->FESpace())) )
831 {
832 pfes->ExchangeFaceNbrData();
833 }
834 #endif
835 }
836
Assemble()837 void FABilinearFormExtension::Assemble()
838 {
839 EABilinearFormExtension::Assemble();
840 FiniteElementSpace &fes = *a->FESpace();
841 int width = fes.GetVSize();
842 int height = fes.GetVSize();
843 bool keep_nbr_block = false;
844 #ifdef MFEM_USE_MPI
845 ParFiniteElementSpace *pfes = nullptr;
846 if ( a->GetFBFI()->Size()>0 &&
847 (pfes = dynamic_cast<ParFiniteElementSpace*>(&fes)) )
848 {
849 pfes->ExchangeFaceNbrData();
850 width += pfes->GetFaceNbrVSize();
851 dg_x.SetSize(width);
852 ParBilinearForm *pb = nullptr;
853 if ((pb = dynamic_cast<ParBilinearForm*>(a)) && (pb->keep_nbr_block))
854 {
855 height += pfes->GetFaceNbrVSize();
856 dg_y.SetSize(height);
857 keep_nbr_block = true;
858 }
859 }
860 #endif
861 if (a->mat) // We reuse the sparse matrix memory
862 {
863 if (fes.IsDGSpace())
864 {
865 const L2ElementRestriction *restE =
866 static_cast<const L2ElementRestriction*>(elem_restrict);
867 const L2FaceRestriction *restF =
868 static_cast<const L2FaceRestriction*>(int_face_restrict_lex);
869 // 1. Fill J and Data
870 // 1.1 Fill J and Data with Elem ea_data
871 restE->FillJAndData(ea_data, *mat);
872 // 1.2 Fill J and Data with Face ea_data_ext
873 if (restF) { restF->FillJAndData(ea_data_ext, *mat, keep_nbr_block); }
874 // 1.3 Shift indirections in I back to original
875 auto I = mat->HostReadWriteI();
876 for (int i = height; i > 0; i--)
877 {
878 I[i] = I[i-1];
879 }
880 I[0] = 0;
881 }
882 else
883 {
884 const ElementRestriction &rest =
885 static_cast<const ElementRestriction&>(*elem_restrict);
886 rest.FillJAndData(ea_data, *mat);
887 }
888 }
889 else // We create, compute the sparsity, and fill the sparse matrix
890 {
891 mat = new SparseMatrix(height, width, 0);
892 if (fes.IsDGSpace())
893 {
894 const L2ElementRestriction *restE =
895 static_cast<const L2ElementRestriction*>(elem_restrict);
896 const L2FaceRestriction *restF =
897 static_cast<const L2FaceRestriction*>(int_face_restrict_lex);
898 // 1. Fill I
899 mat->GetMemoryI().New(height+1, mat->GetMemoryI().GetMemoryType());
900 // 1.1 Increment with restE
901 restE->FillI(*mat);
902 // 1.2 Increment with restF
903 if (restF) { restF->FillI(*mat, keep_nbr_block); }
904 // 1.3 Sum the non-zeros in I
905 auto h_I = mat->HostReadWriteI();
906 int cpt = 0;
907 for (int i = 0; i < height; i++)
908 {
909 const int nnz = h_I[i];
910 h_I[i] = cpt;
911 cpt += nnz;
912 }
913 const int nnz = cpt;
914 h_I[height] = nnz;
915 mat->GetMemoryJ().New(nnz, mat->GetMemoryJ().GetMemoryType());
916 mat->GetMemoryData().New(nnz, mat->GetMemoryData().GetMemoryType());
917 // 2. Fill J and Data
918 // 2.1 Fill J and Data with Elem ea_data
919 restE->FillJAndData(ea_data, *mat);
920 // 2.2 Fill J and Data with Face ea_data_ext
921 if (restF) { restF->FillJAndData(ea_data_ext, *mat, keep_nbr_block); }
922 // 2.3 Shift indirections in I back to original
923 auto I = mat->HostReadWriteI();
924 for (int i = height; i > 0; i--)
925 {
926 I[i] = I[i-1];
927 }
928 I[0] = 0;
929 }
930 else // continuous Galerkin case
931 {
932 const ElementRestriction &rest =
933 static_cast<const ElementRestriction&>(*elem_restrict);
934 rest.FillSparseMatrix(ea_data, *mat);
935 }
936 a->mat = mat;
937 }
938 }
939
DGMult(const Vector & x,Vector & y) const940 void FABilinearFormExtension::DGMult(const Vector &x, Vector &y) const
941 {
942 #ifdef MFEM_USE_MPI
943 const ParFiniteElementSpace *pfes;
944 if ( (pfes = dynamic_cast<const ParFiniteElementSpace*>(testFes)) )
945 {
946 // DG Prolongation
947 ParGridFunction x_gf;
948 x_gf.MakeRef(const_cast<ParFiniteElementSpace*>(pfes),
949 const_cast<Vector&>(x),0);
950 x_gf.ExchangeFaceNbrData();
951 Vector &shared_x = x_gf.FaceNbrData();
952 const int local_size = a->FESpace()->GetVSize();
953 auto dg_x_ptr = dg_x.Write();
954 auto x_ptr = x.Read();
955 MFEM_FORALL(i,local_size,
956 {
957 dg_x_ptr[i] = x_ptr[i];
958 });
959 const int shared_size = shared_x.Size();
960 auto shared_x_ptr = shared_x.Read();
961 MFEM_FORALL(i,shared_size,
962 {
963 dg_x_ptr[local_size+i] = shared_x_ptr[i];
964 });
965 ParBilinearForm *pform = nullptr;
966 if ((pform = dynamic_cast<ParBilinearForm*>(a)) && (pform->keep_nbr_block))
967 {
968 mat->Mult(dg_x, dg_y);
969 // DG Restriction
970 auto dg_y_ptr = dg_y.Read();
971 auto y_ptr = y.ReadWrite();
972 MFEM_FORALL(i,local_size,
973 {
974 y_ptr[i] += dg_y_ptr[i];
975 });
976 }
977 else
978 {
979 mat->Mult(dg_x, y);
980 }
981 }
982 else
983 #endif
984 {
985 mat->Mult(x, y);
986 }
987 }
988
Mult(const Vector & x,Vector & y) const989 void FABilinearFormExtension::Mult(const Vector &x, Vector &y) const
990 {
991 if ( a->GetFBFI()->Size()>0 )
992 {
993 DGMult(x, y);
994 }
995 else
996 {
997 mat->Mult(x, y);
998 }
999 }
1000
DGMultTranspose(const Vector & x,Vector & y) const1001 void FABilinearFormExtension::DGMultTranspose(const Vector &x, Vector &y) const
1002 {
1003 #ifdef MFEM_USE_MPI
1004 const ParFiniteElementSpace *pfes;
1005 if ( (pfes = dynamic_cast<const ParFiniteElementSpace*>(testFes)) )
1006 {
1007 // DG Prolongation
1008 ParGridFunction x_gf;
1009 x_gf.MakeRef(const_cast<ParFiniteElementSpace*>(pfes),
1010 const_cast<Vector&>(x),0);
1011 x_gf.ExchangeFaceNbrData();
1012 Vector &shared_x = x_gf.FaceNbrData();
1013 const int local_size = a->FESpace()->GetVSize();
1014 auto dg_x_ptr = dg_x.Write();
1015 auto x_ptr = x.Read();
1016 MFEM_FORALL(i,local_size,
1017 {
1018 dg_x_ptr[i] = x_ptr[i];
1019 });
1020 const int shared_size = shared_x.Size();
1021 auto shared_x_ptr = shared_x.Read();
1022 MFEM_FORALL(i,shared_size,
1023 {
1024 dg_x_ptr[local_size+i] = shared_x_ptr[i];
1025 });
1026 ParBilinearForm *pb = nullptr;
1027 if ((pb = dynamic_cast<ParBilinearForm*>(a)) && (pb->keep_nbr_block))
1028 {
1029 mat->MultTranspose(dg_x, dg_y);
1030 // DG Restriction
1031 auto dg_y_ptr = dg_y.Read();
1032 auto y_ptr = y.ReadWrite();
1033 MFEM_FORALL(i,local_size,
1034 {
1035 y_ptr[i] += dg_y_ptr[i];
1036 });
1037 }
1038 else
1039 {
1040 mat->MultTranspose(dg_x, y);
1041 }
1042 }
1043 else
1044 #endif
1045 {
1046 mat->MultTranspose(x, y);
1047 }
1048 }
1049
MultTranspose(const Vector & x,Vector & y) const1050 void FABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
1051 {
1052 if ( a->GetFBFI()->Size()>0 )
1053 {
1054 DGMultTranspose(x, y);
1055 }
1056 else
1057 {
1058 mat->MultTranspose(x, y);
1059 }
1060 }
1061
1062
MixedBilinearFormExtension(MixedBilinearForm * form)1063 MixedBilinearFormExtension::MixedBilinearFormExtension(MixedBilinearForm *form)
1064 : Operator(form->Height(), form->Width()), a(form)
1065 {
1066 // empty
1067 }
1068
GetProlongation() const1069 const Operator *MixedBilinearFormExtension::GetProlongation() const
1070 {
1071 return a->GetProlongation();
1072 }
1073
GetRestriction() const1074 const Operator *MixedBilinearFormExtension::GetRestriction() const
1075 {
1076 return a->GetRestriction();
1077 }
1078
GetOutputProlongation() const1079 const Operator *MixedBilinearFormExtension::GetOutputProlongation() const
1080 {
1081 return a->GetOutputProlongation();
1082 }
1083
GetOutputRestriction() const1084 const Operator *MixedBilinearFormExtension::GetOutputRestriction() const
1085 {
1086 return a->GetOutputRestriction();
1087 }
1088
1089 // Data and methods for partially-assembled bilinear forms
1090
PAMixedBilinearFormExtension(MixedBilinearForm * form)1091 PAMixedBilinearFormExtension::PAMixedBilinearFormExtension(
1092 MixedBilinearForm *form)
1093 : MixedBilinearFormExtension(form),
1094 trialFes(form->TrialFESpace()),
1095 testFes(form->TestFESpace()),
1096 elem_restrict_trial(NULL),
1097 elem_restrict_test(NULL)
1098 {
1099 Update();
1100 }
1101
Assemble()1102 void PAMixedBilinearFormExtension::Assemble()
1103 {
1104 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1105 const int integratorCount = integrators.Size();
1106 for (int i = 0; i < integratorCount; ++i)
1107 {
1108 integrators[i]->AssemblePA(*trialFes, *testFes);
1109 }
1110 MFEM_VERIFY(a->GetBBFI()->Size() == 0,
1111 "Partial assembly does not support AddBoundaryIntegrator yet.");
1112 MFEM_VERIFY(a->GetTFBFI()->Size() == 0,
1113 "Partial assembly does not support AddTraceFaceIntegrator yet.");
1114 MFEM_VERIFY(a->GetBTFBFI()->Size() == 0,
1115 "Partial assembly does not support AddBdrTraceFaceIntegrator yet.");
1116 }
1117
Update()1118 void PAMixedBilinearFormExtension::Update()
1119 {
1120 trialFes = a->TrialFESpace();
1121 testFes = a->TestFESpace();
1122 height = testFes->GetVSize();
1123 width = trialFes->GetVSize();
1124 elem_restrict_trial = trialFes->GetElementRestriction(
1125 ElementDofOrdering::LEXICOGRAPHIC);
1126 elem_restrict_test = testFes->GetElementRestriction(
1127 ElementDofOrdering::LEXICOGRAPHIC);
1128 if (elem_restrict_trial)
1129 {
1130 localTrial.UseDevice(true);
1131 localTrial.SetSize(elem_restrict_trial->Height(),
1132 Device::GetMemoryType());
1133 }
1134 if (elem_restrict_test)
1135 {
1136 localTest.UseDevice(true); // ensure 'localY = 0.0' is done on device
1137 localTest.SetSize(elem_restrict_test->Height(), Device::GetMemoryType());
1138 }
1139 }
1140
FormRectangularSystemOperator(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,OperatorHandle & A)1141 void PAMixedBilinearFormExtension::FormRectangularSystemOperator(
1142 const Array<int> &trial_tdof_list,
1143 const Array<int> &test_tdof_list,
1144 OperatorHandle &A)
1145 {
1146 Operator * oper;
1147 Operator::FormRectangularSystemOperator(trial_tdof_list, test_tdof_list,
1148 oper);
1149 A.Reset(oper); // A will own oper
1150 }
1151
FormRectangularLinearSystem(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B)1152 void PAMixedBilinearFormExtension::FormRectangularLinearSystem(
1153 const Array<int> &trial_tdof_list,
1154 const Array<int> &test_tdof_list,
1155 Vector &x, Vector &b,
1156 OperatorHandle &A,
1157 Vector &X, Vector &B)
1158 {
1159 Operator *oper;
1160 Operator::FormRectangularLinearSystem(trial_tdof_list, test_tdof_list, x, b,
1161 oper, X, B);
1162 A.Reset(oper); // A will own oper
1163 }
1164
SetupMultInputs(const Operator * elem_restrict_x,const Vector & x,Vector & localX,const Operator * elem_restrict_y,Vector & y,Vector & localY,const double c) const1165 void PAMixedBilinearFormExtension::SetupMultInputs(
1166 const Operator *elem_restrict_x,
1167 const Vector &x,
1168 Vector &localX,
1169 const Operator *elem_restrict_y,
1170 Vector &y,
1171 Vector &localY,
1172 const double c) const
1173 {
1174 // * G operation: localX = c*local(x)
1175 if (elem_restrict_x)
1176 {
1177 elem_restrict_x->Mult(x, localX);
1178 if (c != 1.0)
1179 {
1180 localX *= c;
1181 }
1182 }
1183 else
1184 {
1185 if (c == 1.0)
1186 {
1187 localX.SyncAliasMemory(x);
1188 }
1189 else
1190 {
1191 localX.Set(c, x);
1192 }
1193 }
1194 if (elem_restrict_y)
1195 {
1196 localY = 0.0;
1197 }
1198 else
1199 {
1200 y.UseDevice(true);
1201 localY.SyncAliasMemory(y);
1202 }
1203 }
1204
Mult(const Vector & x,Vector & y) const1205 void PAMixedBilinearFormExtension::Mult(const Vector &x, Vector &y) const
1206 {
1207 y = 0.0;
1208 AddMult(x, y);
1209 }
1210
AddMult(const Vector & x,Vector & y,const double c) const1211 void PAMixedBilinearFormExtension::AddMult(const Vector &x, Vector &y,
1212 const double c) const
1213 {
1214 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1215 const int iSz = integrators.Size();
1216
1217 // * G operation
1218 SetupMultInputs(elem_restrict_trial, x, localTrial,
1219 elem_restrict_test, y, localTest, c);
1220
1221 // * B^TDB operation
1222 for (int i = 0; i < iSz; ++i)
1223 {
1224 integrators[i]->AddMultPA(localTrial, localTest);
1225 }
1226
1227 // * G^T operation
1228 if (elem_restrict_test)
1229 {
1230 tempY.SetSize(y.Size());
1231 elem_restrict_test->MultTranspose(localTest, tempY);
1232 y += tempY;
1233 }
1234 }
1235
MultTranspose(const Vector & x,Vector & y) const1236 void PAMixedBilinearFormExtension::MultTranspose(const Vector &x,
1237 Vector &y) const
1238 {
1239 y = 0.0;
1240 AddMultTranspose(x, y);
1241 }
1242
AddMultTranspose(const Vector & x,Vector & y,const double c) const1243 void PAMixedBilinearFormExtension::AddMultTranspose(const Vector &x, Vector &y,
1244 const double c) const
1245 {
1246 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1247 const int iSz = integrators.Size();
1248
1249 // * G operation
1250 SetupMultInputs(elem_restrict_test, x, localTest,
1251 elem_restrict_trial, y, localTrial, c);
1252
1253 // * B^TD^TB operation
1254 for (int i = 0; i < iSz; ++i)
1255 {
1256 integrators[i]->AddMultTransposePA(localTest, localTrial);
1257 }
1258
1259 // * G^T operation
1260 if (elem_restrict_trial)
1261 {
1262 tempY.SetSize(y.Size());
1263 elem_restrict_trial->MultTranspose(localTrial, tempY);
1264 y += tempY;
1265 }
1266 }
1267
AssembleDiagonal_ADAt(const Vector & D,Vector & diag) const1268 void PAMixedBilinearFormExtension::AssembleDiagonal_ADAt(const Vector &D,
1269 Vector &diag) const
1270 {
1271 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1272
1273 const int iSz = integrators.Size();
1274
1275 if (elem_restrict_trial)
1276 {
1277 const ElementRestriction* H1elem_restrict_trial =
1278 dynamic_cast<const ElementRestriction*>(elem_restrict_trial);
1279 if (H1elem_restrict_trial)
1280 {
1281 H1elem_restrict_trial->MultUnsigned(D, localTrial);
1282 }
1283 else
1284 {
1285 elem_restrict_trial->Mult(D, localTrial);
1286 }
1287 }
1288
1289 if (elem_restrict_test)
1290 {
1291 localTest = 0.0;
1292 for (int i = 0; i < iSz; ++i)
1293 {
1294 if (elem_restrict_trial)
1295 {
1296 integrators[i]->AssembleDiagonalPA_ADAt(localTrial, localTest);
1297 }
1298 else
1299 {
1300 integrators[i]->AssembleDiagonalPA_ADAt(D, localTest);
1301 }
1302 }
1303 const ElementRestriction* H1elem_restrict_test =
1304 dynamic_cast<const ElementRestriction*>(elem_restrict_test);
1305 if (H1elem_restrict_test)
1306 {
1307 H1elem_restrict_test->MultTransposeUnsigned(localTest, diag);
1308 }
1309 else
1310 {
1311 elem_restrict_test->MultTranspose(localTest, diag);
1312 }
1313 }
1314 else
1315 {
1316 diag.UseDevice(true); // typically this is a large vector, so store on device
1317 diag = 0.0;
1318 for (int i = 0; i < iSz; ++i)
1319 {
1320 if (elem_restrict_trial)
1321 {
1322 integrators[i]->AssembleDiagonalPA_ADAt(localTrial, diag);
1323 }
1324 else
1325 {
1326 integrators[i]->AssembleDiagonalPA_ADAt(D, diag);
1327 }
1328 }
1329 }
1330 }
1331
PADiscreteLinearOperatorExtension(DiscreteLinearOperator * linop)1332 PADiscreteLinearOperatorExtension::PADiscreteLinearOperatorExtension(
1333 DiscreteLinearOperator *linop) :
1334 PAMixedBilinearFormExtension(linop)
1335 {
1336 }
1337
1338 const
GetOutputRestrictionTranspose() const1339 Operator *PADiscreteLinearOperatorExtension::GetOutputRestrictionTranspose()
1340 const
1341 {
1342 return a->GetOutputRestrictionTranspose();
1343 }
1344
Assemble()1345 void PADiscreteLinearOperatorExtension::Assemble()
1346 {
1347 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1348 const int integratorCount = integrators.Size();
1349 for (int i = 0; i < integratorCount; ++i)
1350 {
1351 integrators[i]->AssemblePA(*trialFes, *testFes);
1352 }
1353
1354 test_multiplicity.UseDevice(true);
1355 test_multiplicity.SetSize(elem_restrict_test->Width()); // l-vector
1356 Vector ones(elem_restrict_test->Height()); // e-vector
1357 ones = 1.0;
1358
1359 const ElementRestriction* elem_restrict =
1360 dynamic_cast<const ElementRestriction*>(elem_restrict_test);
1361 if (elem_restrict)
1362 {
1363 elem_restrict->MultTransposeUnsigned(ones, test_multiplicity);
1364 }
1365 else
1366 {
1367 mfem_error("A real ElementRestriction is required in this setting!");
1368 }
1369
1370 auto tm = test_multiplicity.ReadWrite();
1371 MFEM_FORALL(i, test_multiplicity.Size(),
1372 {
1373 tm[i] = 1.0 / tm[i];
1374 });
1375 }
1376
AddMult(const Vector & x,Vector & y,const double c) const1377 void PADiscreteLinearOperatorExtension::AddMult(
1378 const Vector &x, Vector &y, const double c) const
1379 {
1380 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1381 const int iSz = integrators.Size();
1382
1383 // * G operation
1384 SetupMultInputs(elem_restrict_trial, x, localTrial,
1385 elem_restrict_test, y, localTest, c);
1386
1387 // * B^TDB operation
1388 for (int i = 0; i < iSz; ++i)
1389 {
1390 integrators[i]->AddMultPA(localTrial, localTest);
1391 }
1392
1393 // do a kind of "set" rather than "add" in the below
1394 // operation as compared to the BilinearForm case
1395 // * G^T operation (kind of...)
1396 const ElementRestriction* elem_restrict =
1397 dynamic_cast<const ElementRestriction*>(elem_restrict_test);
1398 if (elem_restrict)
1399 {
1400 tempY.SetSize(y.Size());
1401 elem_restrict->MultLeftInverse(localTest, tempY);
1402 y += tempY;
1403 }
1404 else
1405 {
1406 mfem_error("In this setting you need a real ElementRestriction!");
1407 }
1408 }
1409
AddMultTranspose(const Vector & x,Vector & y,const double c) const1410 void PADiscreteLinearOperatorExtension::AddMultTranspose(
1411 const Vector &x, Vector &y, const double c) const
1412 {
1413 Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
1414 const int iSz = integrators.Size();
1415
1416 // do a kind of "set" rather than "add" in the below
1417 // operation as compared to the BilinearForm case
1418 // * G operation (kinda)
1419 Vector xscaled(x);
1420 MFEM_VERIFY(x.Size() == test_multiplicity.Size(), "Input vector of wrong size");
1421 auto xs = xscaled.ReadWrite();
1422 auto tm = test_multiplicity.Read();
1423 MFEM_FORALL(i, x.Size(),
1424 {
1425 xs[i] *= tm[i];
1426 });
1427 SetupMultInputs(elem_restrict_test, xscaled, localTest,
1428 elem_restrict_trial, y, localTrial, c);
1429
1430 // * B^TD^TB operation
1431 for (int i = 0; i < iSz; ++i)
1432 {
1433 integrators[i]->AddMultTransposePA(localTest, localTrial);
1434 }
1435
1436 // * G^T operation
1437 if (elem_restrict_trial)
1438 {
1439 tempY.SetSize(y.Size());
1440 elem_restrict_trial->MultTranspose(localTrial, tempY);
1441 y += tempY;
1442 }
1443 else
1444 {
1445 mfem_error("Trial ElementRestriction not defined");
1446 }
1447 }
1448
FormRectangularSystemOperator(const Array<int> & ess1,const Array<int> & ess2,OperatorHandle & A)1449 void PADiscreteLinearOperatorExtension::FormRectangularSystemOperator(
1450 const Array<int>& ess1, const Array<int>& ess2, OperatorHandle &A)
1451 {
1452 const Operator *Pi = this->GetProlongation();
1453 const Operator *RoT = this->GetOutputRestrictionTranspose();
1454 Operator *rap = SetupRAP(Pi, RoT);
1455
1456 RectangularConstrainedOperator *Arco
1457 = new RectangularConstrainedOperator(rap, ess1, ess2, rap != this);
1458
1459 A.Reset(Arco);
1460 }
1461
1462 } // namespace mfem
1463