MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
blockstaticcond.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "blockstaticcond.hpp"
13
14namespace mfem
15{
16
18 fes_)
19{
20 SetSpaces(fes_);
21
22 Array<int> rvdofs;
23 Array<int> vdofs;
24 Array<int> rdof_edof0;
25 for (int k = 0; k<nblocks; k++)
26 {
27 if (!tr_fes[k]) { continue; }
28 rdof_edof0.SetSize(tr_fes[k]->GetVSize());
29 for (int i = 0; i < mesh->GetNE(); i++)
30 {
31 fes[k]->GetElementVDofs(i, vdofs);
32 tr_fes[k]->GetElementVDofs(i, rvdofs);
33 const int vdim = fes[k]->GetVDim();
34 const int nsd = vdofs.Size()/vdim;
35 const int nsrd = rvdofs.Size()/vdim;
36 for (int vd = 0; vd < vdim; vd++)
37 {
38 for (int j = 0; j < nsrd; j++)
39 {
40 int rvdof = rvdofs[j+nsrd*vd];
41 int vdof = vdofs[j+nsd*vd];
42 if (rvdof < 0)
43 {
44 rvdof = -1-rvdof;
45 vdof = -1-vdof;
46 }
47 MFEM_ASSERT(vdof >= 0, "incompatible volume and trace FE spaces");
48 rdof_edof0[rvdof] = vdof + dof_offsets[k];
49 }
50 }
51 }
52 rdof_edof.Append(rdof_edof0);
53 }
54}
55
56void BlockStaticCondensation::SetSpaces(Array<FiniteElementSpace*> & fes_)
57{
58#ifdef MFEM_USE_MPI
59 ParMesh *pmesh = nullptr;
60 parallel = false;
61 if (dynamic_cast<ParFiniteElementSpace *>(fes_[0]))
62 {
63 parallel = true;
64 }
65#else
66 parallel = false;
67#endif
68 fes=fes_;
69 nblocks = fes.Size();
70 rblocks = 0;
71 tr_fes.SetSize(nblocks);
72 mesh = fes[0]->GetMesh();
73
74 IsTraceSpace.SetSize(nblocks);
75 const FiniteElementCollection * fec;
76 for (int i = 0; i < nblocks; i++)
77 {
78 fec = fes[i]->FEColl();
79 IsTraceSpace[i] =
80 (dynamic_cast<const H1_Trace_FECollection*>(fec) ||
81 dynamic_cast<const ND_Trace_FECollection*>(fec) ||
82 dynamic_cast<const RT_Trace_FECollection*>(fec));
83#ifdef MFEM_USE_MPI
84 if (parallel)
85 {
86 pmesh = dynamic_cast<ParMesh *>(mesh);
87 tr_fes[i] = (fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) ?
88 nullptr : (IsTraceSpace[i]) ? fes[i] :
89 new ParFiniteElementSpace(pmesh, fec->GetTraceCollection(), fes[i]->GetVDim(),
90 fes[i]->GetOrdering());
91 }
92 else
93 {
94 tr_fes[i] = (fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) ?
95 nullptr : (IsTraceSpace[i]) ? fes[i] :
96 new FiniteElementSpace(mesh, fec->GetTraceCollection(), fes[i]->GetVDim(),
97 fes[i]->GetOrdering());
98 }
99#else
100 // skip if it's an L2 space (no trace space to construct)
101 tr_fes[i] = (fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) ?
102 nullptr : (IsTraceSpace[i]) ? fes[i] :
103 new FiniteElementSpace(mesh, fec->GetTraceCollection(), fes[i]->GetVDim(),
104 fes[i]->GetOrdering());
105#endif
106 if (tr_fes[i]) { rblocks++; }
107 }
108 if (parallel)
109 {
110 ess_tdofs.SetSize(rblocks);
111 for (int i = 0; i<rblocks; i++)
112 {
113 ess_tdofs[i] = new Array<int>();
114 }
115 }
116 Init();
117}
118
119void BlockStaticCondensation::ComputeOffsets()
120{
121 dof_offsets.SetSize(nblocks+1);
122 tdof_offsets.SetSize(nblocks+1);
123 dof_offsets[0] = 0;
124 tdof_offsets[0] = 0;
125
126 rdof_offsets.SetSize(rblocks+1);
127 rtdof_offsets.SetSize(rblocks+1);
128 rdof_offsets[0] = 0;
129 rtdof_offsets[0] = 0;
130
131 int j=0;
132 for (int i =0; i<nblocks; i++)
133 {
134 dof_offsets[i+1] = fes[i]->GetVSize();
135 tdof_offsets[i+1] = fes[i]->GetTrueVSize();
136 if (tr_fes[i])
137 {
138 rdof_offsets[j+1] = tr_fes[i]->GetVSize();
139 rtdof_offsets[j+1] = tr_fes[i]->GetTrueVSize();
140 j++;
141 }
142 }
143 rdof_offsets.PartialSum();
144 rtdof_offsets.PartialSum();
145 dof_offsets.PartialSum();
146 tdof_offsets.PartialSum();
147}
148
149
150void BlockStaticCondensation::Init()
151{
152 lmat.SetSize(mesh->GetNE());
153 lvec.SetSize(mesh->GetNE());
154 for (int i = 0; i < mesh->GetNE(); i++)
155 {
156 lmat[i] = nullptr;
157 lvec[i] = nullptr;
158 }
159
160 ComputeOffsets();
161
162 S = new BlockMatrix(rdof_offsets);
163 S->owns_blocks = 1;
164
165 for (int i = 0; i<S->NumRowBlocks(); i++)
166 {
167 int h = rdof_offsets[i+1] - rdof_offsets[i];
168 for (int j = 0; j<S->NumColBlocks(); j++)
169 {
170 int w = rdof_offsets[j+1] - rdof_offsets[j];
171 S->SetBlock(i,j,new SparseMatrix(h, w));
172 }
173 }
174 y = new BlockVector(rdof_offsets);
175 *y = 0.;
176}
177
178void BlockStaticCondensation::GetReducedElementIndicesAndOffsets(int el,
179 Array<int> & trace_ldofs,
180 Array<int> & interior_ldofs,
181 Array<int> & offsets) const
182{
183 int dim = mesh->Dimension();
184 offsets.SetSize(tr_fes.Size()+1); offsets = 0;
185 Array<int> dofs;
186 Array<int> faces, ori;
187 if (dim == 1)
188 {
189 mesh->GetElementVertices(el, faces);
190 }
191 else if (dim == 2)
192 {
193 mesh->GetElementEdges(el, faces, ori);
194 }
195 else if (dim == 3)
196 {
197 mesh->GetElementFaces(el,faces,ori);
198 }
199 else
200 {
201 MFEM_ABORT("BlockStaticCondensation::GetReducedElementIndicesAndOffsets: "
202 "dim > 3 not supported");
203 }
204 int numfaces = faces.Size();
205
206 trace_ldofs.SetSize(0);
207 interior_ldofs.SetSize(0);
208 // construct Array of bubble dofs to be extracted
209 int skip=0;
210 Array<int> tr_dofs;
211 Array<int> int_dofs;
212 for (int i = 0; i<tr_fes.Size(); i++)
213 {
214 int td = 0;
215 int ndof;
216 // if it's an L2 space (bubbles)
217 if (!tr_fes[i])
218 {
219 ndof = fes[i]->GetVDim()*fes[i]->GetFE(el)->GetDof();
220 td = 0;
221 }
222 else if (IsTraceSpace[i])
223 {
224 for (int iface = 0; iface < numfaces; iface++)
225 {
226 td += fes[i]->GetVDim()*fes[i]->GetFaceElement(faces[iface])->GetDof();
227 }
228 ndof = td;
229 }
230 else
231 {
232 Array<int> trace_dofs;
233 ndof = fes[i]->GetVDim()*fes[i]->GetFE(el)->GetDof();
234 tr_fes[i]->GetElementVDofs(el, trace_dofs);
235 td = trace_dofs.Size(); // number of trace dofs
236 }
237 offsets[i+1] = td;
238 tr_dofs.SetSize(td);
239 int_dofs.SetSize(ndof - td);
240 for (int j = 0; j<td; j++)
241 {
242 tr_dofs[j] = skip + j;
243 }
244 for (int j = 0; j<ndof-td; j++)
245 {
246 int_dofs[j] = skip + td + j;
247 }
248 skip+=ndof;
249
250 trace_ldofs.Append(tr_dofs);
251 interior_ldofs.Append(int_dofs);
252 }
253 offsets.PartialSum();
254}
255
256
257void BlockStaticCondensation::GetReducedElementVDofs(int el,
258 Array<int> & rdofs) const
259{
260 Array<int> faces, ori;
261 int dim = mesh->Dimension();
262 if (dim == 1)
263 {
264 mesh->GetElementVertices(el, faces);
265 }
266 else if (dim == 2)
267 {
268 mesh->GetElementEdges(el, faces, ori);
269 }
270 else if (dim == 3)
271 {
272 mesh->GetElementFaces(el,faces,ori);
273 }
274 else
275 {
276 MFEM_ABORT("BlockStaticCondensation::GetReducedElementVDofs: "
277 "dim > 3 not supported");
278 }
279 int numfaces = faces.Size();
280 rdofs.SetSize(0);
281 int skip = 0;
282 for (int i = 0; i<tr_fes.Size(); i++)
283 {
284 if (!tr_fes[i]) { continue; }
285 Array<int> vdofs;
286 if (IsTraceSpace[i])
287 {
288 Array<int> face_vdofs;
289 for (int k = 0; k < numfaces; k++)
290 {
291 int iface = faces[k];
292 tr_fes[i]->GetFaceVDofs(iface, face_vdofs);
293 vdofs.Append(face_vdofs);
294 }
295 }
296 else
297 {
298 tr_fes[i]->GetElementVDofs(el, vdofs);
299 }
300 for (int j=0; j<vdofs.Size(); j++)
301 {
302 vdofs[j] = (vdofs[j]>=0) ? vdofs[j]+rdof_offsets[skip] :
303 vdofs[j]-rdof_offsets[skip];
304 }
305 skip++;
306 rdofs.Append(vdofs);
307 }
308}
309
310void BlockStaticCondensation::GetElementVDofs(int el, Array<int> & vdofs) const
311{
312 Array<int> faces, ori;
313 int dim = mesh->Dimension();
314 if (dim == 1)
315 {
316 mesh->GetElementVertices(el, faces);
317 }
318 else if (dim == 2)
319 {
320 mesh->GetElementEdges(el, faces, ori);
321 }
322 else if (dim == 3)
323 {
324 mesh->GetElementFaces(el,faces,ori);
325 }
326 else
327 {
328 MFEM_ABORT("BlockStaticCondensation::GetElementVDofs: "
329 "dim > 3 not supported");
330 }
331 int numfaces = faces.Size();
332 vdofs.SetSize(0);
333 for (int i = 0; i<tr_fes.Size(); i++)
334 {
335 Array<int> dofs;
336 if (IsTraceSpace[i])
337 {
338 Array<int> face_vdofs;
339 for (int k = 0; k < numfaces; k++)
340 {
341 int iface = faces[k];
342 fes[i]->GetFaceVDofs(iface, face_vdofs);
343 dofs.Append(face_vdofs);
344 }
345 }
346 else
347 {
348 fes[i]->GetElementVDofs(el, dofs);
349 }
350 for (int j=0; j<dofs.Size(); j++)
351 {
352 dofs[j] = (dofs[j]>=0) ? dofs[j]+dof_offsets[i] :
353 dofs[j]-dof_offsets[i];
354 }
355 vdofs.Append(dofs);
356 }
357}
358
359
360void BlockStaticCondensation::GetLocalSchurComplement(int el,
361 const Array<int> & tr_idx,
362 const Array<int> & int_idx,
363 const DenseMatrix & elmat,
364 const Vector & elvect,
365 DenseMatrix & rmat,
366 Vector & rvect)
367{
368 int rdofs = tr_idx.Size();
369 int idofs = int_idx.Size();
370 MFEM_VERIFY(idofs != 0, "Number of interior dofs is zero");
371 MFEM_VERIFY(rdofs != 0, "Number of interface dofs is zero");
372
373 rmat.SetSize(rdofs);
374 rvect.SetSize(rdofs);
375
376 DenseMatrix A_tt, A_ti, A_it, A_ii;
377 Vector y_t, y_i;
378
379 elmat.GetSubMatrix(tr_idx,A_tt);
380 elmat.GetSubMatrix(tr_idx,int_idx, A_ti);
381 elmat.GetSubMatrix(int_idx, tr_idx, A_it);
382 elmat.GetSubMatrix(int_idx, A_ii);
383
384 elvect.GetSubVector(tr_idx, y_t);
385 elvect.GetSubVector(int_idx, y_i);
386
387 DenseMatrixInverse lu(A_ii);
388 lu.Factor();
389 lmat[el] = new DenseMatrix(idofs,rdofs);
390 lvec[el] = new Vector(idofs);
391
392 lu.Mult(A_it,*lmat[el]);
393 lu.Mult(y_i,*lvec[el]);
394
395 // LHS
396 mfem::Mult(A_ti,*lmat[el],rmat);
397
398 rmat.Neg();
399 rmat.Add(1., A_tt);
400
401 // RHS
402 A_ti.Mult(*lvec[el], rvect);
403 rvect.Neg();
404 rvect.Add(1., y_t);
405}
406
407
409 DenseMatrix &elmat,
410 Vector & elvect)
411{
412 // Get Schur Complement
413 Array<int> tr_idx, int_idx;
414 Array<int> offsets;
415 // Get local element idx and offsets for global assembly
416 GetReducedElementIndicesAndOffsets(el, tr_idx,int_idx, offsets);
417
418 DenseMatrix rmat, *rmatptr;
419 Vector rvec, *rvecptr;
420 // Extract the reduced matrices based on tr_idx and int_idx
421 if (int_idx.Size()!=0)
422 {
423 GetLocalSchurComplement(el,tr_idx,int_idx, elmat, elvect, rmat, rvec);
424 rmatptr = &rmat;
425 rvecptr = &rvec;
426 }
427 else
428 {
429 rmatptr = &elmat;
430 rvecptr = &elvect;
431 }
432
433 // Assemble global mat and rhs
434 DofTransformation * doftrans_i, *doftrans_j;
435
436 Array<int> faces, ori;
437 int dim = mesh->Dimension();
438 if (dim == 1)
439 {
440 mesh->GetElementVertices(el, faces);
441 }
442 else if (dim == 2)
443 {
444 mesh->GetElementEdges(el, faces, ori);
445 }
446 else if (dim == 3)
447 {
448 mesh->GetElementFaces(el,faces,ori);
449 }
450 else
451 {
452 MFEM_ABORT("BlockStaticCondensation::AssembleReducedSystem: "
453 "dim > 3 not supported");
454 }
455 int numfaces = faces.Size();
456
457 int skip_i=0;
458 for (int i = 0; i<tr_fes.Size(); i++)
459 {
460 if (!tr_fes[i]) { continue; }
461 Array<int> vdofs_i;
462 doftrans_i = nullptr;
463 if (IsTraceSpace[i])
464 {
465 Array<int> face_vdofs;
466 for (int k = 0; k < numfaces; k++)
467 {
468 int iface = faces[k];
469 tr_fes[i]->GetFaceVDofs(iface, face_vdofs);
470 vdofs_i.Append(face_vdofs);
471 }
472 }
473 else
474 {
475 doftrans_i = tr_fes[i]->GetElementVDofs(el, vdofs_i);
476 }
477 int skip_j=0;
478 for (int j = 0; j<tr_fes.Size(); j++)
479 {
480 if (!tr_fes[j]) { continue; }
481 Array<int> vdofs_j;
482 doftrans_j = nullptr;
483
484 if (IsTraceSpace[j])
485 {
486 Array<int> face_vdofs;
487 for (int k = 0; k < numfaces; k++)
488 {
489 int iface = faces[k];
490 tr_fes[j]->GetFaceVDofs(iface, face_vdofs);
491 vdofs_j.Append(face_vdofs);
492 }
493 }
494 else
495 {
496 doftrans_j = tr_fes[j]->GetElementVDofs(el, vdofs_j);
497 }
498
499 DenseMatrix Ae;
500 rmatptr->GetSubMatrix(offsets[i],offsets[i+1],
501 offsets[j],offsets[j+1], Ae);
502 if (doftrans_i || doftrans_j)
503 {
504 TransformDual(doftrans_i, doftrans_j, Ae);
505 }
506 S->GetBlock(skip_i,skip_j).AddSubMatrix(vdofs_i,vdofs_j, Ae);
507 skip_j++;
508 }
509
510 // assemble rhs
511 real_t * data = rvecptr->GetData();
512 Vector vec1;
513 // ref subvector
514 vec1.SetDataAndSize(&data[offsets[i]],
515 offsets[i+1]-offsets[i]);
516 if (doftrans_i)
517 {
518 doftrans_i->TransformDual(vec1);
519 }
520 y->GetBlock(skip_i).AddElementVector(vdofs_i,vec1);
521 skip_i++;
522 }
523}
524
525void BlockStaticCondensation::BuildProlongation()
526{
527 P = new BlockMatrix(rdof_offsets, rtdof_offsets);
528 R = new BlockMatrix(rtdof_offsets, rdof_offsets);
529 P->owns_blocks = 0;
530 R->owns_blocks = 0;
531 int skip = 0;
532 for (int i = 0; i<nblocks; i++)
533 {
534 if (!tr_fes[i]) { continue; }
535 const SparseMatrix *P_ = tr_fes[i]->GetConformingProlongation();
536 if (P_)
537 {
538 const SparseMatrix *R_ = tr_fes[i]->GetRestrictionMatrix();
539 P->SetBlock(skip,skip,const_cast<SparseMatrix*>(P_));
540 R->SetBlock(skip,skip,const_cast<SparseMatrix*>(R_));
541 }
542 skip++;
543 }
544}
545
546#ifdef MFEM_USE_MPI
547void BlockStaticCondensation::BuildParallelProlongation()
548{
549 MFEM_VERIFY(parallel, "BuildParallelProlongation: wrong code path");
550 pP = new BlockOperator(rdof_offsets, rtdof_offsets);
551 R = new BlockMatrix(rtdof_offsets, rdof_offsets);
552 pP->owns_blocks = 0;
553 R->owns_blocks = 0;
554 int skip = 0;
555 for (int i = 0; i<nblocks; i++)
556 {
557 if (!tr_fes[i]) { continue; }
558 const HypreParMatrix *P_ =
559 dynamic_cast<ParFiniteElementSpace *>(tr_fes[i])->Dof_TrueDof_Matrix();
560 if (P_)
561 {
562 const SparseMatrix *R_ = tr_fes[i]->GetRestrictionMatrix();
563 pP->SetBlock(skip,skip,const_cast<HypreParMatrix*>(P_));
564 R->SetBlock(skip,skip,const_cast<SparseMatrix*>(R_));
565 }
566 skip++;
567 }
568}
569
571{
572 if (!pP) { BuildParallelProlongation(); }
573
574 pS = new BlockOperator(rtdof_offsets);
575 pS_e = new BlockOperator(rtdof_offsets);
576 pS->owns_blocks = 1;
577 pS_e->owns_blocks = 1;
578 HypreParMatrix * A = nullptr;
579 HypreParMatrix * PtAP = nullptr;
580 int skip_i=0;
581 ParFiniteElementSpace * pfes_i = nullptr;
582 ParFiniteElementSpace * pfes_j = nullptr;
583 for (int i = 0; i<nblocks; i++)
584 {
585 if (!tr_fes[i]) { continue; }
586 pfes_i = dynamic_cast<ParFiniteElementSpace*>(tr_fes[i]);
587 HypreParMatrix * Pi = (HypreParMatrix*)(&pP->GetBlock(skip_i,skip_i));
588 int skip_j=0;
589 for (int j = 0; j<nblocks; j++)
590 {
591 if (!tr_fes[j]) { continue; }
592 if (m->IsZeroBlock(skip_i,skip_j)) { continue; }
593 if (skip_i == skip_j)
594 {
595 // Make block diagonal square hypre matrix
596 A = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
597 pfes_i->GetDofOffsets(),&m->GetBlock(skip_i,skip_i));
598 PtAP = RAP(A,Pi);
599 delete A;
600 pS_e->SetBlock(skip_i,skip_i,PtAP->EliminateRowsCols(*ess_tdofs[skip_i]));
601 }
602 else
603 {
604 pfes_j = dynamic_cast<ParFiniteElementSpace*>(tr_fes[j]);
605 HypreParMatrix * Pj = (HypreParMatrix*)(&pP->GetBlock(skip_j,skip_j));
606 A = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
607 pfes_j->GlobalVSize(), pfes_i->GetDofOffsets(),
608 pfes_j->GetDofOffsets(), &m->GetBlock(skip_i,skip_j));
609 PtAP = RAP(Pi,A,Pj);
610 delete A;
611 pS_e->SetBlock(skip_i,skip_j,PtAP->EliminateCols(*ess_tdofs[skip_j]));
612 PtAP->EliminateRows(*ess_tdofs[skip_i]);
613 }
614 pS->SetBlock(skip_i,skip_j,PtAP);
615 skip_j++;
616 }
617 skip_i++;
618 }
619}
620
621#endif
622
623
624void BlockStaticCondensation::ConformingAssemble(int skip_zeros)
625{
626 Finalize(0);
627 if (!P) { BuildProlongation(); }
628
629 BlockMatrix * Pt = Transpose(*P);
630 BlockMatrix * PtA = mfem::Mult(*Pt, *S);
631 delete S;
632 if (S_e)
633 {
634 BlockMatrix *PtAe = mfem::Mult(*Pt, *S_e);
635 delete S_e;
636 S_e = PtAe;
637 }
638 delete Pt;
639 S = mfem::Mult(*PtA, *P);
640 delete PtA;
641
642 if (S_e)
643 {
644 BlockMatrix *PtAeP = mfem::Mult(*S_e, *P);
645 S_e = PtAeP;
646 }
647 height = S->Height();
648 width = S->Width();
649}
650
652{
653 if (S) { S->Finalize(skip_zeros); }
654 if (S_e) { S_e->Finalize(skip_zeros); }
655}
656
658 diag_policy)
659{
660 if (!parallel)
661 {
662 if (!S_e)
663 {
664 bool conforming = true;
665 for (int i = 0; i<nblocks; i++)
666 {
667 if (!tr_fes[i]) { continue; }
668 const SparseMatrix *P_ = tr_fes[i]->GetConformingProlongation();
669 if (P_)
670 {
671 conforming = false;
672 break;
673 }
674 }
675 if (!conforming) { ConformingAssemble(0); }
676 const int remove_zeros = 0;
677 EliminateReducedTrueDofs(ess_rtdof_list, diag_policy);
678 Finalize(remove_zeros);
679 }
680 }
681 else
682 {
683#ifdef MFEM_USE_MPI
684 FillEssTdofLists(ess_rtdof_list);
685 if (S)
686 {
687 const int remove_zeros = 0;
688 Finalize(remove_zeros);
690 delete S; S=nullptr;
691 delete S_e; S_e = nullptr;
692 }
693#endif
694 }
695}
696
697
698void BlockStaticCondensation::ConvertMarkerToReducedTrueDofs(
699 Array<int> & tdof_marker,
700 Array<int> & rtdof_marker)
701{
702 // convert tdof_marker to dof_marker
703 rtdof_marker.SetSize(0);
704 Array<int> tdof_marker0;
705 Array<int> dof_marker0;
706 Array<int> dof_marker;
707 int * data = tdof_marker.GetData();
708 for (int i = 0; i<nblocks; i++)
709 {
710 tdof_marker0.MakeRef(&data[tdof_offsets[i]],tdof_offsets[i+1]-tdof_offsets[i]);
711 const SparseMatrix * R_ = fes[i]->GetRestrictionMatrix();
712 if (!R_)
713 {
714 dof_marker0.MakeRef(tdof_marker0);
715 }
716 else
717 {
718 dof_marker0.SetSize(fes[i]->GetVSize());
719 R_->BooleanMultTranspose(tdof_marker0, dof_marker0);
720 }
721 dof_marker.Append(dof_marker0);
722 }
723
724 int rdofs = rdof_edof.Size();
725 Array<int> rdof_marker(rdofs);
726
727 for (int i = 0; i < rdofs; i++)
728 {
729 rdof_marker[i] = dof_marker[rdof_edof[i]];
730 }
731
732 // convert rdof_marker to rtdof_marker
733 Array<int> rtdof_marker0;
734 Array<int> rdof_marker0;
735 int * rdata = rdof_marker.GetData();
736 int k=0;
737 for (int i = 0; i<nblocks; i++)
738 {
739 if (!tr_fes[i]) { continue; }
740 rdof_marker0.MakeRef(&rdata[rdof_offsets[k]],rdof_offsets[k+1]-rdof_offsets[k]);
741 const SparseMatrix *tr_R = tr_fes[i]->GetRestrictionMatrix();
742 if (!tr_R)
743 {
744 rtdof_marker0.MakeRef(rdof_marker0);
745 }
746 else
747 {
748 rtdof_marker0.SetSize(tr_fes[i]->GetTrueVSize());
749 tr_R->BooleanMult(rdof_marker0, rtdof_marker0);
750 }
751 rtdof_marker.Append(rtdof_marker0);
752 k++;
753 }
754}
755
756void BlockStaticCondensation::FillEssTdofLists(const Array<int> & ess_tdof_list)
757{
758 int j;
759 for (int i = 0; i<ess_tdof_list.Size(); i++)
760 {
761 int tdof = ess_tdof_list[i];
762 for (j = 0; j < rblocks; j++)
763 {
764 if (rtdof_offsets[j+1] > tdof) { break; }
765 }
766 ess_tdofs[j]->Append(tdof-rtdof_offsets[j]);
767 }
768}
769
771 &ess_tdof_list)
772{
773 Array<int> tdof_marker;
774 Array<int> rtdof_marker;
775 FiniteElementSpace::ListToMarker(ess_tdof_list,tdof_offsets.Last(),tdof_marker);
776 ConvertMarkerToReducedTrueDofs(tdof_marker, rtdof_marker);
777 FiniteElementSpace::MarkerToList(rtdof_marker,ess_rtdof_list);
778}
779
781 &ess_rtdof_list_,
783{
784 MFEM_VERIFY(!parallel, "EliminateReducedTrueDofs::Wrong Code path");
785
786 if (S_e == NULL)
787 {
788 Array<int> offsets;
789
790 offsets.MakeRef( (P) ? rtdof_offsets : rdof_offsets);
791
792 S_e = new BlockMatrix(offsets);
793 S_e->owns_blocks = 1;
794 for (int i = 0; i<S_e->NumRowBlocks(); i++)
795 {
796 int h = offsets[i+1] - offsets[i];
797 for (int j = 0; j<S_e->NumColBlocks(); j++)
798 {
799 int w = offsets[j+1] - offsets[j];
800 S_e->SetBlock(i,j,new SparseMatrix(h, w));
801 }
802 }
803 }
804 S->EliminateRowCols(ess_rtdof_list_,S_e,dpolicy);
805}
806
808 Vector &sc_sol) const
809{
810 MFEM_ASSERT(sol.Size() == dof_offsets.Last(), "'sol' has incorrect size");
811 const int nrdofs = rdof_offsets.Last();
812 Vector sol_r;
813 if (!R)
814 {
815 sc_sol.SetSize(nrdofs);
816 sol_r.SetDataAndSize(sc_sol.GetData(), sc_sol.Size());
817 }
818 else
819 {
820 sol_r.SetSize(nrdofs);
821 }
822 for (int i = 0; i < nrdofs; i++)
823 {
824 sol_r(i) = sol(rdof_edof[i]);
825 }
826 if (R)
827 {
828 // wrap vector into a block vector
829 BlockVector blsol_r(sol_r,rdof_offsets);
830 sc_sol.SetSize(R->Height());
831 R->Mult(blsol_r, sc_sol);
832 }
833}
834
836 Vector &B,
837 int copy_interior) const
838{
839 ReduceSolution(x, X);
840 if (!parallel)
841 {
842 if (!P)
843 {
844 S_e->AddMult(X,*y,-1.);
845 S->PartMult(ess_rtdof_list,X,*y);
846 B.MakeRef(*y, 0, y->Size());
847 }
848 else
849 {
850 B.SetSize(P->Width());
851 P->MultTranspose(*y, B);
852 S_e->AddMult(X,B,-1.);
853 S->PartMult(ess_rtdof_list,X,B);
854 }
855 }
856 else
857 {
858#ifdef MFEM_USE_MPI
859 B.SetSize(pP->Width());
860 pP->MultTranspose(*y,B);
861
862 Vector tmp(B.Size());
863 pS_e->Mult(X,tmp);
864 B-=tmp;
865 for (int j = 0; j<rblocks; j++)
866 {
867 if (!ess_tdofs[j]->Size()) { continue; }
868 HypreParMatrix *Ah = (HypreParMatrix *)(&pS->GetBlock(j,j));
869 Vector diag;
870 Ah->GetDiag(diag);
871 for (int i = 0; i < ess_tdofs[j]->Size(); i++)
872 {
873 int tdof = (*ess_tdofs[j])[i];
874 int gdof = tdof + rtdof_offsets[j];
875 B(gdof) = diag(tdof)*X(gdof);
876 }
877 }
878#endif
879 }
880 if (!copy_interior) { X.SetSubVectorComplement(ess_rtdof_list, 0.0); }
881}
882
883
885 Vector &sol) const
886{
887
888 const int nrdofs = rdof_offsets.Last();
889 const int nrtdofs = rtdof_offsets.Last();
890 MFEM_VERIFY(sc_sol.Size() == nrtdofs, "'sc_sol' has incorrect size");
891
892 Vector sol_r;
893 if (!parallel)
894 {
895 if (!P)
896 {
897 sol_r.SetDataAndSize(sc_sol.GetData(), sc_sol.Size());
898 }
899 else
900 {
901 sol_r.SetSize(nrdofs);
902 P->Mult(sc_sol, sol_r);
903 }
904 }
905 else
906 {
907#ifdef MFEM_USE_MPI
908 sol_r.SetSize(nrdofs);
909 pP->Mult(sc_sol, sol_r);
910#endif
911 }
912
913 if (rdof_offsets.Last() == dof_offsets.Last())
914 {
915 sol = sol_r;
916 return;
917 }
918 else
919 {
920 sol.SetSize(dof_offsets.Last());
921 }
922
923 Vector lsr; // element (local) sc solution vector
924 Vector lsi; // element (local) interior solution vector
925 const int NE = mesh->GetNE();
926
927 Array<int> trace_vdofs;
928 Array<int> vdofs;
929 Array<int> tr_offsets;
930 Vector lsol;
931 for (int iel = 0; iel < NE; iel++)
932 {
933 lsol.SetSize(lmat[iel]->Width() + lmat[iel]->Height());
934 GetReducedElementVDofs(iel, trace_vdofs);
935
936 lsr.SetSize(trace_vdofs.Size());
937 sol_r.GetSubVector(trace_vdofs, lsr);
938
939 // complete the interior dofs
940 lsi.SetSize(lmat[iel]->Height());
941 lmat[iel]->Mult(lsr,lsi);
942 lsi.Neg();
943 lsi+=*lvec[iel];
944
945 Array<int> tr_idx,int_idx,idx_offs;
946 GetReducedElementIndicesAndOffsets(iel,tr_idx, int_idx, idx_offs);
947 lsol.SetSubVector(tr_idx,lsr);
948
949 lsol.SetSubVector(int_idx,lsi);
950
951 GetElementVDofs(iel, vdofs);
952 sol.SetSubVector(vdofs,lsol);
953
954 }
955
956}
957
959{
960 delete S_e; S_e = nullptr;
961 delete S; S=nullptr;
962 delete y; y=nullptr;
963
964 if (P) { delete P; } P=nullptr;
965 if (R) { delete R; } R=nullptr;
966
967#ifdef MFEM_USE_MPI
968 if (parallel)
969 {
970 delete pS; pS=nullptr;
971 delete pS_e; pS_e=nullptr;
972 for (int i = 0; i<rblocks; i++)
973 {
974 delete ess_tdofs[i];
975 }
976 delete pP; pP=nullptr;
977 }
978#endif
979
980 for (int i=0; i<lmat.Size(); i++)
981 {
982 delete lmat[i]; lmat[i] = nullptr;
983 delete lvec[i]; lvec[i] = nullptr;
984 }
985}
986
987} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
Definition array.hpp:882
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
T * GetData()
Returns the data.
Definition array.hpp:118
T & Last()
Return the last element in the array.
Definition array.hpp:802
virtual void AddMult(const Vector &x, Vector &y, const real_t val=1.) const
Matrix-Vector Multiplication y = y + val*A*x.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Partial matrix vector multiplication of (*this) with x involving only the rows given by rows....
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
int NumColBlocks() const
Return the number of column blocks.
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
void EliminateRowCols(const Array< int > &vdofs, BlockMatrix *Ae, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the rows and columns corresponding to the entries in vdofs + save the eliminated entries in...
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
int owns_blocks
If owns_blocks the SparseMatrix objects Aij will be deallocated.
virtual void Finalize(int skip_zeros=1)
Finalize all the submatrices.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A'*x.
int NumRowBlocks() const
Return the number of row blocks.
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
void FormSystemMatrix(Operator::DiagonalPolicy diag_policy)
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
BlockStaticCondensation(Array< FiniteElementSpace * > &fes_)
void ParallelAssemble(BlockMatrix *m)
void ReduceSystem(Vector &x, Vector &X, Vector &B, int copy_interior=0) const
Set the reduced solution X and r.h.s B vectors from the full linear system solution x and r....
void ComputeSolution(const Vector &sc_sol, Vector &sol) const
void ReduceSolution(const Vector &sol, Vector &sc_sol) const
void Finalize(int skip_zeros=0)
Finalize the construction of the Schur complement matrix.
void AssembleReducedSystem(int el, DenseMatrix &elmat, Vector &elvect)
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void GetSubMatrix(const Array< int > &idx, DenseMatrix &A) const
void TransformDual(real_t *v) const
Definition doftrans.cpp:81
@ DISCONTINUOUS
Field is discontinuous across element interfaces.
Definition fe_coll.hpp:48
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
Definition fespace.cpp:667
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition fespace.cpp:648
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition hypre.cpp:1557
void EliminateRows(const Array< int > &rows)
Eliminate rows from the diagonal and off-diagonal blocks of the matrix.
Definition hypre.cpp:2395
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition hypre.cpp:2356
HypreParMatrix * EliminateCols(const Array< int > &cols)
Definition hypre.cpp:2381
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1438
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition mesh.cpp:7201
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition mesh.cpp:6985
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition operator.hpp:48
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
HYPRE_BigInt GlobalVSize() const
Definition pfespace.hpp:283
HYPRE_BigInt * GetDofOffsets() const
Definition pfespace.hpp:281
Class for parallel meshes.
Definition pmesh.hpp:34
Data type sparse matrix.
Definition sparsemat.hpp:51
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:175
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition vector.cpp:670
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
Definition vector.cpp:739
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
int dim
Definition ex24.cpp:53
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
void TransformDual(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition doftrans.cpp:160
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
float real_t
Definition config.hpp:43