MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
complexstaticcond.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 "complexstaticcond.hpp"
13
14namespace mfem
15{
16
19 fes_)
20{
21 SetSpaces(fes_);
22
23 Array<int> rvdofs;
24 Array<int> vdofs;
25 Array<int> rdof_edof0;
26 for (int k = 0; k<nblocks; k++)
27 {
28 if (!tr_fes[k]) { continue; }
29 rdof_edof0.SetSize(tr_fes[k]->GetVSize());
30 for (int i = 0; i < mesh->GetNE(); i++)
31 {
32 fes[k]->GetElementVDofs(i, vdofs);
33 tr_fes[k]->GetElementVDofs(i, rvdofs);
34 const int vdim = fes[k]->GetVDim();
35 const int nsd = vdofs.Size()/vdim;
36 const int nsrd = rvdofs.Size()/vdim;
37 for (int vd = 0; vd < vdim; vd++)
38 {
39 for (int j = 0; j < nsrd; j++)
40 {
41 int rvdof = rvdofs[j+nsrd*vd];
42 int vdof = vdofs[j+nsd*vd];
43 if (rvdof < 0)
44 {
45 rvdof = -1-rvdof;
46 vdof = -1-vdof;
47 }
48 MFEM_ASSERT(vdof >= 0, "incompatible volume and trace FE spaces");
49 rdof_edof0[rvdof] = vdof + dof_offsets[k];
50 }
51 }
52 }
53 rdof_edof.Append(rdof_edof0);
54 }
55}
56
57void ComplexBlockStaticCondensation::SetSpaces(Array<FiniteElementSpace*> &
58 fes_)
59{
60#ifdef MFEM_USE_MPI
61 ParMesh *pmesh = nullptr;
62 parallel = false;
63 if (dynamic_cast<ParFiniteElementSpace *>(fes_[0]))
64 {
65 parallel = true;
66 }
67#else
68 parallel = false;
69#endif
70 fes=fes_;
71 nblocks = fes.Size();
72 rblocks = 0;
73 tr_fes.SetSize(nblocks);
74 mesh = fes[0]->GetMesh();
75
76 IsTraceSpace.SetSize(nblocks);
77 const FiniteElementCollection * fec;
78 for (int i = 0; i < nblocks; i++)
79 {
80 fec = fes[i]->FEColl();
81 IsTraceSpace[i] =
82 (dynamic_cast<const H1_Trace_FECollection*>(fec) ||
83 dynamic_cast<const ND_Trace_FECollection*>(fec) ||
84 dynamic_cast<const RT_Trace_FECollection*>(fec));
85#ifdef MFEM_USE_MPI
86 if (parallel)
87 {
88 pmesh = dynamic_cast<ParMesh *>(mesh);
89 tr_fes[i] = (fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) ?
90 nullptr : (IsTraceSpace[i]) ? fes[i] :
91 new ParFiniteElementSpace(pmesh, fec->GetTraceCollection(), fes[i]->GetVDim(),
92 fes[i]->GetOrdering());
93 }
94 else
95 {
96 tr_fes[i] = (fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) ?
97 nullptr : (IsTraceSpace[i]) ? fes[i] :
98 new FiniteElementSpace(mesh, fec->GetTraceCollection(), fes[i]->GetVDim(),
99 fes[i]->GetOrdering());
100 }
101#else
102 // skip if it's an L2 space (no trace space to construct)
103 tr_fes[i] = (fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) ?
104 nullptr : (IsTraceSpace[i]) ? fes[i] :
105 new FiniteElementSpace(mesh, fec->GetTraceCollection(), fes[i]->GetVDim(),
106 fes[i]->GetOrdering());
107#endif
108 if (tr_fes[i]) { rblocks++; }
109 }
110 if (parallel)
111 {
112 ess_tdofs.SetSize(rblocks);
113 for (int i = 0; i<rblocks; i++)
114 {
115 ess_tdofs[i] = new Array<int>();
116 }
117 }
118 Init();
119}
120
121void ComplexBlockStaticCondensation::ComputeOffsets()
122{
123 dof_offsets.SetSize(nblocks+1);
124 tdof_offsets.SetSize(nblocks+1);
125 dof_offsets[0] = 0;
126 tdof_offsets[0] = 0;
127
128 rdof_offsets.SetSize(rblocks+1);
129 rtdof_offsets.SetSize(rblocks+1);
130 rdof_offsets[0] = 0;
131 rtdof_offsets[0] = 0;
132
133 int j=0;
134 for (int i =0; i<nblocks; i++)
135 {
136 dof_offsets[i+1] = fes[i]->GetVSize();
137 tdof_offsets[i+1] = fes[i]->GetTrueVSize();
138 if (tr_fes[i])
139 {
140 rdof_offsets[j+1] = tr_fes[i]->GetVSize();
141 rtdof_offsets[j+1] = tr_fes[i]->GetTrueVSize();
142 j++;
143 }
144 }
145 rdof_offsets.PartialSum();
146 rtdof_offsets.PartialSum();
147 dof_offsets.PartialSum();
148 tdof_offsets.PartialSum();
149}
150
151void ComplexBlockStaticCondensation::Init()
152{
153 lmat.SetSize(mesh->GetNE());
154 lvec.SetSize(mesh->GetNE());
155 for (int i = 0; i < mesh->GetNE(); i++)
156 {
157 lmat[i] = nullptr;
158 lvec[i] = nullptr;
159 }
160
161 ComputeOffsets();
162
163 S_r = new BlockMatrix(rdof_offsets);
164 S_r->owns_blocks = 1;
165 S_i = new BlockMatrix(rdof_offsets);
166 S_i->owns_blocks = 1;
167
168 for (int i = 0; i<S_r->NumRowBlocks(); i++)
169 {
170 int h = rdof_offsets[i+1] - rdof_offsets[i];
171 for (int j = 0; j<S_r->NumColBlocks(); j++)
172 {
173 int w = rdof_offsets[j+1] - rdof_offsets[j];
174 S_r->SetBlock(i,j,new SparseMatrix(h, w));
175 S_i->SetBlock(i,j,new SparseMatrix(h, w));
176 }
177 }
178
179 y = new Vector(2*rdof_offsets.Last());
180 *y=0.;
181 y_r = new BlockVector(*y, rdof_offsets);
182 y_i = new BlockVector(*y, rdof_offsets.Last(), rdof_offsets);
183}
184
185void ComplexBlockStaticCondensation::GetReduceElementIndicesAndOffsets(int el,
186 Array<int> & trace_ldofs,
187 Array<int> & interior_ldofs,
188 Array<int> & offsets) const
189{
190 int dim = mesh->Dimension();
191 offsets.SetSize(tr_fes.Size()+1); offsets = 0;
192 Array<int> dofs;
193 Array<int> faces, ori;
194 if (dim == 1)
195 {
196 mesh->GetElementVertices(el, faces);
197 }
198 else if (dim == 2)
199 {
200 mesh->GetElementEdges(el, faces, ori);
201 }
202 else if (dim == 3)
203 {
204 mesh->GetElementFaces(el,faces,ori);
205 }
206 else
207 {
208 MFEM_ABORT("ComplexBlockStaticCondensation::GetReduceElementIndicesAndOffsets: "
209 "dim > 3 not supported");
210 }
211 int numfaces = faces.Size();
212
213 trace_ldofs.SetSize(0);
214 interior_ldofs.SetSize(0);
215 // construct Array of bubble dofs to be extracted
216 int skip=0;
217 Array<int> tr_dofs;
218 Array<int> int_dofs;
219 for (int i = 0; i<tr_fes.Size(); i++)
220 {
221 int td = 0;
222 int ndof;
223 // if it's an L2 space (bubbles)
224 if (!tr_fes[i])
225 {
226 ndof = fes[i]->GetVDim()*fes[i]->GetFE(el)->GetDof();
227 td = 0;
228 }
229 else if (IsTraceSpace[i])
230 {
231 for (int iface = 0; iface < numfaces; iface++)
232 {
233 td += fes[i]->GetVDim()*fes[i]->GetFaceElement(faces[iface])->GetDof();
234 }
235 ndof = td;
236 }
237 else
238 {
239 Array<int> trace_dofs;
240 ndof = fes[i]->GetVDim()*fes[i]->GetFE(el)->GetDof();
241 tr_fes[i]->GetElementVDofs(el, trace_dofs);
242 td = trace_dofs.Size(); // number of trace dofs
243 }
244 offsets[i+1] = td;
245 tr_dofs.SetSize(td);
246 int_dofs.SetSize(ndof - td);
247 for (int j = 0; j<td; j++)
248 {
249 tr_dofs[j] = skip + j;
250 }
251 for (int j = 0; j<ndof-td; j++)
252 {
253 int_dofs[j] = skip + td + j;
254 }
255 skip+=ndof;
256
257 trace_ldofs.Append(tr_dofs);
258 interior_ldofs.Append(int_dofs);
259 }
260 offsets.PartialSum();
261}
262
263
264void ComplexBlockStaticCondensation::GetReduceElementVDofs(int el,
265 Array<int> & rdofs) const
266{
267 Array<int> faces, ori;
268 int dim = mesh->Dimension();
269 if (dim == 1)
270 {
271 mesh->GetElementVertices(el, faces);
272 }
273 else if (dim == 2)
274 {
275 mesh->GetElementEdges(el, faces, ori);
276 }
277 else if (dim == 3)
278 {
279 mesh->GetElementFaces(el,faces,ori);
280 }
281 else
282 {
283 MFEM_ABORT("ComplexBlockStaticCondensation::GetReduceElementVDofs: "
284 "dim > 3 not supported");
285 }
286 int numfaces = faces.Size();
287 rdofs.SetSize(0);
288 int skip = 0;
289 for (int i = 0; i<tr_fes.Size(); i++)
290 {
291 if (!tr_fes[i]) { continue; }
292 Array<int> vdofs;
293 if (IsTraceSpace[i])
294 {
295 Array<int> face_vdofs;
296 for (int k = 0; k < numfaces; k++)
297 {
298 int iface = faces[k];
299 tr_fes[i]->GetFaceVDofs(iface, face_vdofs);
300 vdofs.Append(face_vdofs);
301 }
302 }
303 else
304 {
305 tr_fes[i]->GetElementVDofs(el, vdofs);
306 }
307 for (int j=0; j<vdofs.Size(); j++)
308 {
309 vdofs[j] = (vdofs[j]>=0) ? vdofs[j]+rdof_offsets[skip] :
310 vdofs[j]-rdof_offsets[skip];
311 }
312 skip++;
313 rdofs.Append(vdofs);
314 }
315}
316void ComplexBlockStaticCondensation::GetElementVDofs(int el,
317 Array<int> & vdofs) const
318{
319 Array<int> faces, ori;
320 int dim = mesh->Dimension();
321 if (dim == 1)
322 {
323 mesh->GetElementVertices(el, faces);
324 }
325 else if (dim == 2)
326 {
327 mesh->GetElementEdges(el, faces, ori);
328 }
329 else if (dim == 3)
330 {
331 mesh->GetElementFaces(el,faces,ori);
332 }
333 else
334 {
335 MFEM_ABORT("ComplexBlockStaticCondensation::GetElementVDofs: "
336 "dim > 3 not supported");
337 }
338 int numfaces = faces.Size();
339 vdofs.SetSize(0);
340 for (int i = 0; i<tr_fes.Size(); i++)
341 {
342 Array<int> dofs;
343 if (IsTraceSpace[i])
344 {
345 Array<int> face_vdofs;
346 for (int k = 0; k < numfaces; k++)
347 {
348 int iface = faces[k];
349 fes[i]->GetFaceVDofs(iface, face_vdofs);
350 dofs.Append(face_vdofs);
351 }
352 }
353 else
354 {
355 fes[i]->GetElementVDofs(el, dofs);
356 }
357 for (int j=0; j<dofs.Size(); j++)
358 {
359 dofs[j] = (dofs[j]>=0) ? dofs[j]+dof_offsets[i] :
360 dofs[j]-dof_offsets[i];
361 }
362 vdofs.Append(dofs);
363 }
364}
365
366
367ComplexDenseMatrix * ComplexBlockStaticCondensation::GetLocalShurComplement(
368 int el,
369 const Array<int> & tr_idx, const Array<int> & int_idx,
370 const ComplexDenseMatrix & elmat,
371 const Vector & elvect_real, const Vector & elvect_imag,
372 Vector & rvect_real, Vector & rvect_imag)
373{
374 int rdofs = tr_idx.Size();
375 int idofs = int_idx.Size();
376 MFEM_VERIFY(idofs != 0, "Number of interior dofs is zero");
377 MFEM_VERIFY(rdofs != 0, "Number of interface dofs is zero");
378
379 DenseMatrix A_tt_real, A_ti_real, A_it_real, A_ii_real;
380 DenseMatrix A_tt_imag, A_ti_imag, A_it_imag, A_ii_imag;
381
382
383 Vector yt(2*rdofs);
384 Vector yi(2*idofs);
385
386 Vector yt_real(yt, 0,rdofs);
387 Vector yt_imag(yt, rdofs, rdofs);
388
389 Vector yi_real(yi, 0, idofs);
390 Vector yi_imag(yi,idofs, idofs);
391
392 // real part of Matrix and vectors
393 elmat.real().GetSubMatrix(tr_idx,A_tt_real);
394 elmat.real().GetSubMatrix(tr_idx,int_idx, A_ti_real);
395 elmat.real().GetSubMatrix(int_idx, tr_idx, A_it_real);
396 elmat.real().GetSubMatrix(int_idx, A_ii_real);
397
398 elvect_real.GetSubVector(tr_idx, yt_real);
399 elvect_real.GetSubVector(int_idx, yi_real);
400
401 // imag part of Matrix and vectors
402 elmat.imag().GetSubMatrix(tr_idx,A_tt_imag);
403 elmat.imag().GetSubMatrix(tr_idx,int_idx, A_ti_imag);
404 elmat.imag().GetSubMatrix(int_idx, tr_idx, A_it_imag);
405 elmat.imag().GetSubMatrix(int_idx, A_ii_imag);
406
407 elvect_imag.GetSubVector(tr_idx, yt_imag);
408 elvect_imag.GetSubVector(int_idx, yi_imag);
409
410 // construct complex
411 ComplexDenseMatrix A_tt(&A_tt_real,&A_tt_imag,false,false);
412 ComplexDenseMatrix A_ti(&A_ti_real,&A_ti_imag,false,false);
413 ComplexDenseMatrix A_it(&A_it_real,&A_it_imag,false,false);
414 ComplexDenseMatrix A_ii(&A_ii_real,&A_ii_imag,false,false);
415
416 ComplexDenseMatrix * invA_ii = A_ii.ComputeInverse();
417
418 // LHS
419 lmat[el] = mfem::Mult(*invA_ii,A_it);
420 ComplexDenseMatrix * rmat = mfem::Mult(A_ti,*lmat[el]);
421 rmat->real().Neg();
422 rmat->imag().Neg();
423 rmat->real().Add(1., A_tt.real());
424 rmat->imag().Add(1., A_tt.imag());
425
426 // RHS
427 lvec[el] = new Vector(2*idofs);
428 invA_ii->Mult(yi,*lvec[el]);
429 delete invA_ii;
430
431 Vector rvect(2*rdofs);
432 A_ti.Mult(*lvec[el], rvect);
433 rvect_real.SetSize(rdofs);
434 rvect_imag.SetSize(rdofs);
435 for (int i = 0; i<rdofs; i++)
436 {
437 rvect_real(i) = yt_real(i) - rvect(i);
438 rvect_imag(i) = yt_imag(i) - rvect(i+rdofs);
439 }
440 return rmat;
441}
442
443
445 ComplexDenseMatrix &elmat,
446 Vector & elvect_r, Vector & elvect_i)
447{
448 // Get Shur Complement
449 Array<int> tr_idx, int_idx;
450 Array<int> offsets;
451 // Get local element idx and offsets for global assembly
452 GetReduceElementIndicesAndOffsets(el, tr_idx,int_idx, offsets);
453
454 ComplexDenseMatrix *rmat = nullptr;
455 Vector rvec_real, *rvecptr_real;
456 Vector rvec_imag, *rvecptr_imag;
457 // Extract the reduced matrices based on tr_idx and int_idx
458 if (int_idx.Size()!=0)
459 {
460 rmat = GetLocalShurComplement(el,tr_idx,int_idx, elmat, elvect_r, elvect_i,
461 rvec_real,rvec_imag);
462 rvecptr_real = &rvec_real;
463 rvecptr_imag = &rvec_imag;
464 }
465 else
466 {
467 rmat = &elmat;
468 rvecptr_real = &elvect_r;
469 rvecptr_imag = &elvect_i;
470 }
471
472 // Assemble global mat and rhs
473 DofTransformation * doftrans_i, *doftrans_j;
474
475
476 Array<int> faces, ori;
477 int dim = mesh->Dimension();
478 if (dim == 1)
479 {
480 mesh->GetElementVertices(el, faces);
481 }
482 else if (dim == 2)
483 {
484 mesh->GetElementEdges(el, faces, ori);
485 }
486 else if (dim == 3)
487 {
488 mesh->GetElementFaces(el,faces,ori);
489 }
490 else
491 {
492 MFEM_ABORT("ComplexBlockStaticCondensation::AssembleReducedSystem: "
493 "dim > 3 not supported");
494 }
495 int numfaces = faces.Size();
496
497 int skip_i=0;
498 for (int i = 0; i<tr_fes.Size(); i++)
499 {
500 if (!tr_fes[i]) { continue; }
501 Array<int> vdofs_i;
502 doftrans_i = nullptr;
503 if (IsTraceSpace[i])
504 {
505 Array<int> face_vdofs;
506 for (int k = 0; k < numfaces; k++)
507 {
508 int iface = faces[k];
509 tr_fes[i]->GetFaceVDofs(iface, face_vdofs);
510 vdofs_i.Append(face_vdofs);
511 }
512 }
513 else
514 {
515 doftrans_i = tr_fes[i]->GetElementVDofs(el, vdofs_i);
516 }
517 int skip_j=0;
518 for (int j = 0; j<tr_fes.Size(); j++)
519 {
520 if (!tr_fes[j]) { continue; }
521 Array<int> vdofs_j;
522 doftrans_j = nullptr;
523
524 if (IsTraceSpace[j])
525 {
526 Array<int> face_vdofs;
527 for (int k = 0; k < numfaces; k++)
528 {
529 int iface = faces[k];
530 tr_fes[j]->GetFaceVDofs(iface, face_vdofs);
531 vdofs_j.Append(face_vdofs);
532 }
533 }
534 else
535 {
536 doftrans_j = tr_fes[j]->GetElementVDofs(el, vdofs_j);
537 }
538
539 DenseMatrix Ae_r, Ae_i;
540 rmat->real().GetSubMatrix(offsets[i],offsets[i+1],
541 offsets[j],offsets[j+1], Ae_r);
542 rmat->imag().GetSubMatrix(offsets[i],offsets[i+1],
543 offsets[j],offsets[j+1], Ae_i);
544 if (doftrans_i || doftrans_j)
545 {
546 TransformDual(doftrans_i, doftrans_j, Ae_r);
547 TransformDual(doftrans_i, doftrans_j, Ae_i);
548 }
549 S_r->GetBlock(skip_i,skip_j).AddSubMatrix(vdofs_i,vdofs_j, Ae_r);
550 S_i->GetBlock(skip_i,skip_j).AddSubMatrix(vdofs_i,vdofs_j, Ae_i);
551 skip_j++;
552 }
553
554 // assemble rhs
555 Vector vec1_r(*rvecptr_real, offsets[i], offsets[i+1]-offsets[i]);
556 Vector vec1_i(*rvecptr_imag, offsets[i], offsets[i+1]-offsets[i]);
557 // ref subvector
558 if (doftrans_i)
559 {
560 doftrans_i->TransformDual(vec1_r);
561 doftrans_i->TransformDual(vec1_i);
562 }
563 y_r->GetBlock(skip_i).AddElementVector(vdofs_i,vec1_r);
564 y_i->GetBlock(skip_i).AddElementVector(vdofs_i,vec1_i);
565 skip_i++;
566 }
567 if (int_idx.Size()!=0) { delete rmat; }
568}
569
570void ComplexBlockStaticCondensation::BuildProlongation()
571{
572 P = new BlockMatrix(rdof_offsets, rtdof_offsets);
573 R = new BlockMatrix(rtdof_offsets, rdof_offsets);
574 P->owns_blocks = 0;
575 R->owns_blocks = 0;
576 int skip = 0;
577 for (int i = 0; i<nblocks; i++)
578 {
579 if (!tr_fes[i]) { continue; }
580 const SparseMatrix *P_ = tr_fes[i]->GetConformingProlongation();
581 if (P_)
582 {
583 const SparseMatrix *R_ = tr_fes[i]->GetRestrictionMatrix();
584 P->SetBlock(skip,skip,const_cast<SparseMatrix*>(P_));
585 R->SetBlock(skip,skip,const_cast<SparseMatrix*>(R_));
586 }
587 skip++;
588 }
589}
590
591#ifdef MFEM_USE_MPI
592void ComplexBlockStaticCondensation::BuildParallelProlongation()
593{
594 MFEM_VERIFY(parallel, "BuildParallelProlongation: wrong code path");
595 pP = new BlockOperator(rdof_offsets, rtdof_offsets);
596 R = new BlockMatrix(rtdof_offsets, rdof_offsets);
597 pP->owns_blocks = 0;
598 R->owns_blocks = 0;
599 int skip = 0;
600 for (int i = 0; i<nblocks; i++)
601 {
602 if (!tr_fes[i]) { continue; }
603 const HypreParMatrix *P_ =
604 dynamic_cast<ParFiniteElementSpace *>(tr_fes[i])->Dof_TrueDof_Matrix();
605 if (P_)
606 {
607 const SparseMatrix *R_ = tr_fes[i]->GetRestrictionMatrix();
608 pP->SetBlock(skip,skip,const_cast<HypreParMatrix*>(P_));
609 R->SetBlock(skip,skip,const_cast<SparseMatrix*>(R_));
610 }
611 skip++;
612 }
613}
614
616 BlockMatrix *m_i)
617{
618
619 if (!pP) { BuildParallelProlongation(); }
620
621 pS_r = new BlockOperator(rtdof_offsets);
622 pS_e_r = new BlockOperator(rtdof_offsets);
623 pS_i = new BlockOperator(rtdof_offsets);
624 pS_e_i = new BlockOperator(rtdof_offsets);
625 pS_r->owns_blocks = 1;
626 pS_i->owns_blocks = 1;
627 pS_e_r->owns_blocks = 1;
628 pS_e_i->owns_blocks = 1;
629 HypreParMatrix * A_r = nullptr;
630 HypreParMatrix * A_i = nullptr;
631 HypreParMatrix * PtAP_r = nullptr;
632 HypreParMatrix * PtAP_i = nullptr;
633 int skip_i=0;
634 ParFiniteElementSpace * pfes_i = nullptr;
635 ParFiniteElementSpace * pfes_j = nullptr;
636 for (int i = 0; i<nblocks; i++)
637 {
638 if (!tr_fes[i]) { continue; }
639 pfes_i = dynamic_cast<ParFiniteElementSpace*>(tr_fes[i]);
640 HypreParMatrix * Pi = (HypreParMatrix*)(&pP->GetBlock(skip_i,skip_i));
641 int skip_j=0;
642 for (int j = 0; j<nblocks; j++)
643 {
644 if (!tr_fes[j]) { continue; }
645 if (m_r->IsZeroBlock(skip_i,skip_j)) { continue; }
646 if (skip_i == skip_j)
647 {
648 // Make block diagonal square hypre matrix
649 A_r = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
650 pfes_i->GetDofOffsets(),&m_r->GetBlock(skip_i,skip_i));
651 PtAP_r = RAP(A_r,Pi);
652 delete A_r;
653
654 pS_e_r->SetBlock(skip_i,skip_i,PtAP_r->EliminateRowsCols(*ess_tdofs[skip_i]));
655
656 A_i = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
657 pfes_i->GetDofOffsets(),&m_i->GetBlock(skip_i,skip_i));
658 PtAP_i = RAP(A_i,Pi);
659 delete A_i;
660 pS_e_i->SetBlock(skip_i,skip_j,PtAP_i->EliminateCols(*ess_tdofs[skip_j]));
661 PtAP_i->EliminateRows(*ess_tdofs[skip_i]);
662 }
663 else
664 {
665 pfes_j = dynamic_cast<ParFiniteElementSpace*>(tr_fes[j]);
666 HypreParMatrix * Pj = (HypreParMatrix*)(&pP->GetBlock(skip_j,skip_j));
667 A_r = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
668 pfes_j->GlobalVSize(), pfes_i->GetDofOffsets(),
669 pfes_j->GetDofOffsets(), &m_r->GetBlock(skip_i,skip_j));
670 PtAP_r = RAP(Pi,A_r,Pj);
671 delete A_r;
672 pS_e_r->SetBlock(skip_i,skip_j,PtAP_r->EliminateCols(*ess_tdofs[skip_j]));
673 PtAP_r->EliminateRows(*ess_tdofs[skip_i]);
674
675 A_i = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
676 pfes_j->GlobalVSize(), pfes_i->GetDofOffsets(),
677 pfes_j->GetDofOffsets(), &m_i->GetBlock(skip_i,skip_j));
678 PtAP_i = RAP(Pi,A_i,Pj);
679 delete A_i;
680 pS_e_i->SetBlock(skip_i,skip_j,PtAP_i->EliminateCols(*ess_tdofs[skip_j]));
681 PtAP_i->EliminateRows(*ess_tdofs[skip_i]);
682
683 }
684 pS_r->SetBlock(skip_i,skip_j,PtAP_r);
685 pS_i->SetBlock(skip_i,skip_j,PtAP_i);
686 skip_j++;
687 }
688 skip_i++;
689 }
690}
691
692#endif
693
694
695void ComplexBlockStaticCondensation::ConformingAssemble(int skip_zeros)
696{
697 Finalize(0);
698 if (!P) { BuildProlongation(); }
699
700 BlockMatrix * Pt = Transpose(*P);
701 BlockMatrix * PtA_r = mfem::Mult(*Pt, *S_r);
702 BlockMatrix * PtA_i = mfem::Mult(*Pt, *S_i);
703 delete S_r;
704 delete S_i;
705 if (S_e_r)
706 {
707 BlockMatrix *PtAe_r = mfem::Mult(*Pt, *S_e_r);
708 BlockMatrix *PtAe_i = mfem::Mult(*Pt, *S_e_i);
709 delete S_e_r;
710 delete S_e_i;
711 S_e_r = PtAe_r;
712 S_e_i = PtAe_i;
713 }
714 delete Pt;
715 S_r = mfem::Mult(*PtA_r, *P);
716 S_i = mfem::Mult(*PtA_i, *P);
717 delete PtA_r;
718 delete PtA_i;
719
720 if (S_e_r)
721 {
722 BlockMatrix *PtAeP_r = mfem::Mult(*S_e_r, *P);
723 BlockMatrix *PtAeP_i = mfem::Mult(*S_e_i, *P);
724 S_e_r = PtAeP_r;
725 S_e_i = PtAeP_i;
726 }
727 height = 2*S_r->Height();
728 width = 2*S_r->Width();
729}
730
732{
733 if (S_r)
734 {
735 S_r->Finalize(skip_zeros);
736 S_i->Finalize(skip_zeros);
737 }
738 if (S_e_r)
739 {
740 S_e_r->Finalize(skip_zeros);
741 S_e_i->Finalize(skip_zeros);
742 }
743}
744
746 diag_policy)
747{
748
749 if (!parallel)
750 {
751 if (!S_e_r)
752 {
753 bool conforming = true;
754 for (int i = 0; i<nblocks; i++)
755 {
756 if (!tr_fes[i]) { continue; }
757 const SparseMatrix *P_ = tr_fes[i]->GetConformingProlongation();
758 if (P_)
759 {
760 conforming = false;
761 break;
762 }
763 }
764 if (!conforming) { ConformingAssemble(0); }
765 const int remove_zeros = 0;
766 EliminateReducedTrueDofs(ess_rtdof_list, diag_policy);
767 Finalize(remove_zeros);
768 }
769 }
770 else
771 {
772#ifdef MFEM_USE_MPI
773 FillEssTdofLists(ess_rtdof_list);
774 if (S_r)
775 {
776 const int remove_zeros = 0;
777 Finalize(remove_zeros);
778 ParallelAssemble(S_r, S_i);
779 delete S_r; S_r=nullptr;
780 delete S_i; S_i=nullptr;
781 delete S_e_r; S_e_r = nullptr;
782 delete S_e_i; S_e_i = nullptr;
783 }
784#endif
785 }
786}
787
788void ComplexBlockStaticCondensation::ConvertMarkerToReducedTrueDofs(
789 Array<int> & tdof_marker,
790 Array<int> & rtdof_marker)
791{
792 // convert tdof_marker to dof_marker
793 rtdof_marker.SetSize(0);
794 Array<int> tdof_marker0;
795 Array<int> dof_marker0;
796 Array<int> dof_marker;
797
798 for (int i = 0; i<nblocks; i++)
799 {
800 tdof_marker0.MakeRef(&tdof_marker[tdof_offsets[i]],
801 tdof_offsets[i+1]-tdof_offsets[i]);
802 const SparseMatrix * R_ = fes[i]->GetRestrictionMatrix();
803 if (!R_)
804 {
805 dof_marker0.MakeRef(tdof_marker0);
806 }
807 else
808 {
809 dof_marker0.SetSize(fes[i]->GetVSize());
810 R_->BooleanMultTranspose(tdof_marker0, dof_marker0);
811 }
812 dof_marker.Append(dof_marker0);
813 }
814
815 int rdofs = rdof_edof.Size();
816 Array<int> rdof_marker(rdofs);
817
818 for (int i = 0; i < rdofs; i++)
819 {
820 rdof_marker[i] = dof_marker[rdof_edof[i]];
821 }
822
823 // convert rdof_marker to rtdof_marker
824 Array<int> rtdof_marker0;
825 Array<int> rdof_marker0;
826 int k=0;
827 for (int i = 0; i<nblocks; i++)
828 {
829 if (!tr_fes[i]) { continue; }
830 rdof_marker0.MakeRef(&rdof_marker[rdof_offsets[k]],
831 rdof_offsets[k+1]-rdof_offsets[k]);
832 const SparseMatrix *tr_R = tr_fes[i]->GetRestrictionMatrix();
833 if (!tr_R)
834 {
835 rtdof_marker0.MakeRef(rdof_marker0);
836 }
837 else
838 {
839 rtdof_marker0.SetSize(tr_fes[i]->GetTrueVSize());
840 tr_R->BooleanMult(rdof_marker0, rtdof_marker0);
841 }
842 rtdof_marker.Append(rtdof_marker0);
843 k++;
844 }
845}
846
847void ComplexBlockStaticCondensation::FillEssTdofLists(const Array<int> &
848 ess_tdof_list)
849{
850 int j;
851 for (int i = 0; i<ess_tdof_list.Size(); i++)
852 {
853 int tdof = ess_tdof_list[i];
854 for (j = 0; j < rblocks; j++)
855 {
856 if (rtdof_offsets[j+1] > tdof) { break; }
857 }
858 ess_tdofs[j]->Append(tdof-rtdof_offsets[j]);
859 }
860}
861
863 &ess_tdof_list)
864{
865 Array<int> tdof_marker;
866 Array<int> rtdof_marker;
867 FiniteElementSpace::ListToMarker(ess_tdof_list,tdof_offsets.Last(),tdof_marker);
868 ConvertMarkerToReducedTrueDofs(tdof_marker, rtdof_marker);
869 FiniteElementSpace::MarkerToList(rtdof_marker,ess_rtdof_list);
870}
871
873 &ess_rtdof_list_,
875{
876
877 MFEM_VERIFY(!parallel, "EliminateReducedTrueDofs::Wrong code path");
878
879 if (S_e_r == NULL)
880 {
881 Array<int> offsets;
882
883 offsets.MakeRef( (P) ? rtdof_offsets : rdof_offsets);
884
885 S_e_r = new BlockMatrix(offsets);
886 S_e_i = new BlockMatrix(offsets);
887 S_e_r->owns_blocks = 1;
888 S_e_i->owns_blocks = 1;
889 for (int i = 0; i<S_e_r->NumRowBlocks(); i++)
890 {
891 int h = offsets[i+1] - offsets[i];
892 for (int j = 0; j<S_e_r->NumColBlocks(); j++)
893 {
894 int w = offsets[j+1] - offsets[j];
895 S_e_r->SetBlock(i,j,new SparseMatrix(h, w));
896 S_e_i->SetBlock(i,j,new SparseMatrix(h, w));
897 }
898 }
899 }
900 S_r->EliminateRowCols(ess_rtdof_list_,S_e_r,dpolicy);
901 S_i->EliminateRowCols(ess_rtdof_list_,S_e_i,
903}
904
906 Vector &sc_sol) const
907{
908 MFEM_ASSERT(sol.Size() == 2*dof_offsets.Last(), "'sol' has incorrect size");
909 const int nrdofs = rdof_offsets.Last();
910
911 Vector sol_r_real;
912 Vector sol_r_imag;
913
914 if (!R)
915 {
916 sc_sol.SetSize(2*nrdofs);
917 sol_r_real.MakeRef(sc_sol, 0, nrdofs);
918 sol_r_imag.MakeRef(sc_sol, nrdofs, nrdofs);
919 }
920 else
921 {
922 sol_r_real.SetSize(nrdofs);
923 sol_r_imag.SetSize(nrdofs);
924 }
925 for (int i = 0; i < nrdofs; i++)
926 {
927 sol_r_real(i) = sol(rdof_edof[i]);
928 sol_r_imag(i) = sol(rdof_edof[i] + dof_offsets.Last());
929 }
930
931 if (R)
932 {
933 int n = R->Height();
934 sc_sol.SetSize(2*n);
935 Vector sc_real(sc_sol, 0, n);
936 Vector sc_imag(sc_sol, n, n);
937
938 // wrap vector into a block vector
939 BlockVector blsol_r_real(sol_r_real,rdof_offsets);
940 BlockVector blsol_r_imag(sol_r_imag,rdof_offsets);
941 R->Mult(blsol_r_real, sc_real);
942 R->Mult(blsol_r_imag, sc_imag);
943 }
944}
945
947 Vector &B,
948 int copy_interior) const
949{
950 ReduceSolution(x, X);
951 Vector X_r(X,0, X.Size()/2);
952 Vector X_i(X, X.Size()/2, X.Size()/2);
953
954 if (!parallel)
955 {
956 if (!P)
957 {
958
959 S_e_r->AddMult(X_r,*y_r,-1.);
960 S_e_i->AddMult(X_i,*y_r,1.);
961 S_e_r->AddMult(X_i,*y_i,-1.);
962 S_e_i->AddMult(X_r,*y_i,-1.);
963
964 S_r->PartMult(ess_rtdof_list,X_r,*y_r);
965 S_r->PartMult(ess_rtdof_list,X_i,*y_i);
966 B.MakeRef(*y, 0, y->Size());
967 }
968 else
969 {
970 B.SetSize(2*P->Width());
971 Vector B_r(B, 0, P->Width());
972 Vector B_i(B, P->Width(), P->Width());
973
974 P->MultTranspose(*y_r, B_r);
975 P->MultTranspose(*y_i, B_i);
976
977 S_e_r->AddMult(X_r,B_r,-1.);
978 S_e_i->AddMult(X_i,B_r,1.);
979 S_e_r->AddMult(X_i,B_i,-1.);
980 S_e_i->AddMult(X_r,B_i,-1.);
981 S_r->PartMult(ess_rtdof_list,X_r,B_r);
982 S_r->PartMult(ess_rtdof_list,X_i,B_i);
983 }
984 }
985 else
986 {
987#ifdef MFEM_USE_MPI
988 int n = pP->Width();
989 B.SetSize(2*n);
990 Vector B_r(B, 0, n);
991 Vector B_i(B, n, n);
992
993 pP->MultTranspose(*y_r,B_r);
994 pP->MultTranspose(*y_i,B_i);
995
996 Vector tmp(B_r.Size());
997 pS_e_r->Mult(X_r,tmp); B_r-=tmp;
998 pS_e_i->Mult(X_i,tmp); B_r+=tmp;
999
1000 pS_e_i->Mult(X_r,tmp); B_i-=tmp;
1001 pS_e_r->Mult(X_i,tmp); B_i-=tmp;
1002
1003 for (int j = 0; j<rblocks; j++)
1004 {
1005 if (!ess_tdofs[j]->Size()) { continue; }
1006 for (int i = 0; i < ess_tdofs[j]->Size(); i++)
1007 {
1008 int tdof = (*ess_tdofs[j])[i];
1009 int gdof = tdof + rtdof_offsets[j];
1010 B_r(gdof) = X_r(gdof);
1011 B_i(gdof) = X_i(gdof);
1012 }
1013 }
1014#endif
1015 }
1016 if (!copy_interior)
1017 {
1018 X_r.SetSubVectorComplement(ess_rtdof_list, 0.0);
1019 X_i.SetSubVectorComplement(ess_rtdof_list, 0.0);
1020 }
1021}
1022
1023
1025 Vector &sol) const
1026{
1027
1028 const int nrdofs = rdof_offsets.Last();
1029 const int nrtdofs = rtdof_offsets.Last();
1030 MFEM_VERIFY(sc_sol.Size() == 2*nrtdofs, "'sc_sol' has incorrect size");
1031
1032 Vector sol_r_real;
1033 Vector sol_r_imag;
1034 if (!parallel)
1035 {
1036 if (!P)
1037 {
1038 sol_r_real.MakeRef(const_cast<Vector &>(sc_sol), 0, sc_sol.Size()/2);
1039 sol_r_imag.MakeRef(const_cast<Vector &>(sc_sol), sc_sol.Size()/2,
1040 sc_sol.Size()/2);
1041 }
1042 else
1043 {
1044 Vector sc_real(const_cast<Vector &>(sc_sol),0, nrtdofs);
1045 Vector sc_imag(const_cast<Vector &>(sc_sol),nrtdofs, nrtdofs);
1046 sol_r_real.SetSize(nrdofs);
1047 sol_r_imag.SetSize(nrdofs);
1048 P->Mult(sc_real, sol_r_real);
1049 P->Mult(sc_imag, sol_r_imag);
1050 }
1051 }
1052 else
1053 {
1054#ifdef MFEM_USE_MPI
1055 Vector sc_real(const_cast<Vector &>(sc_sol),0, nrtdofs);
1056 Vector sc_imag(const_cast<Vector &>(sc_sol),nrtdofs, nrtdofs);
1057 sol_r_real.SetSize(nrdofs);
1058 sol_r_imag.SetSize(nrdofs);
1059 pP->Mult(sc_real, sol_r_real);
1060 pP->Mult(sc_imag, sol_r_imag);
1061#endif
1062 }
1063
1064 sol.SetSize(2*dof_offsets.Last());
1065 Vector sol_real(sol,0,dof_offsets.Last());
1066 Vector sol_imag(sol,dof_offsets.Last(),dof_offsets.Last());
1067
1068 if (rdof_offsets.Last() == dof_offsets.Last())
1069 {
1070 sol_real = sol_r_real;
1071 sol_imag = sol_r_imag;
1072 return;
1073 }
1074
1075 Vector lsr; // element (local) sc solution vector
1076 Vector lsr_real; // element (local) sc solution vector
1077 Vector lsr_imag; // element (local) sc solution vector
1078 Vector lsi; // element (local) interior solution vector
1079 Vector lsi_real; // element (local) interior solution vector
1080 Vector lsi_imag; // element (local) interior solution vector
1081
1082 const int NE = mesh->GetNE();
1083
1084 Array<int> trace_vdofs;
1085 Array<int> vdofs;
1086 Array<int> tr_offsets;
1087 Vector lsol;
1088 Vector lsol_real;
1089 Vector lsol_imag;
1090 for (int iel = 0; iel < NE; iel++)
1091 {
1092 GetReduceElementVDofs(iel, trace_vdofs);
1093
1094 int n = trace_vdofs.Size();
1095 lsr.SetSize(2*n);
1096 lsr_real.MakeRef(lsr, 0, n);
1097 lsr_imag.MakeRef(lsr, n, n);
1098 sol_r_real.GetSubVector(trace_vdofs, lsr_real);
1099 sol_r_imag.GetSubVector(trace_vdofs, lsr_imag);
1100
1101 // complete the interior dofs
1102 int m = lmat[iel]->Height()/2;
1103 lsi.SetSize(2*m);
1104 lsi_real.MakeRef(lsi, 0, m);
1105 lsi_imag.MakeRef(lsi, m, m);
1106 lmat[iel]->Mult(lsr,lsi);
1107 lsi.Neg();
1108 lsi+=*lvec[iel];
1109
1110 Array<int> tr_idx,int_idx,idx_offs;
1111 GetReduceElementIndicesAndOffsets(iel,tr_idx, int_idx, idx_offs);
1112
1113 // complete all the dofs in the element
1114 int k = (lmat[iel]->Width() + lmat[iel]->Height())/2;
1115 lsol.SetSize(2*k);
1116 lsol_real.MakeRef(lsol, 0, k);
1117 lsol_imag.MakeRef(lsol, k, k);
1118
1119 lsol_real.SetSubVector(tr_idx,lsr_real);
1120 lsol_real.SetSubVector(int_idx,lsi_real);
1121 lsol_imag.SetSubVector(tr_idx,lsr_imag);
1122 lsol_imag.SetSubVector(int_idx,lsi_imag);
1123
1124 GetElementVDofs(iel, vdofs);
1125
1126 // complete all the dofs in the global vector
1127 sol_real.SetSubVector(vdofs,lsol_real);
1128 sol_imag.SetSubVector(vdofs,lsol_imag);
1129 }
1130}
1131
1133{
1134 delete S_e_r; S_e_r = nullptr;
1135 delete S_e_i; S_e_i = nullptr;
1136 delete S_r; S_r = nullptr;
1137 delete S_i; S_i = nullptr;
1138 delete S; S=nullptr;
1139 delete y_r; y_r=nullptr;
1140 delete y_i; y_i=nullptr;
1141 delete y; y=nullptr;
1142
1143 if (P) { delete P; } P=nullptr;
1144 if (R) { delete R; } R=nullptr;
1145
1146#ifdef MFEM_USE_MPI
1147 if (parallel)
1148 {
1149 // The Complex Operator (S) is deleted above
1150 delete pS_e_r; pS_e_r=nullptr;
1151 delete pS_e_i; pS_e_i=nullptr;
1152 delete pS_r; pS_r=nullptr;
1153 delete pS_i; pS_i=nullptr;
1154 for (int i = 0; i<rblocks; i++)
1155 {
1156 delete ess_tdofs[i];
1157 }
1158 delete pP; pP = nullptr;
1159 }
1160#endif
1161
1162 for (int i=0; i<lmat.Size(); i++)
1163 {
1164 delete lmat[i]; lmat[i] = nullptr;
1165 delete lvec[i]; lvec[i] = nullptr;
1166 }
1167}
1168
1169}
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
Definition array.hpp:882
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
T & Last()
Return the last element in the array.
Definition array.hpp:802
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.
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
ComplexBlockStaticCondensation(Array< FiniteElementSpace * > &fes_)
void ParallelAssemble(BlockMatrix *m_r, BlockMatrix *m_i)
void Finalize(int skip_zeros=0)
Finalize the construction of the Schur complement matrix.
void AssembleReducedSystem(int el, ComplexDenseMatrix &elmat, Vector &elvect_r, Vector &elvect_i)
void FormSystemMatrix(Operator::DiagonalPolicy diag_policy)
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
void ComputeSolution(const Vector &sc_sol, Vector &sol) const
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 ReduceSolution(const Vector &sol, Vector &sc_sol) const
Specialization of the ComplexOperator built from a pair of Dense Matrices. The purpose of this specia...
virtual DenseMatrix & real()
Real or imaginary part accessor methods.
virtual DenseMatrix & imag()
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 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
@ DIAG_ZERO
Set the diagonal value to zero.
Definition operator.hpp:49
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 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
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)