MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fespace.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// Implementation of FiniteElementSpace
13
14#include "../general/text.hpp"
15#include "../general/forall.hpp"
17#include "fem.hpp"
19
20#include <cmath>
21#include <cstdarg>
22
23using namespace std;
24
25namespace mfem
26{
27
28template <> void Ordering::
29DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
30{
31 // static method
32 int size = dofs.Size();
33 dofs.SetSize(size*vdim);
34 for (int vd = 1; vd < vdim; vd++)
35 {
36 for (int i = 0; i < size; i++)
37 {
38 dofs[i+size*vd] = Map<byNODES>(ndofs, vdim, dofs[i], vd);
39 }
40 }
41}
42
43template <> void Ordering::
44DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
45{
46 // static method
47 int size = dofs.Size();
48 dofs.SetSize(size*vdim);
49 for (int vd = vdim-1; vd >= 0; vd--)
50 {
51 for (int i = 0; i < size; i++)
52 {
53 dofs[i+size*vd] = Map<byVDIM>(ndofs, vdim, dofs[i], vd);
54 }
55 }
56}
57
58
60 : mesh(NULL), fec(NULL), vdim(0), ordering(Ordering::byNODES),
61 ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0),
62 bdofs(NULL),
63 elem_dof(NULL), elem_fos(NULL), bdr_elem_dof(NULL), bdr_elem_fos(NULL),
64 face_dof(NULL),
65 NURBSext(NULL), own_ext(false),
66 cP_is_set(false),
67 Th(Operator::ANY_TYPE),
68 sequence(0), mesh_sequence(0), orders_changed(false), relaxed_hp(false)
69{ }
70
72 Mesh *mesh_,
73 const FiniteElementCollection *fec_)
74{
75 mesh_ = mesh_ ? mesh_ : orig.mesh;
76 fec_ = fec_ ? fec_ : orig.fec;
77
78 NURBSExtension *nurbs_ext = NULL;
79 if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
80 {
81#ifdef MFEM_USE_MPI
82 ParNURBSExtension *pNURBSext =
83 dynamic_cast<ParNURBSExtension *>(orig.NURBSext);
84 if (pNURBSext)
85 {
86 nurbs_ext = new ParNURBSExtension(*pNURBSext);
87 }
88 else
89#endif
90 {
91 nurbs_ext = new NURBSExtension(*orig.NURBSext);
92 }
93 }
94
95 Constructor(mesh_, nurbs_ext, fec_, orig.vdim, orig.ordering);
96}
97
99 const FiniteElementCollection *fec,
100 int vdim, int ordering)
101{ Constructor(mesh, NULL, fec, vdim, ordering); }
102
104 const FiniteElementCollection *fec,
105 int vdim, int ordering)
106{ Constructor(mesh, ext, fec, vdim, ordering); }
107
109 const FiniteElementSpace &fes, const Array<int> *perm)
110{
111 MFEM_VERIFY(cP == NULL, "");
112 MFEM_VERIFY(cR == NULL, "");
113
114 SparseMatrix *perm_mat = NULL, *perm_mat_tr = NULL;
115 if (perm)
116 {
117 // Note: although n and fes.GetVSize() are typically equal, in
118 // variable-order spaces they may differ, since nonconforming edges/faces
119 // my have fictitious DOFs.
120 int n = perm->Size();
121 perm_mat = new SparseMatrix(n, fes.GetVSize());
122 for (int i=0; i<n; ++i)
123 {
124 real_t s;
125 int j = DecodeDof((*perm)[i], s);
126 perm_mat->Set(i, j, s);
127 }
128 perm_mat->Finalize();
129 perm_mat_tr = Transpose(*perm_mat);
130 }
131
132 if (fes.GetConformingProlongation() != NULL)
133 {
134 if (perm) { cP.reset(Mult(*perm_mat, *fes.GetConformingProlongation())); }
135 else { cP.reset(new SparseMatrix(*fes.GetConformingProlongation())); }
136 cP_is_set = true;
137 }
138 else if (perm != NULL)
139 {
140 cP.reset(perm_mat);
141 cP_is_set = true;
142 perm_mat = NULL;
143 }
144 if (fes.GetConformingRestriction() != NULL)
145 {
146 if (perm) { cR.reset(Mult(*fes.GetConformingRestriction(), *perm_mat_tr)); }
147 else { cR.reset(new SparseMatrix(*fes.GetConformingRestriction())); }
148 }
149 else if (perm != NULL)
150 {
151 cR.reset(perm_mat_tr);
152 perm_mat_tr = NULL;
153 }
154
155 delete perm_mat;
156 delete perm_mat_tr;
157}
158
160{
161#ifdef MFEM_USE_MPI
162 MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
163 "Attempting to set serial prolongation operator for "
164 "parallel finite element space.");
165#endif
166
167 if (!cP)
168 {
169 cP = std::unique_ptr<SparseMatrix>(new SparseMatrix(p));
170 }
171 else
172 {
173 *cP = p;
174 }
175 cP_is_set = true;
176}
177
179{
180#ifdef MFEM_USE_MPI
181 MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
182 "Attempting to set serial restriction operator for "
183 "parallel finite element space.");
184#endif
185
186 if (!cR)
187 {
188 cR = std::unique_ptr<SparseMatrix>(new SparseMatrix(r));
189 }
190 else
191 {
192 *cR = r;
193 }
194}
195
197{
198 MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),
199 "Space has not been Updated() after a Mesh change.");
200 MFEM_VERIFY(i >= 0 && i < GetNE(), "Invalid element index");
201 MFEM_VERIFY(p >= 0 && p <= MaxVarOrder, "Order out of range");
202 MFEM_ASSERT(!elem_order.Size() || elem_order.Size() == GetNE(),
203 "Internal error");
204
205 const bool change = elem_order.Size() == 0 || elem_order[i] != p;
206 if (elem_order.Size() == 0) // convert space to variable-order space
207 {
210 }
211
212 if (change)
213 {
214 elem_order[i] = p;
215 orders_changed = true;
216 }
217
218 variableOrder = true;
219}
220
222{
223 MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),
224 "Space has not been Updated() after a Mesh change.");
225 MFEM_VERIFY(i >= 0 && i < GetNE(), "Invalid element index");
226 MFEM_ASSERT(!elem_order.Size() || elem_order.Size() == GetNE(),
227 "Internal error");
228
229 return GetElementOrderImpl(i);
230}
231
233{
234 // (this is an internal version of GetElementOrder without asserts and checks)
235 return elem_order.Size() ? elem_order[i] : fec->GetOrder();
236}
237
238void FiniteElementSpace::GetVDofs(int vd, Array<int>& dofs, int ndofs_) const
239{
240 if (ndofs_ < 0) { ndofs_ = this->ndofs; }
241
243 {
244 for (int i = 0; i < dofs.Size(); i++)
245 {
246 dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs_, vdim, i, vd);
247 }
248 }
249 else
250 {
251 for (int i = 0; i < dofs.Size(); i++)
252 {
253 dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs_, vdim, i, vd);
254 }
255 }
256}
257
258void FiniteElementSpace::DofsToVDofs(Array<int> &dofs, int ndofs_) const
259{
260 if (vdim == 1) { return; }
261 if (ndofs_ < 0) { ndofs_ = this->ndofs; }
262
264 {
265 Ordering::DofsToVDofs<Ordering::byNODES>(ndofs_, vdim, dofs);
266 }
267 else
268 {
269 Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs_, vdim, dofs);
270 }
271}
272
273void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs, int ndofs_) const
274{
275 if (vdim == 1) { return; }
276 if (ndofs_ < 0) { ndofs_ = this->ndofs; }
277
279 {
280 for (int i = 0; i < dofs.Size(); i++)
281 {
282 dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs_, vdim, dofs[i], vd);
283 }
284 }
285 else
286 {
287 for (int i = 0; i < dofs.Size(); i++)
288 {
289 dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs_, vdim, dofs[i], vd);
290 }
291 }
292}
293
294int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
295{
296 if (vdim == 1) { return dof; }
297 if (ndofs_ < 0) { ndofs_ = this->ndofs; }
298
300 {
301 return Ordering::Map<Ordering::byNODES>(ndofs_, vdim, dof, vd);
302 }
303 else
304 {
305 return Ordering::Map<Ordering::byVDIM>(ndofs_, vdim, dof, vd);
306 }
307}
308
309// static function
311{
312 int n = vdofs.Size(), *vdof = vdofs;
313 for (int i = 0; i < n; i++)
314 {
315 int j;
316 if ((j = vdof[i]) < 0)
317 {
318 vdof[i] = -1-j;
319 }
320 }
321}
322
324 DofTransformation &doftrans) const
325{
326 GetElementDofs(i, vdofs, doftrans);
327 DofsToVDofs(vdofs);
328 doftrans.SetVDim(vdim, ordering);
329}
330
333{
335 GetElementVDofs(i, vdofs, DoFTrans);
336 return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
337}
338
340 DofTransformation &doftrans) const
341{
342 GetBdrElementDofs(i, vdofs, doftrans);
343 DofsToVDofs(vdofs);
344 doftrans.SetVDim(vdim, ordering);
345}
346
354
356{
357 GetPatchDofs(i, vdofs);
358 DofsToVDofs(vdofs);
359}
360
362{
363 GetFaceDofs(i, vdofs);
364 DofsToVDofs(vdofs);
365}
366
368{
369 GetEdgeDofs(i, vdofs);
370 DofsToVDofs(vdofs);
371}
372
374{
375 GetVertexDofs(i, vdofs);
376 DofsToVDofs(vdofs);
377}
378
380{
381 GetElementInteriorDofs(i, vdofs);
382 DofsToVDofs(vdofs);
383}
384
386{
387 GetEdgeInteriorDofs(i, vdofs);
388 DofsToVDofs(vdofs);
389}
390
392{
393 if (elem_dof) { return; }
394
395 // TODO: can we call GetElementDofs only once per element?
396 Table *el_dof = new Table;
397 Table *el_fos = (mesh->Dimension() > 2) ? (new Table) : NULL;
398 Array<int> dofs;
399 Array<int> F, Fo;
400 el_dof->MakeI(mesh->GetNE());
401 if (el_fos) { el_fos->MakeI(mesh->GetNE()); }
402 for (int i = 0; i < mesh->GetNE(); i++)
403 {
404 GetElementDofs(i, dofs);
405 el_dof->AddColumnsInRow(i, dofs.Size());
406
407 if (el_fos)
408 {
409 mesh->GetElementFaces(i, F, Fo);
410 el_fos->AddColumnsInRow(i, Fo.Size());
411 }
412 }
413 el_dof->MakeJ();
414 if (el_fos) { el_fos->MakeJ(); }
415 for (int i = 0; i < mesh->GetNE(); i++)
416 {
417 GetElementDofs(i, dofs);
418 el_dof->AddConnections(i, (int *)dofs, dofs.Size());
419
420 if (el_fos)
421 {
422 mesh->GetElementFaces(i, F, Fo);
423 el_fos->AddConnections(i, (int *)Fo, Fo.Size());
424 }
425 }
426 el_dof->ShiftUpI();
427 if (el_fos) { el_fos->ShiftUpI(); }
428 elem_dof = el_dof;
429 elem_fos = el_fos;
430}
431
433{
434 if (bdr_elem_dof) { return; }
435
436 Table *bel_dof = new Table;
437 Table *bel_fos = (mesh->Dimension() == 3) ? (new Table) : NULL;
438 Array<int> dofs;
439 int F, Fo;
440 bel_dof->MakeI(mesh->GetNBE());
441 if (bel_fos) { bel_fos->MakeI(mesh->GetNBE()); }
442 for (int i = 0; i < mesh->GetNBE(); i++)
443 {
444 GetBdrElementDofs(i, dofs);
445 bel_dof->AddColumnsInRow(i, dofs.Size());
446
447 if (bel_fos)
448 {
449 bel_fos->AddAColumnInRow(i);
450 }
451 }
452 bel_dof->MakeJ();
453 if (bel_fos) { bel_fos->MakeJ(); }
454 for (int i = 0; i < mesh->GetNBE(); i++)
455 {
456 GetBdrElementDofs(i, dofs);
457 bel_dof->AddConnections(i, (int *)dofs, dofs.Size());
458
459 if (bel_fos)
460 {
461 mesh->GetBdrElementFace(i, &F, &Fo);
462 bel_fos->AddConnection(i, Fo);
463 }
464 }
465 bel_dof->ShiftUpI();
466 if (bel_fos) { bel_fos->ShiftUpI(); }
467 bdr_elem_dof = bel_dof;
468 bdr_elem_fos = bel_fos;
469}
470
472{
473 // Here, "face" == (dim-1)-dimensional mesh entity.
474
475 if (face_dof) { return; }
476
477 if (NURBSext) { BuildNURBSFaceToDofTable(); return; }
478
479 Table *fc_dof = new Table;
480 Array<int> dofs;
481 fc_dof->MakeI(mesh->GetNumFaces());
482 for (int i = 0; i < fc_dof->Size(); i++)
483 {
484 GetFaceDofs(i, dofs, 0);
485 fc_dof->AddColumnsInRow(i, dofs.Size());
486 }
487 fc_dof->MakeJ();
488 for (int i = 0; i < fc_dof->Size(); i++)
489 {
490 GetFaceDofs(i, dofs, 0);
491 fc_dof->AddConnections(i, (int *)dofs, dofs.Size());
492 }
493 fc_dof->ShiftUpI();
494 face_dof = fc_dof;
495}
496
498{
499 delete elem_dof;
500 delete elem_fos;
501 elem_dof = NULL;
502 elem_fos = NULL;
504}
505
507{
508 Array<int> dof_marker(ndofs);
509
510 dof_marker = -1;
511
512 int *J = elem_dof->GetJ(), nnz = elem_dof->Size_of_connections();
513 for (int k = 0, dof_counter = 0; k < nnz; k++)
514 {
515 const int sdof = J[k]; // signed dof
516 const int dof = (sdof < 0) ? -1-sdof : sdof;
517 int new_dof = dof_marker[dof];
518 if (new_dof < 0)
519 {
520 dof_marker[dof] = new_dof = dof_counter++;
521 }
522 J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
523 }
524}
525
527{
528 if (dof_elem_array.Size()) { return; }
529
531
534 dof_elem_array = -1;
535 for (int i = 0; i < mesh -> GetNE(); i++)
536 {
537 const int *dofs = elem_dof -> GetRow(i);
538 const int n = elem_dof -> RowSize(i);
539 for (int j = 0; j < n; j++)
540 {
541 int dof = DecodeDof(dofs[j]);
542 if (dof_elem_array[dof] < 0)
543 {
544 dof_elem_array[dof] = i;
545 dof_ldof_array[dof] = j;
546 }
547 }
548 }
549}
550
552{
553 if (dof_bdr_elem_array.Size()) { return; }
554
556
560 for (int i = 0; i < mesh -> GetNBE(); i++)
561 {
562 const int *dofs = bdr_elem_dof -> GetRow(i);
563 const int n = bdr_elem_dof -> RowSize(i);
564 for (int j = 0; j < n; j++)
565 {
566 int dof = DecodeDof(dofs[j]);
567 if (dof_bdr_elem_array[dof] < 0)
568 {
569 dof_bdr_elem_array[dof] = i;
570 dof_bdr_ldof_array[dof] = j;
571 }
572 }
573 }
574}
575
576void MarkDofs(const Array<int> &dofs, Array<int> &mark_array)
577{
578 for (auto d : dofs)
579 {
580 mark_array[d >= 0 ? d : -1 - d] = -1;
581 }
582}
583
585 Array<int> &ess_vdofs,
586 int component) const
587{
588 Array<int> dofs;
589 ess_vdofs.SetSize(GetVSize());
590 ess_vdofs = 0;
591 for (int i = 0; i < GetNBE(); i++)
592 {
593 if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
594 {
595 if (component < 0)
596 {
597 // Mark all components.
598 GetBdrElementVDofs(i, dofs);
599 }
600 else
601 {
602 GetBdrElementDofs(i, dofs);
603 for (auto &d : dofs) { d = DofToVDof(d, component); }
604 }
605 MarkDofs(dofs, ess_vdofs);
606 }
607 }
608
609 // mark possible hidden boundary edges in a non-conforming mesh, also
610 // local DOFs affected by boundary elements on other processors
611 if (Nonconforming())
612 {
613 Array<int> bdr_verts, bdr_edges, bdr_faces;
614 mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges,
615 bdr_faces);
616 for (auto v : bdr_verts)
617 {
618 if (component < 0)
619 {
620 GetVertexVDofs(v, dofs);
621 }
622 else
623 {
624 GetVertexDofs(v, dofs);
625 for (auto &d : dofs) { d = DofToVDof(d, component); }
626 }
627 MarkDofs(dofs, ess_vdofs);
628 }
629 for (auto e : bdr_edges)
630 {
631 if (component < 0)
632 {
633 GetEdgeVDofs(e, dofs);
634 }
635 else
636 {
637 GetEdgeDofs(e, dofs);
638 for (auto &d : dofs) { d = DofToVDof(d, component); }
639 }
640 MarkDofs(dofs, ess_vdofs);
641 }
642 for (auto f : bdr_faces)
643 {
644 if (component < 0)
645 {
646 GetEntityVDofs(2, f, dofs);
647 }
648 else
649 {
650 GetEntityDofs(2, f, dofs);
651 for (auto &d : dofs) { d = DofToVDof(d, component); }
652 }
653 MarkDofs(dofs, ess_vdofs);
654 }
655 }
656}
657
659 Array<int> &ess_tdof_list,
660 int component) const
661{
662 Array<int> ess_vdofs, ess_tdofs;
663 GetEssentialVDofs(bdr_attr_is_ess, ess_vdofs, component);
665 if (!R)
666 {
667 ess_tdofs.MakeRef(ess_vdofs);
668 }
669 else
670 {
671 R->BooleanMult(ess_vdofs, ess_tdofs);
672#ifdef MFEM_DEBUG
673 // Verify that in boolean arithmetic: P^T ess_dofs = R ess_dofs
674 Array<int> ess_tdofs2(ess_tdofs.Size());
675 GetConformingProlongation()->BooleanMultTranspose(ess_vdofs, ess_tdofs2);
676
677 int counter = 0;
678 std::string error_msg = "failed dof: ";
679 auto ess_tdofs_ = ess_tdofs.HostRead();
680 auto ess_tdofs2_ = ess_tdofs2.HostRead();
681 for (int i = 0; i < ess_tdofs2.Size(); ++i)
682 {
683 if (bool(ess_tdofs_[i]) != bool(ess_tdofs2_[i]))
684 {
685 error_msg += std::to_string(i) += "(R ";
686 error_msg += std::to_string(bool(ess_tdofs_[i])) += " P^T ";
687 error_msg += std::to_string(bool(ess_tdofs2_[i])) += ") ";
688 counter++;
689 }
690 }
691
692 MFEM_ASSERT(R->Height() == GetConformingProlongation()->Width(), "!");
693 MFEM_ASSERT(R->Width() == GetConformingProlongation()->Height(), "!");
694 MFEM_ASSERT(R->Width() == ess_vdofs.Size(), "!");
695 MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
696 << ' ' << error_msg);
697#endif
698 }
699 MarkerToList(ess_tdofs, ess_tdof_list);
700}
701
703 int component)
704{
705 if (mesh->bdr_attributes.Size())
706 {
707 Array<int> ess_bdr(mesh->bdr_attributes.Max());
708 ess_bdr = 1;
709 GetEssentialTrueDofs(ess_bdr, boundary_dofs, component);
710 }
711 else
712 {
713 boundary_dofs.DeleteAll();
714 }
715}
716
718 int component) const
719{
720 Array<int> dofs;
721 ext_vdofs.SetSize(GetVSize());
722 ext_vdofs = 0;
723
724 Array<int> ext_face_marker;
725 mesh->GetExteriorFaceMarker(ext_face_marker);
726 for (int i = 0; i < ext_face_marker.Size(); i++)
727 {
728 if (ext_face_marker[i])
729 {
730 if (component < 0)
731 {
732 // Mark all components.
733 GetFaceDofs(i, dofs);
734 DofsToVDofs(dofs);
735 }
736 else
737 {
738 GetFaceDofs(i, dofs);
739 for (auto &d : dofs) { d = DofToVDof(d, component); }
740 }
741 MarkDofs(dofs, ext_vdofs);
742 }
743 }
744}
745
747 int component) const
748{
749 Array<int> ext_vdofs, ext_tdofs;
750 GetExteriorVDofs(ext_vdofs, component);
752 if (!R)
753 {
754 ext_tdofs.MakeRef(ext_vdofs);
755 }
756 else
757 {
758 R->BooleanMult(ext_vdofs, ext_tdofs);
759#ifdef MFEM_DEBUG
760 // Verify that in boolean arithmetic: P^T ext_dofs = R ext_dofs
761 Array<int> ext_tdofs2(ext_tdofs.Size());
762 GetConformingProlongation()->BooleanMultTranspose(ext_vdofs, ext_tdofs2);
763
764 int counter = 0;
765 std::string error_msg = "failed dof: ";
766 auto ext_tdofs_ = ext_tdofs.HostRead();
767 auto ext_tdofs2_ = ext_tdofs2.HostRead();
768 for (int i = 0; i < ext_tdofs2.Size(); ++i)
769 {
770 if (bool(ext_tdofs_[i]) != bool(ext_tdofs2_[i]))
771 {
772 error_msg += std::to_string(i) += "(R ";
773 error_msg += std::to_string(bool(ext_tdofs_[i])) += " P^T ";
774 error_msg += std::to_string(bool(ext_tdofs2_[i])) += ") ";
775 counter++;
776 }
777 }
778
779 MFEM_ASSERT(R->Height() == GetConformingProlongation()->Width(), "!");
780 MFEM_ASSERT(R->Width() == GetConformingProlongation()->Height(), "!");
781 MFEM_ASSERT(R->Width() == ext_vdofs.Size(), "!");
782 MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
783 << ' ' << error_msg);
784#endif
785 }
786 MarkerToList(ext_tdofs, ext_tdof_list);
787}
788
789// static method
791 Array<int> &list)
792{
793 int num_marked = 0;
794 marker.HostRead(); // make sure we can read the array on host
795 for (int i = 0; i < marker.Size(); i++)
796 {
797 if (marker[i]) { num_marked++; }
798 }
799 list.SetSize(0);
800 list.HostWrite();
801 list.Reserve(num_marked);
802 for (int i = 0; i < marker.Size(); i++)
803 {
804 if (marker[i]) { list.Append(i); }
805 }
806}
807
808// static method
809void FiniteElementSpace::ListToMarker(const Array<int> &list, int marker_size,
810 Array<int> &marker, int mark_val)
811{
812 list.HostRead(); // make sure we can read the array on host
813 marker.SetSize(marker_size);
814 marker.HostWrite();
815 marker = 0;
816 for (int i = 0; i < list.Size(); i++)
817 {
818 marker[list[i]] = mark_val;
819 }
820}
821
823 Array<int> &cdofs)
824{
826 if (cP) { cP->BooleanMultTranspose(dofs, cdofs); }
827 else { dofs.Copy(cdofs); }
828}
829
831 Array<int> &dofs)
832{
834 if (cR) { cR->BooleanMultTranspose(cdofs, dofs); }
835 else { cdofs.Copy(dofs); }
836}
837
840{
841 int i, j;
842 Array<int> d_vdofs, c_vdofs;
843 SparseMatrix *R;
844
845 R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
846
847 for (i = 0; i < mesh -> GetNE(); i++)
848 {
849 this -> GetElementVDofs (i, d_vdofs);
850 cfes -> GetElementVDofs (i, c_vdofs);
851
852#ifdef MFEM_DEBUG
853 if (d_vdofs.Size() != c_vdofs.Size())
854 {
855 mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
856 }
857#endif
858
859 for (j = 0; j < d_vdofs.Size(); j++)
860 {
861 R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
862 }
863 }
864
865 R -> Finalize();
866
867 return R;
868}
869
872{
873 int i, j;
874 Array<int> d_dofs, c_dofs;
875 SparseMatrix *R;
876
877 R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
878
879 for (i = 0; i < mesh -> GetNE(); i++)
880 {
881 this -> GetElementDofs (i, d_dofs);
882 cfes -> GetElementDofs (i, c_dofs);
883
884#ifdef MFEM_DEBUG
885 if (c_dofs.Size() != 1)
886 mfem_error ("FiniteElementSpace::"
887 "D2Const_GlobalRestrictionMatrix (...)");
888#endif
889
890 for (j = 0; j < d_dofs.Size(); j++)
891 {
892 R -> Set (c_dofs[0], d_dofs[j], 1.0);
893 }
894 }
895
896 R -> Finalize();
897
898 return R;
899}
900
903{
904 SparseMatrix *R;
905 DenseMatrix loc_restr;
906 Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
907
908 int lvdim = lfes->GetVDim();
909 R = new SparseMatrix (lvdim * lfes -> GetNDofs(), lvdim * ndofs);
910
911 Geometry::Type cached_geom = Geometry::INVALID;
912 const FiniteElement *h_fe = NULL;
913 const FiniteElement *l_fe = NULL;
915
916 for (int i = 0; i < mesh -> GetNE(); i++)
917 {
918 this -> GetElementDofs (i, h_dofs);
919 lfes -> GetElementDofs (i, l_dofs);
920
921 // Assuming 'loc_restr' depends only on the Geometry::Type.
923 if (geom != cached_geom)
924 {
925 h_fe = this -> GetFE (i);
926 l_fe = lfes -> GetFE (i);
928 h_fe->Project(*l_fe, T, loc_restr);
929 cached_geom = geom;
930 }
931
932 for (int vd = 0; vd < lvdim; vd++)
933 {
934 l_dofs.Copy(l_vdofs);
935 lfes->DofsToVDofs(vd, l_vdofs);
936
937 h_dofs.Copy(h_vdofs);
938 this->DofsToVDofs(vd, h_vdofs);
939
940 R -> SetSubMatrix (l_vdofs, h_vdofs, loc_restr, 1);
941 }
942 }
943
944 R -> Finalize();
945
946 return R;
947}
948
950 SparseMatrix& deps, Array<int>& master_dofs, Array<int>& slave_dofs,
951 DenseMatrix& I, int skipfirst)
952{
953 for (int i = skipfirst; i < slave_dofs.Size(); i++)
954 {
955 const int sdof = slave_dofs[i];
956 if (!deps.RowSize(sdof)) // not processed yet
957 {
958 for (int j = 0; j < master_dofs.Size(); j++)
959 {
960 const real_t coef = I(i, j);
961 if (std::abs(coef) > 1e-12)
962 {
963 const int mdof = master_dofs[j];
964 if (mdof != sdof && mdof != (-1-sdof))
965 {
966 deps.Add(sdof, mdof, coef);
967 }
968 }
969 }
970 }
971 }
972}
973
975 SparseMatrix &deps, Array<int> &master_dofs, const FiniteElement *master_fe,
976 Array<int> &slave_dofs, int slave_face, const DenseMatrix *pm) const
977{
978 // In variable-order spaces in 3D, we need to only constrain interior face
979 // DOFs (this is done one level up), since edge dependencies can be more
980 // complex and are primarily handled by edge-edge dependencies. The one
981 // exception is edges of slave faces that lie in the interior of the master
982 // face, which are not covered by edge-edge relations. This function finds
983 // such edges and makes them constrained by the master face.
984 // See also https://github.com/mfem/mfem/pull/1423#issuecomment-633916643
985
986 Array<int> V, E, Eo; // TODO: LocalArray
987 mesh->GetFaceVertices(slave_face, V);
988 mesh->GetFaceEdges(slave_face, E, Eo);
989 MFEM_ASSERT(V.Size() == E.Size(), "");
990
991 DenseMatrix I;
993 edge_T.SetFE(&SegmentFE);
994
995 // constrain each edge of the slave face
996 for (int i = 0; i < E.Size(); i++)
997 {
998 int a = i, b = (i+1) % V.Size();
999 if (V[a] > V[b]) { std::swap(a, b); }
1000
1001 DenseMatrix &edge_pm = edge_T.GetPointMat();
1002 edge_pm.SetSize(2, 2);
1003
1004 // copy two points from the face point matrix
1005 real_t mid[2];
1006 for (int j = 0; j < 2; j++)
1007 {
1008 edge_pm(j, 0) = (*pm)(j, a);
1009 edge_pm(j, 1) = (*pm)(j, b);
1010 mid[j] = 0.5*((*pm)(j, a) + (*pm)(j, b));
1011 }
1012
1013 // check that the edge does not coincide with the master face's edge
1014 const real_t eps = 1e-14;
1015 if (mid[0] > eps && mid[0] < 1-eps &&
1016 mid[1] > eps && mid[1] < 1-eps)
1017 {
1018 int order = GetEdgeDofs(E[i], slave_dofs, 0);
1019
1020 const auto *edge_fe = fec->GetFE(Geometry::SEGMENT, order);
1021 edge_fe->GetTransferMatrix(*master_fe, edge_T, I);
1022
1023 AddDependencies(deps, master_dofs, slave_dofs, I, 0);
1024 }
1025 }
1026}
1027
1029 const SparseMatrix& deps)
1030{
1031 const int* dep = deps.GetRowColumns(dof);
1032 int ndep = deps.RowSize(dof);
1033
1034 // are all constraining DOFs finalized?
1035 for (int i = 0; i < ndep; i++)
1036 {
1037 if (!finalized[dep[i]]) { return false; }
1038 }
1039 return true;
1040}
1041
1043 Geometry::Type master_geom,
1044 int variant) const
1045{
1046 // In NC meshes with prisms/tets, a special constraint occurs where a
1047 // prism/tet edge is slave to another element's face (see illustration
1048 // here: https://github.com/mfem/mfem/pull/713#issuecomment-495786362)
1049 // Rather than introduce a new edge-face constraint type, we handle such
1050 // cases as degenerate face-face constraints, where the point-matrix
1051 // rectangle has zero height. This method returns DOFs for the first edge
1052 // of the rectangle, duplicated in the orthogonal direction, to resemble
1053 // DOFs for a quadrilateral face. The extra DOFs are ignored by
1054 // FiniteElementSpace::AddDependencies.
1055
1056 Array<int> edof;
1057 int order = GetEdgeDofs(-1 - index, edof, variant);
1058
1061 int nn = 2*nv + ne;
1062
1063 dofs.SetSize(nn*nn);
1064 if (!dofs.Size()) { return 0; }
1065
1066 dofs = edof[0];
1067
1068 // copy first two vertex DOFs
1069 for (int i = 0; i < nv; i++)
1070 {
1071 dofs[i] = edof[i];
1072 dofs[nv+i] = edof[nv+i];
1073 }
1074 // copy first edge DOFs
1075 int face_vert = Geometry::NumVerts[master_geom];
1076 for (int i = 0; i < ne; i++)
1077 {
1078 dofs[face_vert*nv + i] = edof[2*nv + i];
1079 }
1080
1081 return order;
1082}
1083
1085{
1086 // return the number of vertex and edge DOFs that precede inner DOFs
1087 const int nv = fec->GetNumDof(Geometry::POINT, order);
1088 const int ne = fec->GetNumDof(Geometry::SEGMENT, order);
1089
1090 return Geometry::NumVerts[geom] * (geom == Geometry::SEGMENT ? nv : (nv + ne));
1091}
1092
1094 Geometry::Type master_geom,
1095 int variant) const
1096{
1097 switch (entity)
1098 {
1099 case 0:
1100 GetVertexDofs(index, dofs);
1101 return 0;
1102
1103 case 1:
1104 return GetEdgeDofs(index, dofs, variant);
1105
1106 default:
1107 if (index >= 0)
1108 {
1109 return GetFaceDofs(index, dofs, variant);
1110 }
1111 else
1112 {
1113 return GetDegenerateFaceDofs(index, dofs, master_geom, variant);
1114 }
1115 }
1116}
1117
1119 Geometry::Type master_geom,
1120 int variant) const
1121{
1122 const int n = GetEntityDofs(entity, index, dofs, master_geom, variant);
1123 DofsToVDofs(dofs);
1124 return n;
1125}
1126
1127// Variable-order spaces: enforce minimum rule on conforming edges/faces
1129{
1130 if (!IsVariableOrder()) { return; }
1131
1132 Array<int> master_dofs, slave_dofs;
1133
1135 DenseMatrix I;
1136
1137 for (int entity = 1; entity < mesh->Dimension(); entity++)
1138 {
1139 const Table &ent_dofs = (entity == 1) ? var_edge_dofs : var_face_dofs;
1140 const int num_ent = (entity == 1) ? mesh->GetNEdges() : mesh->GetNFaces();
1141 MFEM_ASSERT(ent_dofs.Size() >= num_ent+1, "");
1142
1143 // add constraints within edges/faces holding multiple DOF sets
1145 for (int i = 0; i < num_ent; i++)
1146 {
1147 if (ent_dofs.RowSize(i) <= 1) { continue; }
1148
1149 Geometry::Type geom =
1150 (entity == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
1151
1152 if (geom != last_geom)
1153 {
1155 last_geom = geom;
1156 }
1157
1158 // get lowest order variant DOFs and FE
1159 const int p = GetEntityDofs(entity, i, master_dofs, geom, 0);
1160 const auto *master_fe = fec->GetFE(geom, p);
1161 if (!master_fe) { break; }
1162
1163 // constrain all higher order DOFs: interpolate lowest order function
1164 for (int variant = 1; ; variant++)
1165 {
1166 const int q = GetEntityDofs(entity, i, slave_dofs, geom, variant);
1167 if (q < 0) { break; }
1168
1169 const auto *slave_fe = fec->GetFE(geom, q);
1170 slave_fe->GetTransferMatrix(*master_fe, T, I);
1171
1172 AddDependencies(deps, master_dofs, slave_dofs, I);
1173 }
1174 }
1175 }
1176}
1177
1179{
1180#ifdef MFEM_USE_MPI
1181 MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
1182 "This method should not be used with a ParFiniteElementSpace!");
1183#endif
1184
1185 if (cP_is_set) { return; }
1186 cP_is_set = true;
1187
1188 if (FEColl()->GetContType() == FiniteElementCollection::DISCONTINUOUS)
1189 {
1190 cP.reset();
1191 cR.reset();
1192 cR_hp.reset();
1193 R_transpose.reset();
1194 return;
1195 }
1196
1197 Array<int> master_dofs, slave_dofs, highest_dofs;
1198
1200 DenseMatrix I;
1201
1202 // For each slave DOF, the dependency matrix will contain a row that
1203 // expresses the slave DOF as a linear combination of its immediate master
1204 // DOFs. Rows of independent DOFs will remain empty.
1205 SparseMatrix deps(ndofs);
1206
1207 // Inverse dependencies for the cR_hp matrix in variable-order spaces:
1208 // For each master edge/face with more DOF sets, the inverse dependency
1209 // matrix contains a row that expresses the master true DOF (lowest order)
1210 // as a linear combination of the highest order set of DOFs.
1211 SparseMatrix inv_deps(ndofs);
1212
1214
1215 // Collect local face/edge dependencies, starting with faces
1216 for (int entity = 2; entity >= 1; entity--)
1217 {
1218 const NCMesh::NCList &list = mesh->ncmesh->GetNCList(entity);
1219 if (!list.masters.Size()) { continue; }
1220
1221 // loop through all master edges/faces, constrain their slave edges/faces
1222 for (const NCMesh::Master &master : list.masters)
1223 {
1224 Geometry::Type master_geom = master.Geom();
1225
1226 const int p = GetEntityDofs(entity, master.index, master_dofs,
1227 master_geom);
1228 if (!master_dofs.Size()) { continue; }
1229
1230 const FiniteElement *master_fe = fec->GetFE(master_geom, p);
1231 if (!master_fe) { continue; }
1232
1233 switch (master_geom)
1234 {
1235 case Geometry::SQUARE: T.SetFE(&QuadrilateralFE); break;
1236 case Geometry::TRIANGLE: T.SetFE(&TriangleFE); break;
1237 case Geometry::SEGMENT: T.SetFE(&SegmentFE); break;
1238 default: MFEM_ABORT("unsupported geometry");
1239 }
1240
1241 for (int si = master.slaves_begin; si < master.slaves_end; si++)
1242 {
1243 const NCMesh::Slave &slave = list.slaves[si];
1244
1245 int q = GetEntityDofs(entity, slave.index, slave_dofs, master_geom);
1246 if (!slave_dofs.Size()) { break; }
1247
1248 const FiniteElement *slave_fe = fec->GetFE(slave.Geom(), q);
1249 list.OrientedPointMatrix(slave, T.GetPointMat());
1250 slave_fe->GetTransferMatrix(*master_fe, T, I);
1251
1252 // variable-order spaces: face edges need to be handled separately
1253 int skipfirst = 0;
1254 if (IsVariableOrder() && entity == 2 && slave.index >= 0)
1255 {
1256 skipfirst = GetNumBorderDofs(master_geom, q);
1257 }
1258
1259 // make each slave DOF dependent on all master DOFs
1260 AddDependencies(deps, master_dofs, slave_dofs, I, skipfirst);
1261
1262 if (skipfirst)
1263 {
1264 // constrain internal edge DOFs if they were skipped
1265 const auto *pm = list.point_matrices[master_geom][slave.matrix];
1266 AddEdgeFaceDependencies(deps, master_dofs, master_fe,
1267 slave_dofs, slave.index, pm);
1268 }
1269 }
1270
1271 // Add inverse dependencies for the cR_hp matrix; if a master has
1272 // more DOF sets, the lowest order set interpolates the highest one.
1273 if (IsVariableOrder())
1274 {
1275 int nvar = GetNVariants(entity, master.index);
1276 if (nvar > 1)
1277 {
1278 const int q = GetEntityDofs(entity, master.index, highest_dofs,
1279 master_geom, nvar-1);
1280 const auto *highest_fe = fec->GetFE(master_geom, q);
1281
1282 T.SetIdentityTransformation(master_geom);
1283 master_fe->GetTransferMatrix(*highest_fe, T, I);
1284
1285 // add dependencies only for the inner dofs
1286 const int skip = GetNumBorderDofs(master_geom, p);
1287 AddDependencies(inv_deps, highest_dofs, master_dofs, I, skip);
1288 }
1289 }
1290 }
1291 }
1292
1293 deps.Finalize();
1294 inv_deps.Finalize();
1295
1296 // DOFs that stayed independent are true DOFs
1297 int n_true_dofs = 0;
1298 for (int i = 0; i < ndofs; i++)
1299 {
1300 if (!deps.RowSize(i)) { n_true_dofs++; }
1301 }
1302
1303 // if all dofs are true dofs leave cP and cR NULL
1304 if (n_true_dofs == ndofs)
1305 {
1306 cP.reset();
1307 cR.reset();
1308 cR_hp.reset();
1309 R_transpose.reset();
1310 return;
1311 }
1312
1313 // create the conforming prolongation matrix cP
1314 cP.reset(new SparseMatrix(ndofs, n_true_dofs));
1315
1316 // create the conforming restriction matrix cR
1317 int *cR_J;
1318 {
1319 int *cR_I = Memory<int>(n_true_dofs+1);
1320 real_t *cR_A = Memory<real_t>(n_true_dofs);
1321 cR_J = Memory<int>(n_true_dofs);
1322 for (int i = 0; i < n_true_dofs; i++)
1323 {
1324 cR_I[i] = i;
1325 cR_A[i] = 1.0;
1326 }
1327 cR_I[n_true_dofs] = n_true_dofs;
1328 cR.reset(new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, ndofs));
1329 }
1330
1331 // In variable-order spaces, create the restriction matrix cR_hp, which is
1332 // similar to cR but has interpolation of the master edge/face DOFs of
1333 // maximum order per edge/face, since the maximum order is on an adjacent
1334 // element (e.g. where projection would be computed).
1335 if (IsVariableOrder())
1336 {
1337 cR_hp.reset(new SparseMatrix(n_true_dofs, ndofs));
1338 }
1339 else
1340 {
1341 cR_hp.reset();
1342 }
1343
1344 Array<bool> finalized(ndofs);
1345 finalized = false;
1346
1347 Array<int> cols;
1348 Vector srow;
1349
1350 // Put identity in the prolongation matrix for true DOFs, and set cR_hp
1351 for (int i = 0, true_dof = 0; i < ndofs; i++)
1352 {
1353 if (!deps.RowSize(i)) // true dof
1354 {
1355 cP->Add(i, true_dof, 1.0);
1356 cR_J[true_dof] = i;
1357 finalized[i] = true;
1358
1359 if (cR_hp)
1360 {
1361 if (inv_deps.RowSize(i))
1362 {
1363 inv_deps.GetRow(i, cols, srow);
1364 cR_hp->AddRow(true_dof, cols, srow);
1365 }
1366 else
1367 {
1368 cR_hp->Add(true_dof, i, 1.0);
1369 }
1370 }
1371
1372 true_dof++;
1373 }
1374 }
1375
1376 // Now calculate cP rows of slave DOFs as combinations of cP rows of their
1377 // master DOFs. It is possible that some slave DOFs depend on DOFs that are
1378 // themselves slaves. Here we resolve such indirect constraints by first
1379 // calculating rows of the cP matrix for DOFs whose master DOF cP rows are
1380 // already known (in the first iteration these are the true DOFs). In the
1381 // second iteration, slaves of slaves can be 'finalized' (given a row in the
1382 // cP matrix), in the third iteration slaves of slaves of slaves, etc.
1383 bool finished;
1384 int n_finalized = n_true_dofs;
1385 do
1386 {
1387 finished = true;
1388 for (int dof = 0; dof < ndofs; dof++)
1389 {
1390 if (!finalized[dof] && DofFinalizable(dof, finalized, deps))
1391 {
1392 const int* dep_col = deps.GetRowColumns(dof);
1393 const real_t* dep_coef = deps.GetRowEntries(dof);
1394 int n_dep = deps.RowSize(dof);
1395
1396 for (int j = 0; j < n_dep; j++)
1397 {
1398 cP->GetRow(dep_col[j], cols, srow);
1399 srow *= dep_coef[j];
1400 cP->AddRow(dof, cols, srow);
1401 }
1402
1403 finalized[dof] = true;
1404 n_finalized++;
1405 finished = false;
1406 }
1407 }
1408 }
1409 while (!finished);
1410
1411 // If everything is consistent (mesh, face orientations, etc.), we should
1412 // be able to finalize all slave DOFs, otherwise it's a serious error.
1413 MFEM_VERIFY(n_finalized == ndofs,
1414 "Error creating cP matrix: n_finalized = "
1415 << n_finalized << ", ndofs = " << ndofs);
1416
1417 cP->Finalize();
1418 if (cR_hp) { cR_hp->Finalize(); }
1419
1420 if (vdim > 1)
1421 {
1424 if (cR_hp) { MakeVDimMatrix(*cR_hp); }
1425 }
1426}
1427
1429{
1430 if (vdim == 1) { return; }
1431
1432 int height = mat.Height();
1433 int width = mat.Width();
1434
1435 SparseMatrix *vmat = new SparseMatrix(vdim*height, vdim*width);
1436
1437 Array<int> dofs, vdofs;
1438 Vector srow;
1439 for (int i = 0; i < height; i++)
1440 {
1441 mat.GetRow(i, dofs, srow);
1442 for (int vd = 0; vd < vdim; vd++)
1443 {
1444 dofs.Copy(vdofs);
1445 DofsToVDofs(vd, vdofs, width);
1446 vmat->SetRow(DofToVDof(i, vd, height), vdofs, srow);
1447 }
1448 }
1449 vmat->Finalize();
1450
1451 mat.Swap(*vmat);
1452 delete vmat;
1453}
1454
1455
1457{
1458 if (Conforming()) { return NULL; }
1460 return cP.get();
1461}
1462
1464{
1465 if (Conforming()) { return NULL; }
1467 if (cR && !R_transpose) { R_transpose.reset(new TransposeOperator(*cR)); }
1468 return cR.get();
1469}
1470
1472{
1473 if (Conforming()) { return NULL; }
1475 return IsVariableOrder() ? cR_hp.get() : cR.get();
1476}
1477
1479{
1480 GetRestrictionOperator(); // Ensure that R_transpose is built
1481 return R_transpose.get();
1482}
1483
1485{
1487 return P ? (P->Width() / vdim) : ndofs;
1488}
1489
1491{
1492 const FiniteElement *fe = GetTypicalFE();
1494 {
1495 return GetVDim();
1496 }
1497 return GetVDim()*std::max(GetMesh()->SpaceDimension(), fe->GetRangeDim());
1498}
1499
1501{
1502 const FiniteElement *fe = GetTypicalFE();
1504 {
1505 return 2 * GetMesh()->SpaceDimension() - 3;
1506 }
1507 return GetVDim()*fe->GetCurlDim();
1508}
1509
1511 ElementDofOrdering e_ordering) const
1512{
1513 // Check if we have a discontinuous space using the FE collection:
1514 if (IsDGSpace())
1515 {
1516 // TODO: when VDIM is 1, we can return IdentityOperator.
1517 if (L2E_nat.Ptr() == NULL)
1518 {
1519 // The input L-vector layout is:
1520 // * ND x NE x VDIM, for Ordering::byNODES, or
1521 // * VDIM x ND x NE, for Ordering::byVDIM.
1522 // The output E-vector layout is: ND x VDIM x NE.
1524 }
1526 }
1527 if (e_ordering == ElementDofOrdering::LEXICOGRAPHIC)
1528 {
1529 if (L2E_lex.Ptr() == NULL)
1530 {
1531 L2E_lex.Reset(new ElementRestriction(*this, e_ordering));
1532 }
1534 }
1535 // e_ordering == ElementDofOrdering::NATIVE
1536 if (L2E_nat.Ptr() == NULL)
1537 {
1538 L2E_nat.Reset(new ElementRestriction(*this, e_ordering));
1539 }
1541}
1542
1544 ElementDofOrdering f_ordering, FaceType type, L2FaceValues mul) const
1545{
1546 const bool is_dg_space = IsDGSpace();
1547 const L2FaceValues m = (is_dg_space && mul==L2FaceValues::DoubleValued) ?
1549 key_face key = std::make_tuple(is_dg_space, f_ordering, type, m);
1550 auto itr = L2F.find(key);
1551 if (itr != L2F.end())
1552 {
1553 return itr->second;
1554 }
1555 else
1556 {
1557 FaceRestriction *res;
1558 if (is_dg_space)
1559 {
1560 if (Conforming())
1561 {
1562 res = new L2FaceRestriction(*this, f_ordering, type, m);
1563 }
1564 else
1565 {
1566 res = new NCL2FaceRestriction(*this, f_ordering, type, m);
1567 }
1568 }
1569 else if (dynamic_cast<const DG_Interface_FECollection*>(fec))
1570 {
1571 res = new L2InterfaceFaceRestriction(*this, f_ordering, type);
1572 }
1573 else
1574 {
1575 res = new ConformingFaceRestriction(*this, f_ordering, type);
1576 }
1577 L2F[key] = res;
1578 return res;
1579 }
1580}
1581
1583 const IntegrationRule &ir) const
1584{
1585 for (int i = 0; i < E2Q_array.Size(); i++)
1586 {
1587 const QuadratureInterpolator *qi = E2Q_array[i];
1588 if (qi->IntRule == &ir) { return qi; }
1589 }
1590
1592 E2Q_array.Append(qi);
1593 return qi;
1594}
1595
1597 const QuadratureSpace &qs) const
1598{
1599 for (int i = 0; i < E2Q_array.Size(); i++)
1600 {
1601 const QuadratureInterpolator *qi = E2Q_array[i];
1602 if (qi->qspace == &qs) { return qi; }
1603 }
1604
1606 E2Q_array.Append(qi);
1607 return qi;
1608}
1609
1612 const IntegrationRule &ir, FaceType type) const
1613{
1614 if (type==FaceType::Interior)
1615 {
1616 for (int i = 0; i < E2IFQ_array.Size(); i++)
1617 {
1619 if (qi->IntRule == &ir) { return qi; }
1620 }
1621
1623 type);
1624 E2IFQ_array.Append(qi);
1625 return qi;
1626 }
1627 else //Boundary
1628 {
1629 for (int i = 0; i < E2BFQ_array.Size(); i++)
1630 {
1632 if (qi->IntRule == &ir) { return qi; }
1633 }
1634
1636 type);
1637 E2BFQ_array.Append(qi);
1638 return qi;
1639 }
1640}
1641
1643 const int coarse_ndofs, const Table &coarse_elem_dof,
1644 const Table *coarse_elem_fos, const DenseTensor localP[]) const
1645{
1646 /// TODO: Implement DofTransformation support
1647
1648 MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
1649
1650 Array<int> dofs, coarse_dofs, coarse_vdofs;
1651 Vector row;
1652
1653 Mesh::GeometryList elem_geoms(*mesh);
1654
1655 SparseMatrix *P;
1656 if (elem_geoms.Size() == 1)
1657 {
1658 const int coarse_ldof = localP[elem_geoms[0]].SizeJ();
1659 P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim, coarse_ldof);
1660 }
1661 else
1662 {
1663 P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim);
1664 }
1665
1666 Array<int> mark(P->Height());
1667 mark = 0;
1668
1670
1671 for (int k = 0; k < mesh->GetNE(); k++)
1672 {
1673 const Embedding &emb = rtrans.embeddings[k];
1675 const DenseMatrix &lP = localP[geom](emb.matrix);
1676 const int fine_ldof = localP[geom].SizeI();
1677
1678 elem_dof->GetRow(k, dofs);
1679 coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
1680
1681 for (int vd = 0; vd < vdim; vd++)
1682 {
1683 coarse_dofs.Copy(coarse_vdofs);
1684 DofsToVDofs(vd, coarse_vdofs, coarse_ndofs);
1685
1686 for (int i = 0; i < fine_ldof; i++)
1687 {
1688 int r = DofToVDof(dofs[i], vd);
1689 int m = (r >= 0) ? r : (-1 - r);
1690
1691 if (!mark[m])
1692 {
1693 lP.GetRow(i, row);
1694 P->SetRow(r, coarse_vdofs, row);
1695 mark[m] = 1;
1696 }
1697 }
1698 }
1699 }
1700
1701 MFEM_ASSERT(mark.Sum() == P->Height(), "Not all rows of P set.");
1702 if (elem_geoms.Size() != 1) { P->Finalize(); }
1703 return P;
1704}
1705
1707 const int coarse_ndofs, const Table &coarse_elem_dof) const
1708{
1709 MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
1710
1711 Array<int> dofs, coarse_dofs, coarse_vdofs;
1712 Vector row;
1713
1714 Mesh::GeometryList elem_geoms(*mesh);
1715
1716 SparseMatrix *P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim);
1717
1718 Array<int> mark(P->Height());
1719 mark = 0;
1720
1722 DenseMatrix lP;
1724 for (int k = 0; k < mesh->GetNE(); k++)
1725 {
1726 const Embedding &emb = rtrans.embeddings[k];
1728
1729 const FiniteElement *fe = GetFE(k);
1730 isotr.SetIdentityTransformation(geom);
1731 const int ldof = fe->GetDof();
1732 lP.SetSize(ldof, ldof);
1733 const DenseTensor &pmats = rtrans.point_matrices[geom];
1734 isotr.SetPointMat(pmats(emb.matrix));
1735 fe->GetLocalInterpolation(isotr, lP);
1736
1737 const int fine_ldof = lP.Height();
1738
1739 elem_dof->GetRow(k, dofs);
1740 coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
1741
1742 for (int vd = 0; vd < vdim; vd++)
1743 {
1744 coarse_dofs.Copy(coarse_vdofs);
1745 DofsToVDofs(vd, coarse_vdofs, coarse_ndofs);
1746
1747 for (int i = 0; i < fine_ldof; i++)
1748 {
1749 const int r = DofToVDof(dofs[i], vd);
1750 int m = (r >= 0) ? r : (-1 - r);
1751
1752 if (!mark[m])
1753 {
1754 lP.GetRow(i, row);
1755 P->SetRow(r, coarse_vdofs, row);
1756 mark[m] = 1;
1757 }
1758 }
1759 }
1760 }
1761
1762 MFEM_VERIFY(mark.Sum() == P->Height(), "Not all rows of P set.");
1763 P->Finalize();
1764 return P;
1765}
1766
1768 Geometry::Type geom, DenseTensor &localP) const
1769{
1770 const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
1771
1773 const DenseTensor &pmats = rtrans.point_matrices[geom];
1774
1775 int nmat = pmats.SizeK();
1776 int ldof = fe->GetDof();
1777
1779 isotr.SetIdentityTransformation(geom);
1780
1781 // calculate local interpolation matrices for all refinement types
1782 localP.SetSize(ldof, ldof, nmat);
1783 for (int i = 0; i < nmat; i++)
1784 {
1785 isotr.SetPointMat(pmats(i));
1786 fe->GetLocalInterpolation(isotr, localP(i));
1787 }
1788}
1789
1791 const Table* old_elem_dof,
1792 const Table* old_elem_fos)
1793{
1794 MFEM_VERIFY(GetNE() >= old_elem_dof->Size(),
1795 "Previous mesh is not coarser.");
1796
1797 Mesh::GeometryList elem_geoms(*mesh);
1798 if (!IsVariableOrder())
1799 {
1801 for (int i = 0; i < elem_geoms.Size(); i++)
1802 {
1803 GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1804 }
1805 return RefinementMatrix_main(old_ndofs, *old_elem_dof, old_elem_fos,
1806 localP);
1807 }
1808 else
1809 {
1810 return VariableOrderRefinementMatrix(old_ndofs, *old_elem_dof);
1811 }
1812}
1813
1815 const FiniteElementSpace* fespace, Table* old_elem_dof, Table* old_elem_fos,
1816 int old_ndofs)
1817 : fespace(fespace),
1818 old_elem_dof(old_elem_dof),
1819 old_elem_fos(old_elem_fos)
1820{
1821 MFEM_VERIFY(fespace->GetNE() >= old_elem_dof->Size(),
1822 "Previous mesh is not coarser.");
1823
1824 width = old_ndofs * fespace->GetVDim();
1825 height = fespace->GetVSize();
1826
1827 Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1828
1829 if (!fespace->IsVariableOrder())
1830 {
1831 for (int i = 0; i < elem_geoms.Size(); i++)
1832 {
1833 fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1834 }
1835 }
1836
1837 ConstructDoFTransArray();
1838}
1839
1841 const FiniteElementSpace *fespace, const FiniteElementSpace *coarse_fes)
1842 : Operator(fespace->GetVSize(), coarse_fes->GetVSize()),
1843 fespace(fespace), old_elem_dof(NULL), old_elem_fos(NULL)
1844{
1845 Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1846
1847 if (!fespace->IsVariableOrder())
1848 {
1849 for (int i = 0; i < elem_geoms.Size(); i++)
1850 {
1851 fespace->GetLocalRefinementMatrices(*coarse_fes, elem_geoms[i],
1852 localP[elem_geoms[i]]);
1853 }
1854 }
1855
1856 // Make a copy of the coarse elem_dof Table.
1857 old_elem_dof = new Table(coarse_fes->GetElementToDofTable());
1858
1859 // Make a copy of the coarse elem_fos Table if it exists.
1860 if (coarse_fes->GetElementToFaceOrientationTable())
1861 {
1862 old_elem_fos = new Table(*coarse_fes->GetElementToFaceOrientationTable());
1863 }
1864
1865 ConstructDoFTransArray();
1866}
1867
1869{
1870 delete old_elem_dof;
1871 delete old_elem_fos;
1872 for (int i=0; i<old_DoFTransArray.Size(); i++)
1873 {
1874 delete old_DoFTransArray[i];
1875 }
1876}
1877
1878void FiniteElementSpace::RefinementOperator::ConstructDoFTransArray()
1879{
1880 old_DoFTransArray.SetSize(Geometry::NUM_GEOMETRIES);
1881 for (int i=0; i<old_DoFTransArray.Size(); i++)
1882 {
1883 old_DoFTransArray[i] = NULL;
1884 }
1885
1886 const FiniteElementCollection *fec_ref = fespace->FEColl();
1887 if (dynamic_cast<const ND_FECollection*>(fec_ref))
1888 {
1889 const FiniteElement *nd_tri =
1891 if (nd_tri)
1892 {
1893 old_DoFTransArray[Geometry::TRIANGLE] =
1894 new ND_TriDofTransformation(nd_tri->GetOrder());
1895 }
1896
1897 const FiniteElement *nd_tet =
1899 if (nd_tet)
1900 {
1901 old_DoFTransArray[Geometry::TETRAHEDRON] =
1902 new ND_TetDofTransformation(nd_tet->GetOrder());
1903 }
1904
1905 const FiniteElement *nd_pri =
1907 if (nd_pri)
1908 {
1909 old_DoFTransArray[Geometry::PRISM] =
1910 new ND_WedgeDofTransformation(nd_pri->GetOrder());
1911 }
1912
1913 const FiniteElement *nd_pyr =
1915 if (nd_pyr)
1916 {
1917 old_DoFTransArray[Geometry::PYRAMID] =
1918 new ND_PyramidDofTransformation(nd_pyr->GetOrder());
1919 }
1920 }
1921}
1922
1924 Vector &y) const
1925{
1926 Mesh* mesh_ref = fespace->GetMesh();
1927 const CoarseFineTransformations &trans_ref =
1928 mesh_ref->GetRefinementTransforms();
1929
1930 Array<int> dofs, vdofs, old_dofs, old_vdofs, old_Fo;
1931
1932 int rvdim = fespace->GetVDim();
1933 int old_ndofs = width / rvdim;
1934
1935 Vector subY, subX;
1936
1937 DenseMatrix eP;
1939
1940 for (int k = 0; k < mesh_ref->GetNE(); k++)
1941 {
1942 const Embedding &emb = trans_ref.embeddings[k];
1943 const Geometry::Type geom = mesh_ref->GetElementBaseGeometry(k);
1944 if (fespace->IsVariableOrder())
1945 {
1946 const FiniteElement *fe = fespace->GetFE(k);
1947 isotr.SetIdentityTransformation(geom);
1948 const int ldof = fe->GetDof();
1949 eP.SetSize(ldof, ldof);
1950 const DenseTensor &pmats = trans_ref.point_matrices[geom];
1951 isotr.SetPointMat(pmats(emb.matrix));
1952 fe->GetLocalInterpolation(isotr, eP);
1953 }
1954 const DenseMatrix &lP = (fespace->IsVariableOrder()) ? eP : localP[geom](
1955 emb.matrix);
1956
1957 subY.SetSize(lP.Height());
1958
1959 DofTransformation *doftrans = fespace->GetElementDofs(k, dofs);
1960 old_elem_dof->GetRow(emb.parent, old_dofs);
1961
1962 if (!doftrans)
1963 {
1964 for (int vd = 0; vd < rvdim; vd++)
1965 {
1966 dofs.Copy(vdofs);
1967 fespace->DofsToVDofs(vd, vdofs);
1968 old_dofs.Copy(old_vdofs);
1969 fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
1970
1971 x.GetSubVector(old_vdofs, subX);
1972 lP.Mult(subX, subY);
1973 y.SetSubVector(vdofs, subY);
1974 }
1975 }
1976 else
1977 {
1978 old_elem_fos->GetRow(emb.parent, old_Fo);
1979 old_DoFTrans.SetDofTransformation(*old_DoFTransArray[geom]);
1980 old_DoFTrans.SetFaceOrientations(old_Fo);
1981
1982 doftrans->SetVDim();
1983 for (int vd = 0; vd < rvdim; vd++)
1984 {
1985 dofs.Copy(vdofs);
1986 fespace->DofsToVDofs(vd, vdofs);
1987 old_dofs.Copy(old_vdofs);
1988 fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
1989
1990 x.GetSubVector(old_vdofs, subX);
1991 old_DoFTrans.InvTransformPrimal(subX);
1992 lP.Mult(subX, subY);
1993 doftrans->TransformPrimal(subY);
1994 y.SetSubVector(vdofs, subY);
1995 }
1996 doftrans->SetVDim(rvdim, fespace->GetOrdering());
1997 }
1998 }
1999}
2000
2002 Vector &y) const
2003{
2004 y = 0.0;
2005
2006 Mesh* mesh_ref = fespace->GetMesh();
2007 const CoarseFineTransformations &trans_ref =
2008 mesh_ref->GetRefinementTransforms();
2009
2010 Array<char> processed(fespace->GetVSize());
2011 processed = 0;
2012
2013 Array<int> f_dofs, c_dofs, f_vdofs, c_vdofs, old_Fo;
2014
2015 int rvdim = fespace->GetVDim();
2016 int old_ndofs = width / rvdim;
2017
2018 Vector subY, subX, subYt;
2019
2020 DenseMatrix eP;
2022 const FiniteElement *fe = nullptr;
2023
2024 for (int k = 0; k < mesh_ref->GetNE(); k++)
2025 {
2026 const Embedding &emb = trans_ref.embeddings[k];
2027 const Geometry::Type geom = mesh_ref->GetElementBaseGeometry(k);
2028
2029 if (fespace->IsVariableOrder())
2030 {
2031 fe = fespace->GetFE(k);
2032 isotr.SetIdentityTransformation(geom);
2033 const int ldof = fe->GetDof();
2034 eP.SetSize(ldof);
2035 const DenseTensor &pmats = trans_ref.point_matrices[geom];
2036 isotr.SetPointMat(pmats(emb.matrix));
2037 fe->GetLocalInterpolation(isotr, eP);
2038 }
2039
2040 const DenseMatrix &lP = (fespace->IsVariableOrder()) ? eP : localP[geom](
2041 emb.matrix);
2042
2043 DofTransformation *doftrans = fespace->GetElementDofs(k, f_dofs);
2044 old_elem_dof->GetRow(emb.parent, c_dofs);
2045
2046 if (!doftrans)
2047 {
2048 subY.SetSize(lP.Width());
2049
2050 for (int vd = 0; vd < rvdim; vd++)
2051 {
2052 f_dofs.Copy(f_vdofs);
2053 fespace->DofsToVDofs(vd, f_vdofs);
2054 c_dofs.Copy(c_vdofs);
2055 fespace->DofsToVDofs(vd, c_vdofs, old_ndofs);
2056
2057 x.GetSubVector(f_vdofs, subX);
2058 for (int p = 0; p < f_dofs.Size(); ++p)
2059 {
2060 if (processed[DecodeDof(f_dofs[p])])
2061 {
2062 subX[p] = 0.0;
2063 }
2064 }
2065 lP.MultTranspose(subX, subY);
2066 y.AddElementVector(c_vdofs, subY);
2067 }
2068 }
2069 else
2070 {
2071 subYt.SetSize(lP.Width());
2072
2073 old_elem_fos->GetRow(emb.parent, old_Fo);
2074 old_DoFTrans.SetDofTransformation(*old_DoFTransArray[geom]);
2075 old_DoFTrans.SetFaceOrientations(old_Fo);
2076
2077 doftrans->SetVDim();
2078 for (int vd = 0; vd < rvdim; vd++)
2079 {
2080 f_dofs.Copy(f_vdofs);
2081 fespace->DofsToVDofs(vd, f_vdofs);
2082 c_dofs.Copy(c_vdofs);
2083 fespace->DofsToVDofs(vd, c_vdofs, old_ndofs);
2084
2085 x.GetSubVector(f_vdofs, subX);
2086 doftrans->InvTransformDual(subX);
2087 for (int p = 0; p < f_dofs.Size(); ++p)
2088 {
2089 if (processed[DecodeDof(f_dofs[p])])
2090 {
2091 subX[p] = 0.0;
2092 }
2093 }
2094 lP.MultTranspose(subX, subYt);
2095 old_DoFTrans.TransformDual(subYt);
2096 y.AddElementVector(c_vdofs, subYt);
2097 }
2098 doftrans->SetVDim(rvdim, fespace->GetOrdering());
2099 }
2100
2101 for (int p = 0; p < f_dofs.Size(); ++p)
2102 {
2103 processed[DecodeDof(f_dofs[p])] = 1;
2104 }
2105 }
2106}
2107
2108namespace internal
2109{
2110
2111// Used in GetCoarseToFineMap() below.
2112struct RefType
2113{
2114 Geometry::Type geom;
2115 int num_children;
2116 const Pair<int,int> *children;
2117
2118 RefType(Geometry::Type g, int n, const Pair<int,int> *c)
2119 : geom(g), num_children(n), children(c) { }
2120
2121 bool operator<(const RefType &other) const
2122 {
2123 if (geom < other.geom) { return true; }
2124 if (geom > other.geom) { return false; }
2125 if (num_children < other.num_children) { return true; }
2126 if (num_children > other.num_children) { return false; }
2127 for (int i = 0; i < num_children; i++)
2128 {
2129 if (children[i].one < other.children[i].one) { return true; }
2130 if (children[i].one > other.children[i].one) { return false; }
2131 }
2132 return false; // everything is equal
2133 }
2134};
2135
2136void GetCoarseToFineMap(const CoarseFineTransformations &cft,
2137 const mfem::Mesh &fine_mesh,
2138 Table &coarse_to_fine,
2139 Array<int> &coarse_to_ref_type,
2140 Table &ref_type_to_matrix,
2141 Array<Geometry::Type> &ref_type_to_geom)
2142{
2143 const int fine_ne = cft.embeddings.Size();
2144 int coarse_ne = -1;
2145 for (int i = 0; i < fine_ne; i++)
2146 {
2147 coarse_ne = std::max(coarse_ne, cft.embeddings[i].parent);
2148 }
2149 coarse_ne++;
2150
2151 coarse_to_ref_type.SetSize(coarse_ne);
2152 coarse_to_fine.SetDims(coarse_ne, fine_ne);
2153
2154 Array<int> cf_i(coarse_to_fine.GetI(), coarse_ne+1);
2155 Array<Pair<int,int> > cf_j(fine_ne);
2156 cf_i = 0;
2157 for (int i = 0; i < fine_ne; i++)
2158 {
2159 cf_i[cft.embeddings[i].parent+1]++;
2160 }
2161 cf_i.PartialSum();
2162 MFEM_ASSERT(cf_i.Last() == cf_j.Size(), "internal error");
2163 for (int i = 0; i < fine_ne; i++)
2164 {
2165 const Embedding &e = cft.embeddings[i];
2166 cf_j[cf_i[e.parent]].one = e.matrix; // used as sort key below
2167 cf_j[cf_i[e.parent]].two = i;
2168 cf_i[e.parent]++;
2169 }
2170 std::copy_backward(cf_i.begin(), cf_i.end()-1, cf_i.end());
2171 cf_i[0] = 0;
2172 for (int i = 0; i < coarse_ne; i++)
2173 {
2174 std::sort(&cf_j[cf_i[i]], cf_j.GetData() + cf_i[i+1]);
2175 }
2176 for (int i = 0; i < fine_ne; i++)
2177 {
2178 coarse_to_fine.GetJ()[i] = cf_j[i].two;
2179 }
2180
2181 using std::map;
2182 using std::pair;
2183
2184 map<RefType,int> ref_type_map;
2185 for (int i = 0; i < coarse_ne; i++)
2186 {
2187 const int num_children = cf_i[i+1]-cf_i[i];
2188 MFEM_ASSERT(num_children > 0, "");
2189 const int fine_el = cf_j[cf_i[i]].two;
2190 // Assuming the coarse and the fine elements have the same geometry:
2191 const Geometry::Type geom = fine_mesh.GetElementBaseGeometry(fine_el);
2192 const RefType ref_type(geom, num_children, &cf_j[cf_i[i]]);
2193 pair<map<RefType,int>::iterator,bool> res =
2194 ref_type_map.insert(
2195 pair<const RefType,int>(ref_type, (int)ref_type_map.size()));
2196 coarse_to_ref_type[i] = res.first->second;
2197 }
2198
2199 ref_type_to_matrix.MakeI((int)ref_type_map.size());
2200 ref_type_to_geom.SetSize((int)ref_type_map.size());
2201 for (map<RefType,int>::iterator it = ref_type_map.begin();
2202 it != ref_type_map.end(); ++it)
2203 {
2204 ref_type_to_matrix.AddColumnsInRow(it->second, it->first.num_children);
2205 ref_type_to_geom[it->second] = it->first.geom;
2206 }
2207
2208 ref_type_to_matrix.MakeJ();
2209 for (map<RefType,int>::iterator it = ref_type_map.begin();
2210 it != ref_type_map.end(); ++it)
2211 {
2212 const RefType &rt = it->first;
2213 for (int j = 0; j < rt.num_children; j++)
2214 {
2215 ref_type_to_matrix.AddConnection(it->second, rt.children[j].one);
2216 }
2217 }
2218 ref_type_to_matrix.ShiftUpI();
2219}
2220
2221} // namespace internal
2222
2223
2224/// TODO: Implement DofTransformation support
2226 const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes,
2227 BilinearFormIntegrator *mass_integ)
2228 : Operator(c_fes->GetVSize(), f_fes->GetVSize()),
2229 fine_fes(f_fes)
2230{
2231 MFEM_VERIFY(c_fes->GetOrdering() == f_fes->GetOrdering() &&
2232 c_fes->GetVDim() == f_fes->GetVDim(),
2233 "incompatible coarse and fine FE spaces");
2234
2236 Mesh *f_mesh = f_fes->GetMesh();
2237 const CoarseFineTransformations &rtrans = f_mesh->GetRefinementTransforms();
2238
2239 Mesh::GeometryList elem_geoms(*f_mesh);
2241 for (int gi = 0; gi < elem_geoms.Size(); gi++)
2242 {
2243 const Geometry::Type geom = elem_geoms[gi];
2244 DenseTensor &lP = localP[geom], &lM = localM[geom];
2245 const FiniteElement *fine_fe =
2246 f_fes->fec->FiniteElementForGeometry(geom);
2247 const FiniteElement *coarse_fe =
2248 c_fes->fec->FiniteElementForGeometry(geom);
2249 const DenseTensor &pmats = rtrans.point_matrices[geom];
2250
2251 lP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), pmats.SizeK());
2252 lM.SetSize(fine_fe->GetDof(), fine_fe->GetDof(), pmats.SizeK());
2253 emb_tr.SetIdentityTransformation(geom);
2254 for (int i = 0; i < pmats.SizeK(); i++)
2255 {
2256 emb_tr.SetPointMat(pmats(i));
2257 // Get the local interpolation matrix for this refinement type
2258 fine_fe->GetTransferMatrix(*coarse_fe, emb_tr, lP(i));
2259 // Get the local mass matrix for this refinement type
2260 mass_integ->AssembleElementMatrix(*fine_fe, emb_tr, lM(i));
2261 }
2262 }
2263
2264 Table ref_type_to_matrix;
2265 internal::GetCoarseToFineMap(rtrans, *f_mesh, coarse_to_fine,
2266 coarse_to_ref_type, ref_type_to_matrix,
2267 ref_type_to_geom);
2268 MFEM_ASSERT(coarse_to_fine.Size() == c_fes->GetNE(), "");
2269
2270 const int total_ref_types = ref_type_to_geom.Size();
2271 int num_ref_types[Geometry::NumGeom], num_fine_elems[Geometry::NumGeom];
2272 Array<int> ref_type_to_coarse_elem_offset(total_ref_types);
2273 ref_type_to_fine_elem_offset.SetSize(total_ref_types);
2274 std::fill(num_ref_types, num_ref_types+Geometry::NumGeom, 0);
2275 std::fill(num_fine_elems, num_fine_elems+Geometry::NumGeom, 0);
2276 for (int i = 0; i < total_ref_types; i++)
2277 {
2278 Geometry::Type g = ref_type_to_geom[i];
2279 ref_type_to_coarse_elem_offset[i] = num_ref_types[g];
2280 ref_type_to_fine_elem_offset[i] = num_fine_elems[g];
2281 num_ref_types[g]++;
2282 num_fine_elems[g] += ref_type_to_matrix.RowSize(i);
2283 }
2284 DenseTensor localPtMP[Geometry::NumGeom];
2285 for (int g = 0; g < Geometry::NumGeom; g++)
2286 {
2287 if (num_ref_types[g] == 0) { continue; }
2288 const int fine_dofs = localP[g].SizeI();
2289 const int coarse_dofs = localP[g].SizeJ();
2290 localPtMP[g].SetSize(coarse_dofs, coarse_dofs, num_ref_types[g]);
2291 localR[g].SetSize(coarse_dofs, fine_dofs, num_fine_elems[g]);
2292 }
2293 for (int i = 0; i < total_ref_types; i++)
2294 {
2295 Geometry::Type g = ref_type_to_geom[i];
2296 DenseMatrix &lPtMP = localPtMP[g](ref_type_to_coarse_elem_offset[i]);
2297 int lR_offset = ref_type_to_fine_elem_offset[i]; // offset in localR[g]
2298 const int *mi = ref_type_to_matrix.GetRow(i);
2299 const int nm = ref_type_to_matrix.RowSize(i);
2300 lPtMP = 0.0;
2301 for (int s = 0; s < nm; s++)
2302 {
2303 DenseMatrix &lP = localP[g](mi[s]);
2304 DenseMatrix &lM = localM[g](mi[s]);
2305 DenseMatrix &lR = localR[g](lR_offset+s);
2306 MultAtB(lP, lM, lR); // lR = lP^T lM
2307 mfem::AddMult(lR, lP, lPtMP); // lPtMP += lP^T lM lP
2308 }
2309 DenseMatrixInverse lPtMP_inv(lPtMP);
2310 for (int s = 0; s < nm; s++)
2311 {
2312 DenseMatrix &lR = localR[g](lR_offset+s);
2313 lPtMP_inv.Mult(lR); // lR <- (P^T M P)^{-1} P^T M
2314 }
2315 }
2316
2317 // Make a copy of the coarse element-to-dof Table.
2318 coarse_elem_dof = new Table(c_fes->GetElementToDofTable());
2319}
2320
2325
2327 Vector &y) const
2328{
2329 Array<int> c_vdofs, f_vdofs;
2330 Vector loc_x, loc_y;
2331 DenseMatrix loc_x_mat, loc_y_mat;
2332 const int fine_vdim = fine_fes->GetVDim();
2333 const int coarse_ndofs = height/fine_vdim;
2334 for (int coarse_el = 0; coarse_el < coarse_to_fine.Size(); coarse_el++)
2335 {
2336 coarse_elem_dof->GetRow(coarse_el, c_vdofs);
2337 fine_fes->DofsToVDofs(c_vdofs, coarse_ndofs);
2338 loc_y.SetSize(c_vdofs.Size());
2339 loc_y = 0.0;
2340 loc_y_mat.UseExternalData(loc_y.GetData(), c_vdofs.Size()/fine_vdim,
2341 fine_vdim);
2342 const int ref_type = coarse_to_ref_type[coarse_el];
2343 const Geometry::Type geom = ref_type_to_geom[ref_type];
2344 const int *fine_elems = coarse_to_fine.GetRow(coarse_el);
2345 const int num_fine_elems = coarse_to_fine.RowSize(coarse_el);
2346 const int lR_offset = ref_type_to_fine_elem_offset[ref_type];
2347 for (int s = 0; s < num_fine_elems; s++)
2348 {
2349 const DenseMatrix &lR = localR[geom](lR_offset+s);
2350 fine_fes->GetElementVDofs(fine_elems[s], f_vdofs);
2351 x.GetSubVector(f_vdofs, loc_x);
2352 loc_x_mat.UseExternalData(loc_x.GetData(), f_vdofs.Size()/fine_vdim,
2353 fine_vdim);
2354 mfem::AddMult(lR, loc_x_mat, loc_y_mat);
2355 }
2356 y.SetSubVector(c_vdofs, loc_y);
2357 }
2358}
2359
2361 DenseTensor &localR) const
2362{
2363 const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
2364
2365 const CoarseFineTransformations &dtrans =
2367 const DenseTensor &pmats = dtrans.point_matrices[geom];
2368
2369 const int nmat = pmats.SizeK();
2370 const int ldof = fe->GetDof();
2371
2373 isotr.SetIdentityTransformation(geom);
2374
2375 // calculate local restriction matrices for all refinement types
2376 localR.SetSize(ldof, ldof, nmat);
2377 for (int i = 0; i < nmat; i++)
2378 {
2379 isotr.SetPointMat(pmats(i));
2380 fe->GetLocalRestriction(isotr, localR(i));
2381 }
2382}
2383
2385 const Table* old_elem_dof,
2386 const Table* old_elem_fos)
2387{
2388 /// TODO: Implement DofTransformation support
2389
2390 MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
2391 MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
2392 MFEM_VERIFY(ndofs <= old_ndofs, "Previous space is not finer.");
2393
2394 Array<int> dofs, old_dofs, old_vdofs;
2395 Vector row;
2396
2397 Mesh::GeometryList elem_geoms(*mesh);
2398
2400 if (!IsVariableOrder())
2401 {
2402 for (int i = 0; i < elem_geoms.Size(); i++)
2403 {
2404 GetLocalDerefinementMatrices(elem_geoms[i], localR[elem_geoms[i]]);
2405 }
2406 }
2407
2408 SparseMatrix *R = new SparseMatrix(ndofs*vdim, old_ndofs*vdim);
2409
2410 Array<int> mark(R->Height());
2411 mark = 0;
2412
2413 const CoarseFineTransformations &dtrans =
2415
2416 MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
2417
2419 int num_marked = 0;
2420 const FiniteElement *fe = nullptr;
2421 DenseMatrix localRVO; //for variable-order only
2422 for (int k = 0; k < dtrans.embeddings.Size(); k++)
2423 {
2424 const Embedding &emb = dtrans.embeddings[k];
2426
2427 if (IsVariableOrder())
2428 {
2429 fe = GetFE(emb.parent);
2430 const DenseTensor &pmats = dtrans.point_matrices[geom];
2431 const int ldof = fe->GetDof();
2432
2434 isotr.SetIdentityTransformation(geom);
2435
2436 localRVO.SetSize(ldof, ldof);
2437 isotr.SetPointMat(pmats(emb.matrix));
2438 // Local restriction is size ldofxldof assuming that the parent and
2439 // child are of same polynomial order.
2440 fe->GetLocalRestriction(isotr, localRVO);
2441 }
2442 DenseMatrix &lR = IsVariableOrder() ? localRVO : localR[geom](emb.matrix);
2443
2444 elem_dof->GetRow(emb.parent, dofs);
2445 old_elem_dof->GetRow(k, old_dofs);
2446 MFEM_VERIFY(old_dofs.Size() == dofs.Size(),
2447 "Parent and child must have same #dofs.");
2448
2449 for (int vd = 0; vd < vdim; vd++)
2450 {
2451 old_dofs.Copy(old_vdofs);
2452 DofsToVDofs(vd, old_vdofs, old_ndofs);
2453
2454 for (int i = 0; i < lR.Height(); i++)
2455 {
2456 if (!std::isfinite(lR(i, 0))) { continue; }
2457
2458 int r = DofToVDof(dofs[i], vd);
2459 int m = (r >= 0) ? r : (-1 - r);
2460
2461 if (is_dg || !mark[m])
2462 {
2463 lR.GetRow(i, row);
2464 R->SetRow(r, old_vdofs, row);
2465
2466 mark[m] = 1;
2467 num_marked++;
2468 }
2469 }
2470 }
2471 }
2472
2473 if (!is_dg && !IsVariableOrder())
2474 {
2475 MFEM_VERIFY(num_marked == R->Height(),
2476 "internal error: not all rows of R were set.");
2477 }
2478
2479 R->Finalize(); // no-op if fixed width
2480 return R;
2481}
2482
2484 const FiniteElementSpace &coarse_fes, Geometry::Type geom,
2485 DenseTensor &localP) const
2486{
2487 // Assumptions: see the declaration of the method.
2488
2489 const FiniteElement *fine_fe = fec->FiniteElementForGeometry(geom);
2490 const FiniteElement *coarse_fe =
2491 coarse_fes.fec->FiniteElementForGeometry(geom);
2492
2494 const DenseTensor &pmats = rtrans.point_matrices[geom];
2495
2496 int nmat = pmats.SizeK();
2497
2499 isotr.SetIdentityTransformation(geom);
2500
2501 // Calculate the local interpolation matrices for all refinement types
2502 localP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), nmat);
2503 for (int i = 0; i < nmat; i++)
2504 {
2505 isotr.SetPointMat(pmats(i));
2506 fine_fe->GetTransferMatrix(*coarse_fe, isotr, localP(i));
2507 }
2508}
2509
2511 const FiniteElementCollection *fec_,
2512 int vdim_, int ordering_)
2513{
2514 mesh = mesh_;
2515 fec = fec_;
2516 vdim = vdim_;
2517 ordering = (Ordering::Type) ordering_;
2518
2519 elem_dof = NULL;
2520 elem_fos = NULL;
2521 face_dof = NULL;
2522
2523 sequence = 0;
2524 orders_changed = false;
2525 relaxed_hp = false;
2526
2528
2529 const NURBSFECollection *nurbs_fec =
2530 dynamic_cast<const NURBSFECollection *>(fec_);
2531
2532 if (nurbs_fec)
2533 {
2534 MFEM_VERIFY(mesh_->NURBSext, "NURBS FE space requires a NURBS mesh.");
2535
2536 if (NURBSext_ == NULL)
2537 {
2538 NURBSext = mesh_->NURBSext;
2539 own_ext = 0;
2540 }
2541 else
2542 {
2543 NURBSext = NURBSext_;
2544 own_ext = 1;
2545 }
2546 UpdateNURBS();
2547 cP.reset();
2548 cR.reset();
2549 cR_hp.reset();
2550 R_transpose.reset();
2551 cP_is_set = false;
2552
2554 }
2555 else
2556 {
2557 NURBSext = NULL;
2558 own_ext = 0;
2559 Construct();
2560 }
2561
2563}
2564
2566{
2568
2570 for (int i=0; i<DoFTransArray.Size(); i++)
2571 {
2572 DoFTransArray[i] = NULL;
2573 }
2574 if (mesh->Dimension() < 3) { return; }
2575 if (dynamic_cast<const ND_FECollection*>(fec))
2576 {
2577 const FiniteElement *nd_tri =
2579 if (nd_tri)
2580 {
2582 new ND_TriDofTransformation(nd_tri->GetOrder());
2583 }
2584
2585 const FiniteElement *nd_tet =
2587 if (nd_tet)
2588 {
2590 new ND_TetDofTransformation(nd_tet->GetOrder());
2591 }
2592
2593 const FiniteElement *nd_pri =
2595 if (nd_pri)
2596 {
2598 new ND_WedgeDofTransformation(nd_pri->GetOrder());
2599 }
2600
2601 const FiniteElement *nd_pyr =
2603 if (nd_pyr)
2604 {
2607 }
2608 }
2609}
2610
2612{
2613 if (NURBSext && !own_ext)
2614 {
2615 mfem_error("FiniteElementSpace::StealNURBSext");
2616 }
2617 own_ext = 0;
2618
2619 return NURBSext;
2620}
2621
2623{
2624 MFEM_VERIFY(NURBSext, "NURBSExt not defined.");
2625
2626 nvdofs = 0;
2627 nedofs = 0;
2628 nfdofs = 0;
2629 nbdofs = 0;
2630 bdofs = NULL;
2631
2632 delete face_dof;
2633 face_dof = NULL;
2635
2636 // Depending on the element type create the appropriate extensions
2637 // for the individual components.
2638 dynamic_cast<const NURBSFECollection *>(fec)->Reset();
2639
2640 if (dynamic_cast<const NURBS_HDivFECollection *>(fec))
2641 {
2642 VNURBSext.SetSize(mesh->Dimension());
2643 for (int d = 0; d < mesh->Dimension(); d++)
2644 {
2646 }
2647 }
2648
2649 if (dynamic_cast<const NURBS_HCurlFECollection *>(fec))
2650 {
2651 VNURBSext.SetSize(mesh->Dimension());
2652 for (int d = 0; d < mesh->Dimension(); d++)
2653 {
2655 }
2656 }
2657
2658 // If required: concatenate the dof tables of the individual components into
2659 // one dof table for the vector fespace.
2660 if (VNURBSext.Size() == 2)
2661 {
2662 int offset1 = VNURBSext[0]->GetNDof();
2663 ndofs = VNURBSext[0]->GetNDof() + VNURBSext[1]->GetNDof();
2664
2665 // Merge Tables
2666 elem_dof = new Table(*VNURBSext[0]->GetElementDofTable(),
2667 *VNURBSext[1]->GetElementDofTable(),offset1 );
2668
2669 bdr_elem_dof = new Table(*VNURBSext[0]->GetBdrElementDofTable(),
2670 *VNURBSext[1]->GetBdrElementDofTable(),offset1);
2671 }
2672 else if (VNURBSext.Size() == 3)
2673 {
2674 int offset1 = VNURBSext[0]->GetNDof();
2675 int offset2 = offset1 + VNURBSext[1]->GetNDof();
2676 ndofs = offset2 + VNURBSext[2]->GetNDof();
2677
2678 // Merge Tables
2679 elem_dof = new Table(*VNURBSext[0]->GetElementDofTable(),
2680 *VNURBSext[1]->GetElementDofTable(),offset1,
2681 *VNURBSext[2]->GetElementDofTable(),offset2);
2682
2683 bdr_elem_dof = new Table(*VNURBSext[0]->GetBdrElementDofTable(),
2684 *VNURBSext[1]->GetBdrElementDofTable(),offset1,
2685 *VNURBSext[2]->GetBdrElementDofTable(),offset2);
2686 }
2687 else
2688 {
2689 ndofs = NURBSext->GetNDof();
2692 }
2694 sequence++;
2695}
2696
2698{
2699 if (face_dof) { return; }
2700
2701 const int dim = mesh->Dimension();
2702
2703 // Find bdr to face mapping
2705 face_to_be = -1;
2706 for (int b = 0; b < GetNBE(); b++)
2707 {
2709 face_to_be[f] = b;
2710 }
2711
2712 // Loop over faces in correct order, to prevent a sort
2713 // Sort will destroy orientation info in ordering of dofs
2714 Array<Connection> face_dof_list;
2715 Array<int> row;
2716 for (int f = 0; f < GetNF(); f++)
2717 {
2718 int b = face_to_be[f];
2719 if (b == -1) { continue; }
2720 // FIXME: this assumes that the boundary element and the face element have
2721 // the same orientation.
2722 if (dim > 1)
2723 {
2724 const Element *fe = mesh->GetFace(f);
2725 const Element *be = mesh->GetBdrElement(b);
2726 const int nv = be->GetNVertices();
2727 const int *fv = fe->GetVertices();
2728 const int *bv = be->GetVertices();
2729 for (int i = 0; i < nv; i++)
2730 {
2731 MFEM_VERIFY(fv[i] == bv[i],
2732 "non-matching face and boundary elements detected!");
2733 }
2734 }
2735 GetBdrElementDofs(b, row);
2736 Connection conn(f,0);
2737 for (int i = 0; i < row.Size(); i++)
2738 {
2739 conn.to = row[i];
2740 face_dof_list.Append(conn);
2741 }
2742 }
2743 face_dof = new Table(GetNF(), face_dof_list);
2744}
2745
2747{
2748 // This method should be used only for non-NURBS spaces.
2749 MFEM_VERIFY(!NURBSext, "internal error");
2750
2751 // Variable-order space needs a nontrivial P matrix + also ghost elements
2752 // in parallel, we thus require the mesh to be NC.
2753 MFEM_VERIFY(!IsVariableOrder() || Nonconforming(),
2754 "Variable-order space requires a nonconforming mesh.");
2755
2756 elem_dof = NULL;
2757 elem_fos = NULL;
2758 bdr_elem_dof = NULL;
2759 bdr_elem_fos = NULL;
2760 face_dof = NULL;
2761
2762 ndofs = 0;
2763 nvdofs = nedofs = nfdofs = nbdofs = 0;
2764 bdofs = NULL;
2765
2766 cP.reset();
2767 cR.reset();
2768 cR_hp.reset();
2769 cP_is_set = false;
2770 R_transpose.reset();
2771 // 'Th' is initialized/destroyed before this method is called.
2772
2773 int dim = mesh->Dimension();
2774 int order = fec->GetOrder();
2775
2776 MFEM_VERIFY((mesh->GetNumGeometries(dim) > 0) || (mesh->GetNE() == 0),
2777 "Mesh was not correctly finalized.");
2778
2779 bool mixed_elements = (mesh->GetNumGeometries(dim) > 1);
2780 bool mixed_faces = (dim > 2 && mesh->GetNumGeometries(2) > 1);
2781
2782 Array<VarOrderBits> edge_orders, face_orders, edge_elem_orders,
2783 face_elem_orders;
2784
2785 if (IsVariableOrder())
2786 {
2787 // for variable-order spaces, calculate orders of edges and faces
2788 CalcEdgeFaceVarOrders(edge_orders, face_orders, edge_elem_orders,
2789 face_elem_orders, skip_edge, skip_face);
2790 }
2791 else if (mixed_faces)
2792 {
2793 // for mixed faces we also create the var_face_dofs table, see below
2794 face_orders.SetSize(mesh->GetNFaces());
2795 face_orders = (VarOrderBits(1) << order);
2796 }
2797
2798 // assign vertex DOFs
2799 if (mesh->GetNV())
2800 {
2801 nvdofs = mesh->GetNV() * fec->GetNumDof(Geometry::POINT, order);
2802 }
2803
2804 // assign edge DOFs
2805 if (mesh->GetNEdges())
2806 {
2807 if (IsVariableOrder())
2808 {
2809 nedofs = MakeDofTable(1, edge_orders, var_edge_dofs, &var_edge_orders);
2810 MakeDofTable(1, edge_elem_orders, loc_var_edge_dofs,
2812 // Set lnedofs from the last row of loc_var_edge_dofs
2813 Array<int> lastRow;
2815 MFEM_ASSERT(lastRow.Size() == 1, "");
2816 lnedofs = lastRow[0];
2817 }
2818 else
2819 {
2820 // the simple case: all edges are of the same order
2822 var_edge_dofs.Clear(); // ensure any old var_edge_dof table is dumped.
2823 }
2824 }
2825
2826 // assign face DOFs
2827 if (mesh->GetNFaces())
2828 {
2829 if (IsVariableOrder() || mixed_faces)
2830 {
2831 // NOTE: for simplicity, we also use Table var_face_dofs for mixed faces
2832 nfdofs = MakeDofTable(2, face_orders, var_face_dofs,
2833 IsVariableOrder() ? &var_face_orders : NULL);
2834 uni_fdof = -1;
2835
2836 if (IsVariableOrder())
2837 {
2838 MakeDofTable(2, face_elem_orders, loc_var_face_dofs,
2840 // Set lnfdofs from the last row of loc_var_face_dofs
2841 Array<int> lastRow;
2843 MFEM_ASSERT(lastRow.Size() == 1, "");
2844 lnfdofs = lastRow[0];
2845 }
2846 }
2847 else
2848 {
2849 // the simple case: all faces are of the same geometry and order
2852 var_face_dofs.Clear(); // ensure any old var_face_dof table is dumped.
2853 }
2854 }
2855
2856 // assign internal ("bubble") DOFs
2857 if (mesh->GetNE() && dim > 0)
2858 {
2859 if (IsVariableOrder() || mixed_elements)
2860 {
2861 bdofs = new int[mesh->GetNE()+1];
2862 bdofs[0] = 0;
2863 for (int i = 0; i < mesh->GetNE(); i++)
2864 {
2865 int p = GetElementOrderImpl(i);
2867 bdofs[i+1] = nbdofs;
2868 }
2869 }
2870 else
2871 {
2872 // the simple case: all elements are the same
2873 bdofs = NULL;
2875 nbdofs = mesh->GetNE() * fec->GetNumDof(geom, order);
2876 }
2877 }
2878
2880
2882
2883 // record the current mesh sequence number to detect refinement etc.
2885
2886 // increment our sequence number to let GridFunctions know they need updating
2887 sequence++;
2888
2889 // DOFs are now assigned according to current element orders
2890 orders_changed = false;
2891
2892 // Do not build elem_dof Table here: in parallel it has to be constructed
2893 // later.
2894}
2895
2896void DofMapHelper(int entity, const Table & var_ent_dofs,
2897 const Table & loc_var_ent_dofs,
2898 const Array<char> & var_ent_orders,
2899 const Array<char> & loc_var_ent_orders,
2900 Array<int> & all2local, int & ndof_all, int & ndof_loc)
2901{
2902 const int osall0 = var_ent_dofs.GetI()[entity];
2903 const int osall1 = var_ent_dofs.GetI()[entity + 1];
2904
2905 const int osloc0 = loc_var_ent_dofs.GetI()[entity];
2906 const int osloc1 = loc_var_ent_dofs.GetI()[entity + 1];
2907
2908 // loc_var_ent_orders must be a subset of var_ent_orders
2909 int j = osall0;
2910 for (int i=osloc0; i<osloc1; ++i) // Loop over local variants
2911 {
2912 const int order = loc_var_ent_orders[i];
2913 // Find the variant in var_ent_orders with the same order
2914 int na = var_ent_dofs.GetJ()[j + 1] - var_ent_dofs.GetJ()[j];
2915 while (var_ent_orders[j] != order && j < osall1 - 1)
2916 {
2917 j++;
2918 ndof_all += na;
2919 na = var_ent_dofs.GetJ()[j + 1] - var_ent_dofs.GetJ()[j];
2920 }
2921
2922 MFEM_ASSERT(var_ent_orders[j] == order, "");
2923
2924 const int n = loc_var_ent_dofs.GetJ()[i + 1] - loc_var_ent_dofs.GetJ()[i];
2925
2926 MFEM_ASSERT(n == na &&
2927 n == var_ent_dofs.GetJ()[j + 1] - var_ent_dofs.GetJ()[j], "");
2928
2929 for (int k=0; k<n; ++k) { all2local[ndof_all + k] = ndof_loc + k; }
2930
2931 ndof_loc += n;
2932 ndof_all += na;
2933 j++;
2934 }
2935
2936 // Reach the end of all variants for ndof_all
2937 while (j < osall1)
2938 {
2939 const int na = var_ent_dofs.GetJ()[j + 1] - var_ent_dofs.GetJ()[j];
2940 ndof_all += na;
2941 j++;
2942 }
2943}
2944
2946{
2947 if (!IsVariableOrder()) { return; }
2948
2949 // Set a map from all DOFs to local DOFs
2951 all2local = -1;
2952
2953 // Vertex DOFs simply have the identity mapping
2954 for (int i=0; i<nvdofs; ++i)
2955 {
2956 all2local[i] = i;
2957 }
2958
2959 // Redefine local edge DOFs
2960 int ndof_all = nvdofs;
2961 int ndof_loc = nvdofs;
2962 if (mesh->GetNEdges())
2963 {
2964 for (int edge=0; edge<mesh->GetNEdges(); ++edge)
2965 {
2967 loc_var_edge_orders, all2local, ndof_all, ndof_loc);
2968 }
2969
2970 MFEM_ASSERT(ndof_loc - nvdofs == lnedofs, "");
2971 nedofs = lnedofs;
2972 }
2973
2974 // Redefine local face DOFs
2975 if (mesh->GetNFaces())
2976 {
2977 for (int face=0; face<mesh->GetNFaces(); ++face)
2978 {
2980 loc_var_face_orders, all2local, ndof_all, ndof_loc);
2981 }
2982
2983 MFEM_ASSERT(ndof_loc - nvdofs - lnedofs == lnfdofs, "");
2984 nfdofs = lnfdofs;
2985 }
2986
2987 // The remaining DOFs simply have the identity mapping
2988 for (int i=ndof_all; i<ndofs; ++i)
2989 {
2990 all2local[i] = ndof_loc + i - ndof_all;
2991 }
2992
2994}
2995
2997{
2998 MFEM_ASSERT(bits != 0, "invalid bit mask");
2999 for (int order = 0; bits != 0; order++, bits >>= 1)
3000 {
3001 if (bits & 1) { return order; }
3002 }
3003 return 0;
3004}
3005
3006// For the serial FiniteElementSpace, there are no ghost elements, and this
3007// function just sets the sizes of edge_orders and face_orders, initializing to
3008// 0.
3010 Array<VarOrderBits> &edge_orders,
3011 Array<VarOrderBits> &face_orders) const
3012{
3013 edge_orders.SetSize(mesh->GetNEdges());
3014 face_orders.SetSize(mesh->GetNFaces());
3015
3016 edge_orders = 0;
3017 face_orders = 0;
3018}
3019
3021 Array<VarOrderBits> &edge_orders, Array<VarOrderBits> &face_orders,
3022 Array<VarOrderBits> &edge_elem_orders, Array<VarOrderBits> &face_elem_orders,
3023 Array<bool> &skip_edges, Array<bool> &skip_faces) const
3024{
3025 MFEM_ASSERT(Nonconforming(), "");
3026
3027 const bool localVar = elem_order.Size() == mesh->GetNE();
3028 const int baseOrder = fec->GetOrder();
3029
3030 ApplyGhostElementOrdersToEdgesAndFaces(edge_orders, face_orders);
3031
3032 edge_elem_orders.SetSize(mesh->GetNEdges());
3033 face_elem_orders.SetSize(mesh->GetNFaces());
3034
3035 edge_elem_orders = 0;
3036 face_elem_orders = 0;
3037
3040
3043
3044 // Calculate initial edge/face orders, as required by incident elements.
3045 // For each edge/face we accumulate in a bit-mask the orders of elements
3046 // sharing the edge/face.
3047 Array<int> E, F, ori;
3048 for (int i = 0; i < mesh->GetNE(); i++)
3049 {
3050 const int order = localVar ? elem_order[i] : baseOrder;
3051 MFEM_ASSERT(order <= MaxVarOrder, "");
3052 const VarOrderBits mask = (VarOrderBits(1) << order);
3053
3054 mesh->GetElementEdges(i, E, ori);
3055 for (int j = 0; j < E.Size(); j++)
3056 {
3057 edge_orders[E[j]] |= mask;
3058 edge_elem_orders[E[j]] |= mask;
3059
3060 if (order < edge_min_nghb_order[E[j]])
3061 {
3062 edge_min_nghb_order[E[j]] = order;
3063 }
3064 }
3065
3066 if (mesh->Dimension() > 2)
3067 {
3068 mesh->GetElementFaces(i, F, ori);
3069 for (int j = 0; j < F.Size(); j++)
3070 {
3071 face_orders[F[j]] |= mask;
3072 face_elem_orders[F[j]] |= mask;
3073
3074 if (order < face_min_nghb_order[F[j]])
3075 {
3076 face_min_nghb_order[F[j]] = order;
3077 }
3078 }
3079 }
3080 }
3081
3082 if (relaxed_hp)
3083 {
3084 // for relaxed conformity we don't need the masters to match the minimum
3085 // orders of the slaves, we can stop now
3086 return;
3087 }
3088
3089 // Iterate while minimum orders propagate by master/slave relations
3090 // (and new orders also propagate from faces to incident edges).
3091 // See https://github.com/mfem/mfem/pull/1423#issuecomment-638930559
3092 // for an illustration of why this is necessary in hp meshes.
3093 bool done;
3094 do
3095 {
3096 std::set<int> changedEdges;
3097 std::set<int> changedFaces;
3098
3099 const int numEdges = mesh->GetNEdges();
3100
3101 // Propagate from slave edges to master edges
3102 const NCMesh::NCList &edge_list = mesh->ncmesh->GetEdgeList();
3103 for (const NCMesh::Master &master : edge_list.masters)
3104 {
3105 VarOrderBits slave_orders = 0;
3106 for (int i = master.slaves_begin; i < master.slaves_end; i++)
3107 {
3108 slave_orders |= edge_orders[edge_list.slaves[i].index];
3109 }
3110
3111 if (slave_orders == 0)
3112 {
3113 continue;
3114 }
3115
3116 const int min_order_slaves = MinOrder(slave_orders);
3117 if (edge_orders[master.index] == 0 ||
3118 min_order_slaves < MinOrder(edge_orders[master.index]))
3119 {
3120 edge_orders[master.index] |= VarOrderBits(1) << min_order_slaves;
3121 changedEdges.insert(master.index);
3122 }
3123
3124 // Also apply the minimum order to all the slave edges, since they must
3125 // interpolate the master edge, which has the minimum order.
3126 const VarOrderBits min_mask = VarOrderBits(1) << MinOrder(
3127 edge_orders[master.index]);
3128 for (int i = master.slaves_begin; i < master.slaves_end; i++)
3129 {
3130 if (edge_list.slaves[i].index >= numEdges)
3131 {
3132 continue; // Skip ghost edges
3133 }
3134
3135 const VarOrderBits eo0 = edge_orders[edge_list.slaves[i].index];
3136 edge_orders[edge_list.slaves[i].index] |= min_mask;
3137 if (eo0 != edge_orders[edge_list.slaves[i].index])
3138 {
3139 changedEdges.insert(edge_list.slaves[i].index);
3140 }
3141 }
3142 }
3143
3144 // Propagate from slave faces(+edges) to master faces.
3145 const int numFaces = mesh->GetNumFaces();
3146
3147 const NCMesh::NCList &face_list = mesh->ncmesh->GetFaceList();
3148
3149 for (const NCMesh::Master &master : face_list.masters)
3150 {
3151 VarOrderBits slave_orders = 0;
3152
3153 for (int i = master.slaves_begin; i < master.slaves_end; i++)
3154 {
3155 const NCMesh::Slave &slave = face_list.slaves[i];
3156
3157 if (slave.index >= 0)
3158 {
3159 // Note that master.index >= numFaces occurs for ghost master faces.
3160
3161 slave_orders |= face_orders[slave.index];
3162
3163 if (slave.index >= numFaces)
3164 {
3165 continue; // Skip ghost faces
3166 }
3167
3168 mesh->GetFaceEdges(slave.index, E, ori);
3169 for (int j = 0; j < E.Size(); j++)
3170 {
3171 slave_orders |= edge_orders[E[j]];
3172 }
3173 }
3174 else
3175 {
3176 // degenerate face (i.e., edge-face constraint)
3177 slave_orders |= edge_orders[-1 - slave.index];
3178 }
3179 }
3180
3181 if (slave_orders == 0)
3182 {
3183 continue;
3184 }
3185
3186 const int min_order_slaves = MinOrder(slave_orders);
3187 if (face_orders[master.index] == 0 ||
3188 min_order_slaves < MinOrder(face_orders[master.index]))
3189 {
3190 face_orders[master.index] |= VarOrderBits(1) << min_order_slaves;
3191 changedFaces.insert(master.index);
3192 }
3193
3194 // Also apply the minimum order to all the slave faces, since they must
3195 // interpolate the master face, which has the minimum order.
3196 const VarOrderBits min_mask =
3197 VarOrderBits(1) << MinOrder(face_orders[master.index]);
3198 for (int i = master.slaves_begin; i < master.slaves_end; i++)
3199 {
3200 const NCMesh::Slave &slave = face_list.slaves[i];
3201
3202 if (slave.index >= 0 && slave.index < numFaces) // Skip ghost faces
3203 {
3204 const VarOrderBits fo0 = face_orders[slave.index];
3205 face_orders[slave.index] |= min_mask;
3206 if (fo0 != face_orders[slave.index])
3207 {
3208 changedFaces.insert(slave.index);
3209 }
3210 }
3211 }
3212 }
3213
3214 // Make sure edges support (new) orders required by incident faces.
3215 for (int i = 0; i < mesh->GetNFaces(); i++)
3216 {
3217 mesh->GetFaceEdges(i, E, ori);
3218 for (int j = 0; j < E.Size(); j++)
3219 {
3220 const VarOrderBits eo0 = edge_orders[E[j]];
3221 edge_orders[E[j]] |= face_orders[i];
3222 if (eo0 != edge_orders[E[j]])
3223 {
3224 changedEdges.insert(E[j]);
3225 }
3226 }
3227 }
3228
3229 // In the parallel case, OrderPropagation communicates orders on updated
3230 // edges and faces.
3231 done = OrderPropagation(changedEdges, changedFaces,
3232 edge_orders, face_orders);
3233 }
3234 while (!done);
3235
3236 GhostFaceOrderToEdges(face_orders, edge_orders);
3237
3238 // Some ghost edges and faces (3D) may not have any orders applied, since we
3239 // only communicate orders of neighboring ghost elements. Such ghost entities
3240 // are marked here, to be skipped by BuildParallelConformingInterpolation as
3241 // master entities constraining slave entity DOFs.
3242
3243 skip_edges.SetSize(edge_orders.Size());
3244 skip_edges = false;
3245
3246 skip_faces.SetSize(face_orders.Size());
3247 skip_faces = false;
3248
3249 for (int i=0; i<edge_orders.Size(); ++i)
3250 {
3251 if (edge_orders[i] == 0)
3252 {
3253 skip_edges[i] = true;
3254 }
3255 }
3256
3257 for (int i=0; i<face_orders.Size(); ++i)
3258 {
3259 if (face_orders[i] == 0)
3260 {
3261 skip_faces[i] = true;
3262 }
3263 }
3264}
3265
3267 const Array<VarOrderBits> &entity_orders,
3268 Table &entity_dofs,
3269 Array<char> *var_ent_order)
3270{
3271 // The tables var_edge_dofs and var_face_dofs hold DOF assignments for edges
3272 // and faces of a variable-order space, in which each edge/face may host
3273 // several DOF sets, called DOF set variants. Example: an edge 'i' shared by
3274 // 4 hexes of orders 2, 3, 4, 5 will hold four DOF sets, each starting at
3275 // indices e.g. 100, 101, 103, 106, respectively. These numbers are stored
3276 // in row 'i' of var_edge_dofs. Variant zero is always the lowest order DOF
3277 // set, followed by consecutive ranges of higher order DOFs. Variable-order
3278 // faces are handled similarly by var_face_dofs. The tables are empty for
3279 // constant-order spaces.
3280
3281 int num_ent = entity_orders.Size();
3282 int total_dofs = 0;
3283 int total_dofs_nonghost = 0;
3284
3285 Array<Connection> list;
3286 list.Reserve(2*num_ent);
3287
3288 if (var_ent_order)
3289 {
3290 var_ent_order->SetSize(0);
3291 var_ent_order->Reserve(num_ent);
3292 }
3293
3294 int nonGhost = num_ent;
3295 if (IsVariableOrder())
3296 {
3297 nonGhost -= (ent_dim == 1) ? NumGhostEdges() : NumGhostFaces();
3298 }
3299
3300 // assign DOFs according to order bit masks
3301 for (int i = 0; i < num_ent; i++)
3302 {
3303 auto geom = Geometry::SEGMENT; // ent_dim == 1 case
3304 if (ent_dim != 1)
3305 {
3306 // TODO: put this logic in mesh->GetFaceGeometry?
3307 if (i >= nonGhost) // if ghost
3308 {
3309 geom = mesh->ncmesh->GetFaceGeometry(i);
3310 }
3311 else
3312 {
3313 geom = mesh->GetFaceGeometry(i);
3314 }
3315 }
3316
3317 VarOrderBits bits = entity_orders[i];
3318 for (int order = 0; bits != 0; order++, bits >>= 1)
3319 {
3320 if (bits & 1)
3321 {
3322 const int dofs = fec->GetNumDof(geom, order);
3323 list.Append(Connection(i, total_dofs));
3324 total_dofs += dofs;
3325 if (i < nonGhost) { total_dofs_nonghost += dofs; }
3326 if (var_ent_order) { var_ent_order->Append(order); }
3327 }
3328 }
3329 }
3330
3331 // append a dummy row as terminator
3332 list.Append(Connection(num_ent, total_dofs));
3333
3334 // build the table
3335 entity_dofs.MakeFromList(num_ent+1, list);
3336 return total_dofs_nonghost;
3337}
3338
3339int FiniteElementSpace::FindDofs(const Table &var_dof_table,
3340 int row, int ndof) const
3341{
3342 const int *beg = var_dof_table.GetRow(row);
3343 const int *end = var_dof_table.GetRow(row + 1); // terminator, see above
3344
3345 while (beg < end)
3346 {
3347 // return the appropriate range of DOFs
3348 if ((beg[1] - beg[0]) == ndof) { return beg[0]; }
3349 beg++;
3350 }
3351
3352 MFEM_ABORT("DOFs not found for ndof = " << ndof);
3353 return 0;
3354}
3355
3356int FiniteElementSpace::GetEdgeOrder(int edge, int variant) const
3357{
3358 if (!IsVariableOrder()) { return fec->GetOrder(); }
3359
3360 if (edge >= var_edge_dofs.Size())
3361 {
3362 return ghost_edge_orders[edge - var_edge_dofs.Size()];
3363 }
3364
3365 const int* beg = var_edge_dofs.GetRow(edge);
3366 const int* end = var_edge_dofs.GetRow(edge + 1);
3367 if (variant >= end - beg) { return -1; } // past last variant
3368
3369 return var_edge_orders[var_edge_dofs.GetI()[edge] + variant];
3370}
3371
3372int FiniteElementSpace::GetFaceOrder(int face, int variant) const
3373{
3374 if (!IsVariableOrder())
3375 {
3376 // face order can be different from fec->GetOrder()
3377 Geometry::Type geom = mesh->GetFaceGeometry(face);
3378 return fec->FiniteElementForGeometry(geom)->GetOrder();
3379 }
3380
3381 if (face >= var_face_dofs.Size())
3382 {
3383 return ghost_face_orders[face - var_face_dofs.Size()];
3384 }
3385
3386 const int* beg = var_face_dofs.GetRow(face);
3387 const int* end = var_face_dofs.GetRow(face + 1);
3388 if (variant >= end - beg) { return -1; } // past last variant
3389
3390 return var_face_orders[var_face_dofs.GetI()[face] + variant];
3391}
3392
3393int FiniteElementSpace::GetNVariants(int entity, int index) const
3394{
3395 MFEM_ASSERT(IsVariableOrder(), "");
3396 const Table &dof_table = (entity == 1) ? var_edge_dofs : var_face_dofs;
3397
3398 MFEM_ASSERT(index >= 0 && index < dof_table.Size(), "");
3399 return dof_table.GetRow(index + 1) - dof_table.GetRow(index);
3400}
3401
3402static const char* msg_orders_changed =
3403 "Element orders changed, you need to Update() the space first.";
3404
3406 DofTransformation &doftrans) const
3407{
3408 MFEM_VERIFY(!orders_changed, msg_orders_changed);
3409
3410 if (elem_dof)
3411 {
3412 elem_dof->GetRow(elem, dofs);
3413
3415 {
3416 Array<int> Fo;
3417 elem_fos -> GetRow (elem, Fo);
3418 doftrans.SetDofTransformation(
3420 doftrans.SetFaceOrientations(Fo);
3421 doftrans.SetVDim();
3422 }
3423 return;
3424 }
3425
3426 Array<int> V, E, Eo, F, Fo; // TODO: LocalArray
3427
3428 const int dim = mesh->Dimension();
3429 const auto geom = mesh->GetElementGeometry(elem);
3430 const int order = GetElementOrderImpl(elem);
3431
3432 const int nv = fec->GetNumDof(Geometry::POINT, order);
3433 const int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
3434 const int nb = (dim > 0) ? fec->GetNumDof(geom, order) : 0;
3435
3436 if (nv) { mesh->GetElementVertices(elem, V); }
3437 if (ne) { mesh->GetElementEdges(elem, E, Eo); }
3438
3439 int nfd = 0;
3440 if (dim > 2 && fec->HasFaceDofs(geom, order))
3441 {
3442 mesh->GetElementFaces(elem, F, Fo);
3443 for (int i = 0; i < F.Size(); i++)
3444 {
3445 nfd += fec->GetNumDof(mesh->GetFaceGeometry(F[i]), order);
3446 }
3448 {
3449 doftrans.SetDofTransformation(
3451 doftrans.SetFaceOrientations(Fo);
3452 doftrans.SetVDim();
3453 }
3454 }
3455
3456 dofs.SetSize(0);
3457 dofs.Reserve(nv*V.Size() + ne*E.Size() + nfd + nb);
3458
3459 if (nv) // vertex DOFs
3460 {
3461 for (int i = 0; i < V.Size(); i++)
3462 {
3463 for (int j = 0; j < nv; j++)
3464 {
3465 dofs.Append(V[i]*nv + j);
3466 }
3467 }
3468 }
3469
3470 if (ne) // edge DOFs
3471 {
3472 for (int i = 0; i < E.Size(); i++)
3473 {
3474 int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
3475 const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
3476
3477 for (int j = 0; j < ne; j++)
3478 {
3479 dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
3480 }
3481 }
3482 }
3483
3484 if (nfd) // face DOFs
3485 {
3486 for (int i = 0; i < F.Size(); i++)
3487 {
3488 auto fgeom = mesh->GetFaceGeometry(F[i]);
3489 int nf = fec->GetNumDof(fgeom, order);
3490
3491 int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F[i], nf) : F[i]*nf;
3492 const int *ind = fec->GetDofOrdering(fgeom, order, Fo[i]);
3493
3494 for (int j = 0; j < nf; j++)
3495 {
3496 dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
3497 }
3498 }
3499 }
3500
3501 if (nb) // interior ("bubble") DOFs
3502 {
3503 int bbase = bdofs ? bdofs[elem] : elem*nb;
3504 bbase += nvdofs + nedofs + nfdofs;
3505
3506 for (int j = 0; j < nb; j++)
3507 {
3508 dofs.Append(bbase + j);
3509 }
3510 }
3511}
3512
3514 Array<int> &dofs) const
3515{
3517 GetElementDofs(elem, dofs, DoFTrans);
3518 return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
3519}
3520
3522 DofTransformation &doftrans) const
3523{
3524 MFEM_VERIFY(!orders_changed, msg_orders_changed);
3525
3526 if (bdr_elem_dof)
3527 {
3528 bdr_elem_dof->GetRow(bel, dofs);
3529
3531 {
3532 Array<int> Fo;
3533 bdr_elem_fos -> GetRow (bel, Fo);
3534 doftrans.SetDofTransformation(
3536 doftrans.SetFaceOrientations(Fo);
3537 doftrans.SetVDim();
3538 }
3539 return;
3540 }
3541
3542 Array<int> V, E, Eo; // TODO: LocalArray
3543 int F, oF;
3544
3545 int dim = mesh->Dimension();
3546 auto geom = mesh->GetBdrElementGeometry(bel);
3547 int order = fec->GetOrder();
3548
3549 if (elem_order.Size()) // determine order from adjacent element
3550 {
3551 int elem, info;
3552 mesh->GetBdrElementAdjacentElement(bel, elem, info);
3553 order = elem_order[elem];
3554 }
3555
3556 int nv = fec->GetNumDof(Geometry::POINT, order);
3557 int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
3558 int nf = (dim > 2) ? fec->GetNumDof(geom, order) : 0;
3559
3560 if (nv) { mesh->GetBdrElementVertices(bel, V); }
3561 if (ne) { mesh->GetBdrElementEdges(bel, E, Eo); }
3562 if (nf)
3563 {
3564 mesh->GetBdrElementFace(bel, &F, &oF);
3565
3567 {
3568 mfem::Array<int> Fo(1);
3569 Fo[0] = oF;
3570 doftrans.SetDofTransformation(
3572 doftrans.SetFaceOrientations(Fo);
3573 doftrans.SetVDim();
3574 }
3575 }
3576
3577 dofs.SetSize(0);
3578 dofs.Reserve(nv*V.Size() + ne*E.Size() + nf);
3579
3580 if (nv) // vertex DOFs
3581 {
3582 for (int i = 0; i < V.Size(); i++)
3583 {
3584 for (int j = 0; j < nv; j++)
3585 {
3586 dofs.Append(V[i]*nv + j);
3587 }
3588 }
3589 }
3590
3591 if (ne) // edge DOFs
3592 {
3593 for (int i = 0; i < E.Size(); i++)
3594 {
3595 int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
3596 const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
3597
3598 for (int j = 0; j < ne; j++)
3599 {
3600 dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
3601 }
3602 }
3603 }
3604
3605 if (nf) // face DOFs
3606 {
3607 int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F, nf) : F*nf;
3608 const int *ind = fec->GetDofOrdering(geom, order, oF);
3609
3610 for (int j = 0; j < nf; j++)
3611 {
3612 dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
3613 }
3614 }
3615}
3616
3624
3626 int variant) const
3627{
3628 MFEM_VERIFY(!orders_changed, msg_orders_changed);
3629
3630 // If face_dof is already built, use it.
3631 // If it is not and we have a NURBS space, build the face_dof and use it.
3632 if ((face_dof && variant == 0) ||
3633 (NURBSext && (BuildNURBSFaceToDofTable(), true)))
3634 {
3635 face_dof->GetRow(face, dofs);
3636 return fec->GetOrder();
3637 }
3638
3639 int order, nf, fbase;
3640 int dim = mesh->Dimension();
3641 auto fgeom = (dim > 2) ? mesh->GetFaceGeometry(face) : Geometry::INVALID;
3642
3643 if (var_face_dofs.Size() > 0) // variable orders or *mixed* faces
3644 {
3645 const int* beg = var_face_dofs.GetRow(face);
3646 const int* end = var_face_dofs.GetRow(face + 1);
3647 if (variant >= end - beg) { return -1; } // past last face DOFs
3648
3649 fbase = beg[variant];
3650 nf = beg[variant+1] - fbase;
3651
3652 order = !IsVariableOrder() ? fec->GetOrder() :
3653 var_face_orders[var_face_dofs.GetI()[face] + variant];
3654 MFEM_ASSERT(fec->GetNumDof(fgeom, order) == nf, [&]()
3655 {
3656 std::stringstream msg;
3657 msg << "fec->GetNumDof(" << (fgeom == Geometry::SQUARE ? "square" : "triangle")
3658 << ", " << order << ") = " << fec->GetNumDof(fgeom, order) << " nf " << nf;
3659 msg << " face " << face << " variant " << variant << std::endl;
3660 return msg.str();
3661 }());
3662 }
3663 else
3664 {
3665 if (variant > 0) { return -1; }
3666 order = fec->GetOrder();
3667 nf = (dim > 2) ? fec->GetNumDof(fgeom, order) : 0;
3668 fbase = face*nf;
3669 }
3670
3671 // for 1D, 2D and 3D faces
3672 int nv = fec->GetNumDof(Geometry::POINT, order);
3673 int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
3674
3675 Array<int> V, E, Eo;
3676 if (nv) { mesh->GetFaceVertices(face, V); }
3677 if (ne) { mesh->GetFaceEdges(face, E, Eo); }
3678
3679 dofs.SetSize(0);
3680 dofs.Reserve(V.Size() * nv + E.Size() * ne + nf);
3681
3682 if (nv) // vertex DOFs
3683 {
3684 for (int i = 0; i < V.Size(); i++)
3685 {
3686 for (int j = 0; j < nv; j++)
3687 {
3688 dofs.Append(V[i]*nv + j);
3689 }
3690 }
3691 }
3692 if (ne) // edge DOFs
3693 {
3694 for (int i = 0; i < E.Size(); i++)
3695 {
3696 int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
3697 const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
3698
3699 for (int j = 0; j < ne; j++)
3700 {
3701 dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
3702 }
3703 }
3704 }
3705 for (int j = 0; j < nf; j++)
3706 {
3707 dofs.Append(nvdofs + nedofs + fbase + j);
3708 }
3709
3710 return order;
3711}
3712
3714 int variant) const
3715{
3716 MFEM_VERIFY(!orders_changed, msg_orders_changed);
3717
3718 int order, ne, base;
3719 if (IsVariableOrder())
3720 {
3721 const int* beg = var_edge_dofs.GetRow(edge);
3722 const int* end = var_edge_dofs.GetRow(edge + 1);
3723 if (variant >= end - beg) { return -1; } // past last edge DOFs
3724
3725 base = beg[variant];
3726 ne = beg[variant+1] - base;
3727
3728 order = var_edge_orders[var_edge_dofs.GetI()[edge] + variant];
3729 MFEM_ASSERT(fec->GetNumDof(Geometry::SEGMENT, order) == ne, "");
3730 }
3731 else
3732 {
3733 if (variant > 0) { return -1; }
3734 order = fec->GetOrder();
3735 ne = fec->GetNumDof(Geometry::SEGMENT, order);
3736 base = edge*ne;
3737 }
3738
3739 Array<int> V; // TODO: LocalArray
3740 int nv = fec->GetNumDof(Geometry::POINT, order);
3741 if (nv) { mesh->GetEdgeVertices(edge, V); }
3742
3743 dofs.SetSize(0);
3744 dofs.Reserve(2*nv + ne);
3745
3746 for (int i = 0; i < 2; i++)
3747 {
3748 for (int j = 0; j < nv; j++)
3749 {
3750 dofs.Append(V[i]*nv + j);
3751 }
3752 }
3753 for (int j = 0; j < ne; j++)
3754 {
3755 dofs.Append(nvdofs + base + j);
3756 }
3757
3758 return order;
3759}
3760
3762{
3764 dofs.SetSize(nv);
3765 for (int j = 0; j < nv; j++)
3766 {
3767 dofs[j] = i*nv+j;
3768 }
3769}
3770
3772{
3773 MFEM_VERIFY(!orders_changed, msg_orders_changed);
3774
3776 int base = bdofs ? bdofs[i] : i*nb;
3777
3778 dofs.SetSize(nb);
3779 base += nvdofs + nedofs + nfdofs;
3780 for (int j = 0; j < nb; j++)
3781 {
3782 dofs[j] = base + j;
3783 }
3784}
3785
3791
3793{
3794 MFEM_VERIFY(!IsVariableOrder(), "not implemented");
3795
3796 int nf, base;
3797 if (var_face_dofs.Size() > 0) // mixed faces
3798 {
3799 base = var_face_dofs.GetRow(i)[0];
3800 nf = var_face_dofs.GetRow(i)[1] - base;
3801 }
3802 else
3803 {
3804 auto geom = mesh->GetTypicalFaceGeometry();
3805 nf = fec->GetNumDof(geom, fec->GetOrder());
3806 base = i*nf;
3807 }
3808
3809 dofs.SetSize(nf);
3810 for (int j = 0; j < nf; j++)
3811 {
3812 dofs[j] = nvdofs + nedofs + base + j;
3813 }
3814}
3815
3817{
3818 MFEM_VERIFY(!IsVariableOrder(), "not implemented");
3819
3821 dofs.SetSize (ne);
3822 for (int j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
3823 {
3824 dofs[j] = k;
3825 }
3826}
3827
3829{
3830 MFEM_ASSERT(NURBSext,
3831 "FiniteElementSpace::GetPatchDofs needs a NURBSExtension");
3832 NURBSext->GetPatchDofs(patch, dofs);
3833}
3834
3836{
3837 if (i < 0 || i >= mesh->GetNE())
3838 {
3839 if (mesh->GetNE() == 0)
3840 {
3841 MFEM_ABORT("Empty MPI partitions are not permitted!");
3842 }
3843 MFEM_ABORT("Invalid element id:" << i << "; minimum allowed:" << 0 <<
3844 ", maximum allowed:" << mesh->GetNE()-1);
3845 }
3846
3847 const FiniteElement *FE =
3849
3850 if (NURBSext)
3851 {
3852 NURBSext->LoadFE(i, FE);
3853 }
3854 else
3855 {
3856#ifdef MFEM_DEBUG
3857 // consistency check: fec->GetOrder() and FE->GetOrder() should return
3858 // the same value (for standard, constant-order spaces)
3859 if (!IsVariableOrder() && FE->GetDim() > 0)
3860 {
3861 MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
3862 "internal error: " <<
3863 FE->GetOrder() << " != " << fec->GetOrder());
3864 }
3865#endif
3866 }
3867
3868 return FE;
3869}
3870
3872{
3873 if (mesh->GetNE() > 0) { return GetFE(0); }
3874
3876 const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
3877 MFEM_VERIFY(fe != nullptr, "Could not determine a typical FE!");
3878 return fe;
3879}
3880
3882{
3883 int order = fec->GetOrder();
3884
3885 if (IsVariableOrder()) // determine order from adjacent element
3886 {
3887 int elem, info;
3888 mesh->GetBdrElementAdjacentElement(i, elem, info);
3889 order = GetElementOrderImpl(elem);
3890 }
3891
3892 const FiniteElement *BE;
3893 switch (mesh->Dimension())
3894 {
3895 case 1:
3896 BE = fec->GetFE(Geometry::POINT, order);
3897 break;
3898 case 2:
3899 BE = fec->GetFE(Geometry::SEGMENT, order);
3900 break;
3901 case 3:
3902 default:
3903 BE = fec->GetFE(mesh->GetBdrElementGeometry(i), order);
3904 }
3905
3906 if (NURBSext)
3907 {
3908 NURBSext->LoadBE(i, BE);
3909 }
3910
3911 return BE;
3912}
3913
3915{
3916 MFEM_VERIFY(!IsVariableOrder(), "not implemented");
3917
3918 const FiniteElement *fe;
3919 switch (mesh->Dimension())
3920 {
3921 case 1:
3923 break;
3924 case 2:
3926 break;
3927 case 3:
3928 default:
3930 }
3931
3932 if (NURBSext)
3933 {
3934 // Ensure 'face_to_be' is built:
3936 MFEM_ASSERT(face_to_be[i] >= 0,
3937 "NURBS mesh: only boundary faces are supported!");
3938 NURBSext->LoadBE(face_to_be[i], fe);
3939 }
3940
3941 return fe;
3942}
3943
3945 int variant) const
3946{
3947 MFEM_ASSERT(mesh->Dimension() > 1, "No edges with mesh dimension < 2");
3948
3949 int eo = IsVariableOrder() ? GetEdgeOrder(i, variant) : fec->GetOrder();
3950 return fec->GetFE(Geometry::SEGMENT, eo);
3951}
3952
3954 int i, Geometry::Type geom_type) const
3955{
3956 return fec->GetTraceFE(geom_type, GetElementOrder(i));
3957}
3958
3963
3968
3970{
3971 R_transpose.reset();
3972 cR.reset();
3973 cR_hp.reset();
3974 cP.reset();
3975 Th.Clear();
3976 L2E_nat.Clear();
3977 L2E_lex.Clear();
3978 for (int i = 0; i < E2Q_array.Size(); i++)
3979 {
3980 delete E2Q_array[i];
3981 }
3982 E2Q_array.SetSize(0);
3983 for (auto &x : L2F)
3984 {
3985 delete x.second;
3986 }
3987 L2F.clear();
3988 for (int i = 0; i < E2IFQ_array.Size(); i++)
3989 {
3990 delete E2IFQ_array[i];
3991 }
3992 E2IFQ_array.SetSize(0);
3993 for (int i = 0; i < E2BFQ_array.Size(); i++)
3994 {
3995 delete E2BFQ_array[i];
3996 }
3997 E2BFQ_array.SetSize(0);
3998
4000
4005
4006 for (int i = 0; i < VNURBSext.Size(); i++)
4007 {
4008 delete VNURBSext[i];
4009 }
4010
4011 if (NURBSext)
4012 {
4013 if (own_ext) { delete NURBSext; }
4014 delete face_dof;
4016 if (VNURBSext.Size() > 0 )
4017 {
4018 delete elem_dof;
4019 delete bdr_elem_dof;
4020 }
4021 }
4022 else
4023 {
4024 delete elem_dof;
4025 delete elem_fos;
4026 delete bdr_elem_dof;
4027 delete bdr_elem_fos;
4028 delete face_dof;
4029 delete [] bdofs;
4030 }
4032
4033
4034}
4035
4037{
4038 for (int i = 0; i < DoFTransArray.Size(); i++)
4039 {
4040 delete DoFTransArray[i];
4041 }
4042 DoFTransArray.SetSize(0);
4043}
4044
4046 const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
4047{
4048 // Assumptions: see the declaration of the method.
4049
4050 if (T.Type() == Operator::MFEM_SPARSEMAT)
4051 {
4052 if (!IsVariableOrder())
4053 {
4054 Mesh::GeometryList elem_geoms(*mesh);
4055
4057 for (int i = 0; i < elem_geoms.Size(); i++)
4058 {
4059 GetLocalRefinementMatrices(coarse_fes, elem_geoms[i],
4060 localP[elem_geoms[i]]);
4061 }
4062 T.Reset(RefinementMatrix_main(coarse_fes.GetNDofs(),
4063 coarse_fes.GetElementToDofTable(),
4064 coarse_fes.
4066 localP));
4067 }
4068 else
4069 {
4071 coarse_fes.GetElementToDofTable()));
4072 }
4073 }
4074 else
4075 {
4076 T.Reset(new RefinementOperator(this, &coarse_fes));
4077 }
4078}
4079
4081 const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
4082{
4083 const SparseMatrix *coarse_P = coarse_fes.GetConformingProlongation();
4084
4085 Operator::Type req_type = T.Type();
4086 GetTransferOperator(coarse_fes, T);
4087
4088 if (req_type == Operator::MFEM_SPARSEMAT)
4089 {
4091 {
4092 T.Reset(mfem::Mult(*cR, *T.As<SparseMatrix>()));
4093 }
4094 if (coarse_P)
4095 {
4096 T.Reset(mfem::Mult(*T.As<SparseMatrix>(), *coarse_P));
4097 }
4098 }
4099 else
4100 {
4101 const int RP_case = bool(GetConformingRestriction()) + 2*bool(coarse_P);
4102 if (RP_case == 0) { return; }
4103 const bool owner = T.OwnsOperator();
4104 T.SetOperatorOwner(false);
4105 switch (RP_case)
4106 {
4107 case 1:
4108 T.Reset(new ProductOperator(cR.get(), T.Ptr(), false, owner));
4109 break;
4110 case 2:
4111 T.Reset(new ProductOperator(T.Ptr(), coarse_P, owner, false));
4112 break;
4113 case 3:
4115 cR.get(), T.Ptr(), coarse_P, false, owner, false));
4116 break;
4117 }
4118 }
4119}
4120
4122{
4123 Array<char> new_order(mesh->GetNE());
4124 switch (mesh->GetLastOperation())
4125 {
4126 case Mesh::REFINE:
4127 {
4129 for (int i = 0; i < mesh->GetNE(); i++)
4130 {
4131 new_order[i] = elem_order[cf_tr.embeddings[i].parent];
4132 }
4133 break;
4134 }
4135 case Mesh::DEREFINE:
4136 {
4137 const CoarseFineTransformations &cf_tr =
4139 Table coarse_to_fine;
4140 cf_tr.MakeCoarseToFineTable(coarse_to_fine);
4141 Array<int> tabrow;
4142 for (int i = 0; i < coarse_to_fine.Size(); i++)
4143 {
4144 coarse_to_fine.GetRow(i, tabrow);
4145 // For now we require all children to be of same polynomial order.
4146 new_order[i] = elem_order[tabrow[0]];
4147 }
4148 break;
4149 }
4150 default:
4151 MFEM_ABORT("not implemented yet");
4152 }
4153
4154 mfem::Swap(elem_order, new_order);
4155}
4156
4157void FiniteElementSpace::Update(bool want_transform)
4158{
4159 lastUpdatePRef = false;
4160
4161 if (!orders_changed)
4162 {
4163 if (mesh->GetSequence() == mesh_sequence)
4164 {
4165 return; // mesh and space are in sync, no-op
4166 }
4167 if (want_transform && mesh->GetSequence() != mesh_sequence + 1)
4168 {
4169 MFEM_ABORT("Error in update sequence. Space needs to be updated after "
4170 "each mesh modification.");
4171 }
4172 }
4173 else
4174 {
4175 if (mesh->GetSequence() != mesh_sequence)
4176 {
4177 MFEM_ABORT("Updating space after both mesh change and element order "
4178 "change is not supported. Please update separately after "
4179 "each change.");
4180 }
4181 }
4182
4183 if (NURBSext)
4184 {
4185 UpdateNURBS();
4186 return;
4187 }
4188
4189 Table* old_elem_dof = NULL;
4190 Table* old_elem_fos = NULL;
4191 int old_ndofs;
4192 bool old_orders_changed = orders_changed;
4193
4194 // save old DOF table
4195 if (want_transform)
4196 {
4197 old_elem_dof = elem_dof;
4198 old_elem_fos = elem_fos;
4199 elem_dof = NULL;
4200 elem_fos = NULL;
4201 old_ndofs = ndofs;
4202 }
4203
4204 // update the 'elem_order' array if the mesh has changed
4206 {
4208 }
4209
4210 Destroy(); // calls Th.Clear()
4211 Construct();
4213
4214 if (want_transform)
4215 {
4216 MFEM_VERIFY(!old_orders_changed, "Interpolation for element order change "
4217 "is not implemented yet, sorry.");
4218
4219 // calculate appropriate GridFunction transformation
4220 switch (mesh->GetLastOperation())
4221 {
4222 case Mesh::REFINE:
4223 {
4225 {
4226 Th.Reset(new RefinementOperator(this, old_elem_dof,
4227 old_elem_fos, old_ndofs));
4228 // The RefinementOperator takes ownership of 'old_elem_dof', so
4229 // we no longer own it:
4230 old_elem_dof = NULL;
4231 old_elem_fos = NULL;
4232 }
4233 else
4234 {
4235 // calculate fully assembled matrix
4236 Th.Reset(RefinementMatrix(old_ndofs, old_elem_dof,
4237 old_elem_fos));
4238 }
4239 break;
4240 }
4241
4242 case Mesh::DEREFINE:
4243 {
4245 Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof, old_elem_fos));
4246 if (IsVariableOrder())
4247 {
4248 if (cP && cR_hp)
4249 {
4250 Th.SetOperatorOwner(false);
4251 Th.Reset(new TripleProductOperator(cP.get(), cR_hp.get(), Th.Ptr(),
4252 false, false, true));
4253 }
4254 }
4255 else
4256 {
4257 if (cP && cR)
4258 {
4259 Th.SetOperatorOwner(false);
4260 Th.Reset(new TripleProductOperator(cP.get(), cR.get(), Th.Ptr(),
4261 false, false, true));
4262 }
4263 }
4264 break;
4265 }
4266
4267 default:
4268 break;
4269 }
4270
4271 delete old_elem_dof;
4272 delete old_elem_fos;
4273 }
4274}
4275
4277 bool want_transfer)
4278{
4279 MFEM_VERIFY(PRefinementSupported(),
4280 "p-refinement is not supported in this space");
4281
4282 if (want_transfer)
4283 {
4285 for (int i = 0; i<mesh->GetNE(); i++)
4286 {
4287 fesPrev->SetElementOrder(i, GetElementOrder(i));
4288 }
4289 fesPrev->Update(false);
4290 }
4291
4292 for (auto ref : refs)
4293 {
4294 SetElementOrder(ref.index, GetElementOrder(ref.index) + ref.delta);
4295 }
4296
4297 Update(false);
4298
4299 if (want_transfer)
4300 {
4301 PTh.reset(new PRefinementTransferOperator(*fesPrev, *this));
4302 }
4303
4304 lastUpdatePRef = true;
4305}
4306
4308{
4309 // Check whether the space type is L2 or H1
4310 if (!dynamic_cast<const L2_FECollection*>(fec) &&
4311 !dynamic_cast<const H1_FECollection*>(fec))
4312 {
4313 return false;
4314 }
4315
4316 // Check whether the mesh is purely quadrilateral or hexahedral.
4317 const int dim = mesh->Dimension();
4319 mesh->GetGeometries(dim, geoms);
4320 if (geoms.Size() != 1) { return false; }
4321 if (dim == 2 && geoms[0] != Geometry::Type::SQUARE) { return false; }
4322 else if (dim == 3 && geoms[0] != Geometry::Type::CUBE) { return false; }
4323
4324 return true;
4325}
4326
4328{
4329 mesh = new_mesh;
4330}
4331
4333 Vector &fes_node_pos,
4334 int fes_nodes_ordering) const
4335{
4336 Mesh *m = GetMesh();
4337 const int NE = m->GetNE();
4338
4339 if (NE == 0) { fes_node_pos.SetSize(0); return; }
4340
4341 const int dim = m->Dimension();
4342 Array<int> dofs;
4343 Vector e_xyz;
4344 fes_node_pos.SetSize(GetNDofs() * dim);
4345 const FiniteElementSpace *mesh_fes = m->GetNodalFESpace();
4346 FiniteElementSpace vector_fes(m, FEColl(), dim, fes_nodes_ordering);
4347
4348 for (int e = 0; e < NE; e++)
4349 {
4350 mesh_fes->GetElementVDofs(e, dofs);
4351 const int mdof_cnt = dofs.Size() / dim;
4352 mesh_nodes.GetSubVector(dofs, e_xyz); //e_xyz is ordered by nodes here
4353
4354 auto ir = GetFE(e)->GetNodes();
4355 const int fdof_cnt = ir.GetNPoints();
4356 Vector mesh_shape(mdof_cnt), gf_xyz(fdof_cnt * dim);
4357 for (int q = 0; q < fdof_cnt; q++)
4358 {
4359 mesh_fes->GetFE(e)->CalcShape(ir.IntPoint(q), mesh_shape);
4360 for (int d = 0; d < dim; d++)
4361 {
4362 Vector x(e_xyz.GetData() + d*mdof_cnt, mdof_cnt);
4363 gf_xyz(d*fdof_cnt + q) = x * mesh_shape; // order by nodes
4364 }
4365 }
4366
4367 // reuse/resize dofs.
4368 vector_fes.GetElementVDofs(e, dofs);
4369 fes_node_pos.SetSubVector(dofs, gf_xyz);
4370 }
4371}
4372
4373void FiniteElementSpace::Save(std::ostream &os) const
4374{
4375 int fes_format = 90; // the original format, v0.9
4376 bool nurbs_unit_weights = false;
4377
4378 // Determine the format that should be used.
4379 if (!NURBSext)
4380 {
4381 // TODO: if this is a variable-order FE space, use fes_format = 100.
4382 }
4383 else
4384 {
4385 const NURBSFECollection *nurbs_fec =
4386 dynamic_cast<const NURBSFECollection *>(fec);
4387 MFEM_VERIFY(nurbs_fec, "invalid FE collection");
4388 nurbs_fec->SetOrder(NURBSext->GetOrder());
4389 const real_t eps = 5e-14;
4390 nurbs_unit_weights = (NURBSext->GetWeights().Min() >= 1.0-eps &&
4391 NURBSext->GetWeights().Max() <= 1.0+eps);
4393 (NURBSext != mesh->NURBSext && !nurbs_unit_weights) ||
4394 (NURBSext->GetMaster().Size() != 0 ))
4395 {
4396 fes_format = 100; // v1.0 format
4397 }
4398 }
4399
4400 os << (fes_format == 90 ?
4401 "FiniteElementSpace\n" : "MFEM FiniteElementSpace v1.0\n")
4402 << "FiniteElementCollection: " << fec->Name() << '\n'
4403 << "VDim: " << vdim << '\n'
4404 << "Ordering: " << ordering << '\n';
4405
4406 if (fes_format == 100) // v1.0
4407 {
4408 if (!NURBSext)
4409 {
4410 // TODO: this is a variable-order FE space --> write 'element_orders'.
4411 }
4412 else if (NURBSext != mesh->NURBSext)
4413 {
4415 {
4416 os << "NURBS_order\n" << NURBSext->GetOrder() << '\n';
4417 }
4418 else
4419 {
4420 os << "NURBS_orders\n";
4421 // 1 = do not write the size, just the entries:
4422 NURBSext->GetOrders().Save(os, 1);
4423 }
4424 // If periodic BCs are given, write connectivity
4425 if (NURBSext->GetMaster().Size() != 0 )
4426 {
4427 os <<"NURBS_periodic\n";
4428 NURBSext->GetMaster().Save(os);
4429 NURBSext->GetSlave().Save(os);
4430 }
4431 // If the weights are not unit, write them to the output:
4432 if (!nurbs_unit_weights)
4433 {
4434 os << "NURBS_weights\n";
4435 NURBSext->GetWeights().Print(os, 1);
4436 }
4437 }
4438 os << "End: MFEM FiniteElementSpace v1.0\n";
4439 }
4440}
4441
4442std::shared_ptr<const PRefinementTransferOperator>
4444
4445void FiniteElementSpace
4446::GetEssentialBdrEdgesFaces(const Array<int> &bdr_attr_is_ess,
4447 std::set<int> & edges, std::set<int> & faces) const
4448{
4449 const int dim = mesh->Dimension();
4450 MFEM_VERIFY(dim == 2 || dim == 3, "");
4451
4452 for (int i = 0; i < GetNBE(); i++)
4453 {
4454 if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
4455 {
4456 int f, o;
4457 mesh->GetBdrElementFace(i, &f, &o);
4458
4459 if (dim == 3)
4460 {
4461 faces.insert(f);
4462 Array<int> edges_i, cor;
4463 mesh->GetBdrElementEdges(i, edges_i, cor);
4464 for (auto edge : edges_i)
4465 {
4466 edges.insert(edge);
4467 }
4468 }
4469 else
4470 {
4471 edges.insert(f);
4472 }
4473 }
4474 }
4475
4476 if (Nonconforming())
4477 {
4478 Array<int> bdr_verts, bdr_edges, bdr_faces;
4479 mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges,
4480 bdr_faces);
4481
4482 for (auto e : bdr_edges)
4483 {
4484 edges.insert(e);
4485 }
4486
4487 for (auto f : bdr_faces)
4488 {
4489 faces.insert(f);
4490 }
4491 }
4492}
4493
4495{
4496 string buff;
4497 int fes_format = 0, ord;
4499
4500 Destroy();
4501
4502 input >> std::ws;
4503 getline(input, buff); // 'FiniteElementSpace'
4504 filter_dos(buff);
4505 if (buff == "FiniteElementSpace") { fes_format = 90; /* v0.9 */ }
4506 else if (buff == "MFEM FiniteElementSpace v1.0") { fes_format = 100; }
4507 else { MFEM_ABORT("input stream is not a FiniteElementSpace!"); }
4508 getline(input, buff, ' '); // 'FiniteElementCollection:'
4509 input >> std::ws;
4510 getline(input, buff);
4511 filter_dos(buff);
4512 r_fec = FiniteElementCollection::New(buff.c_str());
4513 getline(input, buff, ' '); // 'VDim:'
4514 input >> vdim;
4515 getline(input, buff, ' '); // 'Ordering:'
4516 input >> ord;
4517
4518 NURBSFECollection *nurbs_fec = dynamic_cast<NURBSFECollection*>(r_fec);
4519 if (nurbs_fec) { nurbs_fec->SetDim(m->Dimension()); }
4520 NURBSExtension *nurbs_ext = NULL;
4521 if (fes_format == 90) // original format, v0.9
4522 {
4523 if (nurbs_fec)
4524 {
4525 MFEM_VERIFY(m->NURBSext, "NURBS FE collection requires a NURBS mesh!");
4526 const int order = nurbs_fec->GetOrder();
4527 if (order != m->NURBSext->GetOrder() &&
4529 {
4530 nurbs_ext = new NURBSExtension(m->NURBSext, order);
4531 }
4532 }
4533 }
4534 else if (fes_format == 100) // v1.0
4535 {
4536 while (1)
4537 {
4538 skip_comment_lines(input, '#');
4539 MFEM_VERIFY(input.good(), "error reading FiniteElementSpace v1.0");
4540 getline(input, buff);
4541 filter_dos(buff);
4542 if (buff == "NURBS_order" || buff == "NURBS_orders")
4543 {
4544 MFEM_VERIFY(nurbs_fec,
4545 buff << ": NURBS FE collection is required!");
4546 MFEM_VERIFY(m->NURBSext, buff << ": NURBS mesh is required!");
4547 MFEM_VERIFY(!nurbs_ext, buff << ": order redefinition!");
4548 if (buff == "NURBS_order")
4549 {
4550 int order;
4551 input >> order;
4552 nurbs_ext = new NURBSExtension(m->NURBSext, order);
4553 }
4554 else
4555 {
4556 Array<int> orders;
4557 orders.Load(m->NURBSext->GetNKV(), input);
4558 nurbs_ext = new NURBSExtension(m->NURBSext, orders);
4559 }
4560 }
4561 else if (buff == "NURBS_periodic")
4562 {
4563 Array<int> master, slave;
4564 master.Load(input);
4565 slave.Load(input);
4566 nurbs_ext->ConnectBoundaries(master,slave);
4567 }
4568 else if (buff == "NURBS_weights")
4569 {
4570 MFEM_VERIFY(nurbs_ext, "NURBS_weights: NURBS_orders have to be "
4571 "specified before NURBS_weights!");
4572 nurbs_ext->GetWeights().Load(input, nurbs_ext->GetNDof());
4573 }
4574 else if (buff == "element_orders")
4575 {
4576 MFEM_VERIFY(!nurbs_fec, "section element_orders cannot be used "
4577 "with a NURBS FE collection");
4578 MFEM_ABORT("element_orders: not implemented yet!");
4579 }
4580 else if (buff == "End: MFEM FiniteElementSpace v1.0")
4581 {
4582 break;
4583 }
4584 else
4585 {
4586 MFEM_ABORT("unknown section: " << buff);
4587 }
4588 }
4589 }
4590
4591 Constructor(m, nurbs_ext, r_fec, vdim, ord);
4592
4593 return r_fec;
4594}
4595
4602
4603} // namespace mfem
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void Load(std::istream &in, int fmt=0)
Read an Array from the stream in using format fmt. The format fmt can be:
Definition array.cpp:53
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition array.hpp:341
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition array.hpp:165
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
int Size() const
Return the logical size of the array.
Definition array.hpp:147
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
Definition array.hpp:943
void DeleteAll()
Delete the whole array.
Definition array.hpp:925
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
void Save(std::ostream &out, int fmt=0) const
Save the Array to the stream out using the format fmt. The format fmt can be:
Definition array.cpp:40
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition array.hpp:935
T Sum() const
Return the sum of all the array entries using the '+'' operator for class 'T'.
Definition array.cpp:115
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
Definition array.hpp:349
Abstract base class BilinearFormIntegrator.
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Given a particular Finite Element computes the element matrix elmat.
Operator that extracts face degrees of freedom for H1, ND, or RT FiniteElementSpaces.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:125
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:172
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
void UseExternalData(real_t *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition densemat.hpp:80
void GetRow(int r, Vector &row) const
Rank 3 tensor (array of matrices)
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
int SizeJ() const
int SizeI() const
int SizeK() const
void SetFaceOrientations(const Array< int > &Fo)
Configure the transformation using face orientations for the current element.
Definition doftrans.hpp:169
void SetDofTransformation(const StatelessDofTransformation &dof_trans)
Set or change the nested StatelessDofTransformation object.
Definition doftrans.hpp:176
const StatelessDofTransformation * GetDofTransformation() const
Return the nested StatelessDofTransformation object.
Definition doftrans.hpp:186
void SetVDim(int vdim=1, int ordering=0)
Set or change the vdim and ordering parameter.
Definition doftrans.hpp:190
void TransformPrimal(real_t *v) const
Definition doftrans.cpp:17
Abstract base class that defines an interface for element restrictions.
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
virtual int GetNVertices() const =0
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
const IntegrationRule * IntRule
Not owned.
Base class for operators that extracts Face degrees of freedom.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name.
Definition fe_coll.cpp:124
const int * GetDofOrdering(Geometry::Type geom, int p, int ori) const
Variable order version of DofOrderForOrientation().
Definition fe_coll.hpp:230
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition fe_coll.hpp:240
virtual int GetContType() const =0
int HasFaceDofs(Geometry::Type geom, int p) const
Definition fe_coll.cpp:100
virtual int DofForGeometry(Geometry::Type GeomType) const =0
int GetNumDof(Geometry::Type geom, int p) const
Variable order version of DofForGeometry().
Definition fe_coll.hpp:218
virtual const FiniteElement * TraceFiniteElementForGeometry(Geometry::Type GeomType) const
Definition fe_coll.hpp:97
const FiniteElement * GetTraceFE(Geometry::Type geom, int p) const
Variable order version of TraceFiniteElementForGeometry().
Definition fe_coll.hpp:206
virtual const char * Name() const
Definition fe_coll.hpp:79
const FiniteElement * GetFE(Geometry::Type geom, int p) const
Variable order version of FiniteElementForGeometry().
Definition fe_coll.hpp:195
@ DISCONTINUOUS
Field is discontinuous across element interfaces.
Definition fe_coll.hpp:48
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
Definition fespace.cpp:2326
DerefinementOperator(const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes, BilinearFormIntegrator *mass_integ)
TODO: Implement DofTransformation support.
Definition fespace.cpp:2225
GridFunction interpolation operator applicable after mesh refinement.
Definition fespace.hpp:532
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition fespace.cpp:2001
RefinementOperator(const FiniteElementSpace *fespace, Table *old_elem_dof, Table *old_elem_fos, int old_ndofs)
Definition fespace.cpp:1814
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition fespace.cpp:1923
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
std::unique_ptr< FiniteElementSpace > fesPrev
Definition fespace.hpp:344
void Save(std::ostream &out) const
Save finite element space to output stream out.
Definition fespace.cpp:4373
virtual void ApplyGhostElementOrdersToEdgesAndFaces(Array< VarOrderBits > &edge_orders, Array< VarOrderBits > &face_orders) const
Helper function for ParFiniteElementSpace.
Definition fespace.cpp:3009
int GetEntityVDofs(int entity, int index, Array< int > &dofs, Geometry::Type master_geom=Geometry::INVALID, int variant=0) const
Helper to get vertex, edge or face VDOFs (entity=0,1,2 resp.).
Definition fespace.cpp:1118
void GetVDofs(int vd, Array< int > &dofs, int ndofs=-1) const
Returns the indices of all of the VDofs for the specified dimension 'vd'.
Definition fespace.cpp:238
DofTransformation DoFTrans
Definition fespace.hpp:328
static int EncodeDof(int entity_base, int idx)
Helper to encode a sign flip into a DOF index (for Hcurl/Hdiv shapes).
Definition fespace.hpp:1164
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1463
void ReorderElementToDofTable()
Reorder the scalar DOFs based on the element ordering.
Definition fespace.cpp:506
Array< char > var_face_orders
Definition fespace.hpp:296
int GetVectorDim() const
Return the total dimension of a vector in the space.
Definition fespace.cpp:1490
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:709
void SetRestriction(const SparseMatrix &r)
Definition fespace.cpp:178
void BuildNURBSFaceToDofTable() const
Generates partial face_dof table for a NURBS space.
Definition fespace.cpp:2697
void BuildDofToBdrArrays() const
Initialize internal data that enables the use of the methods GetBdrElementForDof() and GetBdrLocalDof...
Definition fespace.cpp:551
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I, int skipfirst=0)
Definition fespace.cpp:949
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element,...
Definition fespace.hpp:1287
Array< int > dof_ldof_array
Definition fespace.hpp:315
Array< StatelessDofTransformation * > DoFTransArray
Definition fespace.hpp:327
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
Definition fespace.cpp:258
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
Definition fespace.cpp:3816
Array< FaceQuadratureInterpolator * > E2BFQ_array
Definition fespace.hpp:373
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition fespace.cpp:3881
Array< bool > skip_face
Definition fespace.hpp:305
Array< FaceQuadratureInterpolator * > E2IFQ_array
Definition fespace.hpp:372
Array< int > face_min_nghb_order
Definition fespace.hpp:301
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
Definition fespace.cpp:3513
int GetNumElementInteriorDofs(int i) const
Returns the number of degrees of freedom associated with the interior of the specified element.
Definition fespace.cpp:3786
std::shared_ptr< PRefinementTransferOperator > PTh
Definition fespace.hpp:349
virtual void GetExteriorTrueDofs(Array< int > &exterior_dofs, int component=-1) const
Get a list of all true dofs on the exterior of the mesh, exterior_dofs. For spaces with 'vdim' > 1,...
Definition fespace.cpp:746
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
Definition fespace.cpp:367
NURBSExtension * NURBSext
Definition fespace.hpp:319
virtual int GetFaceDofs(int face, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition fespace.cpp:3625
virtual void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer true-dof data from coarse_fes,...
Definition fespace.cpp:4080
static void AdjustVDofs(Array< int > &vdofs)
Remove the orientation information encoded into an array of dofs Some basis function types have a rel...
Definition fespace.cpp:310
virtual void GetExteriorVDofs(Array< int > &exterior_vdofs, int component=-1) const
Mark degrees of freedom associated with exterior faces of the mesh. For spaces with 'vdim' > 1,...
Definition fespace.cpp:717
SparseMatrix * VariableOrderRefinementMatrix(const int coarse_ndofs, const Table &coarse_elem_dof) const
Definition fespace.cpp:1706
int GetEdgeOrder(int edge, int variant=0) const
Definition fespace.cpp:3356
Array< char > loc_var_face_orders
Definition fespace.hpp:297
bool Nonconforming() const
Definition fespace.hpp:686
void GetVertexVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified vertices.
Definition fespace.cpp:373
void GetVertexDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the specified vertices.
Definition fespace.cpp:3761
OperatorHandle L2E_lex
Definition fespace.hpp:355
void BuildFaceToDofTable() const
Definition fespace.cpp:471
int GetFaceOrder(int face, int variant=0) const
Returns the polynomial degree of the i'th face finite element.
Definition fespace.cpp:3372
int GetNumBorderDofs(Geometry::Type geom, int order) const
Definition fespace.cpp:1084
FiniteElementSpace()
Default constructor: the object is invalid until initialized using the method Load().
Definition fespace.cpp:59
Array< int > dof_elem_array
Definition fespace.hpp:314
int FindFaceDof(int face, int ndof) const
Similar to FindEdgeDof, but used for mixed meshes too.
Definition fespace.hpp:483
static int MinOrder(VarOrderBits bits)
Return the minimum order (least significant bit set) in the bit mask.
Definition fespace.cpp:2996
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:809
MFEM_DEPRECATED void RebuildElementToDofTable()
(
Definition fespace.cpp:497
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:845
bool orders_changed
True if at least one element order changed (variable-order space only).
Definition fespace.hpp:384
SparseMatrix * DerefinementMatrix(int old_ndofs, const Table *old_elem_dof, const Table *old_elem_fos)
Calculate GridFunction restriction matrix after mesh derefinement.
Definition fespace.cpp:2384
friend class PRefinementTransferOperator
Definition fespace.hpp:246
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
Definition fespace.cpp:1767
void GetEdgeInteriorVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
Definition fespace.cpp:385
void GetTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer GridFunction data from coarse_fes,...
Definition fespace.cpp:4045
SparseMatrix * D2C_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Generate the global restriction matrix from a discontinuous FE space to the continuous FE space of th...
Definition fespace.cpp:839
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:658
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition fespace.hpp:900
virtual void UpdateMeshPointer(Mesh *new_mesh)
Definition fespace.cpp:4327
int GetDegenerateFaceDofs(int index, Array< int > &dofs, Geometry::Type master_geom, int variant) const
Definition fespace.cpp:1042
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition fespace.cpp:1582
Array< char > ghost_edge_orders
Definition fespace.hpp:298
int GetBdrAttribute(int i) const
Definition fespace.hpp:938
Array< int > dof_bdr_elem_array
Definition fespace.hpp:316
int GetEntityDofs(int entity, int index, Array< int > &dofs, Geometry::Type master_geom=Geometry::INVALID, int variant=0) const
Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
Definition fespace.cpp:1093
virtual const Operator * GetRestrictionOperator() const
An abstract operator that performs the same action as GetRestrictionMatrix.
Definition fespace.hpp:746
static constexpr int MaxVarOrder
Definition fespace.hpp:292
virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes, const Array< int > *perm)
Copies the prolongation and restriction matrices from fes.
Definition fespace.cpp:108
Array< char > var_edge_orders
Definition fespace.hpp:296
std::unique_ptr< Operator > R_transpose
Operator computing the action of the transpose of the restriction.
Definition fespace.hpp:340
bool Conforming() const
Definition fespace.hpp:685
int MakeDofTable(int ent_dim, const Array< VarOrderBits > &entity_orders, Table &entity_dofs, Array< char > *var_ent_order)
Definition fespace.cpp:3266
void UpdateElementOrders()
Resize the elem_order array on mesh change.
Definition fespace.cpp:4121
void MakeVDimMatrix(SparseMatrix &mat) const
Replicate 'mat' in the vector dimension, according to vdim ordering mode.
Definition fespace.cpp:1428
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition fespace.hpp:897
SparseMatrix * H2L_GlobalRestrictionMatrix(FiniteElementSpace *lfes)
Construct the restriction matrix from the FE space given by (*this) to the lower degree FE space give...
Definition fespace.cpp:902
SparseMatrix * D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Generate the global restriction matrix from a discontinuous FE space to the piecewise constant FE spa...
Definition fespace.cpp:871
const FiniteElementCollection * fec
Associated FE collection (not owned).
Definition fespace.hpp:255
void GetNodePositions(const Vector &mesh_nodes, Vector &fes_node_pos, int fes_nodes_ordering=Ordering::byNODES) const
Compute the space's node positions w.r.t. given mesh positions. The function uses FiniteElement::GetN...
Definition fespace.cpp:4332
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller.
Definition fespace.cpp:4494
void BuildBdrElementToDofTable() const
Definition fespace.cpp:432
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:332
Array< QuadratureInterpolator * > E2Q_array
Definition fespace.hpp:371
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3835
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:876
NURBSExtension * StealNURBSext()
Definition fespace.cpp:2611
void GetPatchDofs(int patch, Array< int > &dofs) const
Returns indices of degrees of freedom for NURBS patch index patch. Cartesian ordering is used,...
Definition fespace.cpp:3828
void BuildDofToArrays_() const
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
Definition fespace.cpp:526
const Operator * GetRestrictionTransposeOperator() const
Return an operator that performs the transpose of GetRestrictionOperator.
Definition fespace.cpp:1478
Array< char > elem_order
Definition fespace.hpp:272
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Definition fespace.cpp:1028
Array< int > face_to_be
Definition fespace.hpp:325
void GetLocalDerefinementMatrices(Geometry::Type geom, DenseTensor &localR) const
Definition fespace.cpp:2360
void BuildConformingInterpolation() const
Calculate the cP and cR matrices for a nonconforming mesh.
Definition fespace.cpp:1178
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:891
int vdim
Vector dimension (number of unknowns per degree of freedom).
Definition fespace.hpp:258
Table var_face_dofs
NOTE: also used for spaces with mixed faces.
Definition fespace.hpp:281
OperatorHandle L2E_nat
The element restriction operators, see GetElementRestriction().
Definition fespace.hpp:355
bool lastUpdatePRef
Flag to indicate whether the last update was for p-refinement.
Definition fespace.hpp:352
int * bdofs
internal DOFs of elements if mixed/var-order; NULL otherwise
Definition fespace.hpp:276
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1510
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
Definition fespace.cpp:3944
Array< bool > skip_edge
Definition fespace.hpp:305
Array< char > ghost_face_orders
Definition fespace.hpp:298
int GetEdgeDofs(int edge, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
Definition fespace.cpp:3713
std::tuple< bool, ElementDofOrdering, FaceType, L2FaceValues > key_face
The face restriction operators, see GetFaceRestriction().
Definition fespace.hpp:357
std::unique_ptr< SparseMatrix > cR
Conforming restriction matrix such that cR.cP=I.
Definition fespace.hpp:335
void GetElementInteriorVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the interior of the specified element.
Definition fespace.cpp:379
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition fespace.cpp:361
Array< char > loc_var_edge_orders
Definition fespace.hpp:297
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition fespace.cpp:4157
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
Definition fespace.hpp:347
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof, const Table *old_elem_fos)
Definition fespace.cpp:1790
std::uint64_t VarOrderBits
Bit-mask representing a set of orders needed by an edge/face.
Definition fespace.hpp:291
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
Definition fespace.cpp:1642
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
Definition fespace.hpp:266
std::unique_ptr< SparseMatrix > cP
Definition fespace.hpp:333
void AddEdgeFaceDependencies(SparseMatrix &deps, Array< int > &master_dofs, const FiniteElement *master_fe, Array< int > &slave_dofs, int slave_face, const DenseMatrix *pm) const
Definition fespace.cpp:974
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
Definition fespace.cpp:3953
const FiniteElement * GetTypicalTraceElement() const
Return a "typical" trace element.
Definition fespace.cpp:3959
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:790
Array< NURBSExtension * > VNURBSext
Definition fespace.hpp:323
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
Definition fespace.cpp:221
void SetVarOrderLocalDofs()
Sets all2local. See documentation of all2local for details.
Definition fespace.cpp:2945
Mesh * mesh
The mesh that FE space lives on (not owned).
Definition fespace.hpp:252
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
Definition fespace.cpp:3914
Ordering::Type ordering
Definition fespace.hpp:263
void ConvertToConformingVDofs(const Array< int > &dofs, Array< int > &cdofs)
For a partially conforming FE space, convert a marker array (nonzero entries are true) on the partial...
Definition fespace.cpp:822
void CalcEdgeFaceVarOrders(Array< VarOrderBits > &edge_orders, Array< VarOrderBits > &face_orders, Array< VarOrderBits > &edge_elem_orders, Array< VarOrderBits > &face_elem_orders, Array< bool > &skip_edges, Array< bool > &skip_faces) const
Definition fespace.cpp:3020
virtual bool OrderPropagation(const std::set< int > &edges, const std::set< int > &faces, Array< VarOrderBits > &edge_orders, Array< VarOrderBits > &face_orders) const
Returns true if order propagation is done, for variable-order spaces.
Definition fespace.hpp:457
std::shared_ptr< const PRefinementTransferOperator > GetPrefUpdateOperator()
Definition fespace.cpp:4443
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
Definition fespace.cpp:196
const SparseMatrix * GetHpConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1471
const SparseMatrix * GetConformingProlongation() const
Definition fespace.cpp:1456
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:878
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
int GetElementOrderImpl(int i) const
Return element order: internal version of GetElementOrder without checks.
Definition fespace.cpp:232
int GetNVariants(int entity, int index) const
Return number of possible DOF variants for edge/face (var. order spaces).
Definition fespace.cpp:3393
Array< int > dof_bdr_ldof_array
Definition fespace.hpp:317
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:848
const Table * GetElementToFaceOrientationTable() const
Definition fespace.hpp:1283
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
Definition fespace.cpp:3617
void GetBoundaryTrueDofs(Array< int > &boundary_dofs, int component=-1)
Get a list of all boundary true dofs, boundary_dofs. For spaces with 'vdim' > 1, the 'component' para...
Definition fespace.cpp:702
int GetVDim() const
Returns the vector dimension of the finite element space.
Definition fespace.hpp:841
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...
Definition fespace.cpp:1611
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Definition fespace.hpp:1510
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified face.
Definition fespace.cpp:3792
virtual void PRefineAndUpdate(const Array< pRefinement > &refs, bool want_transfer=true)
Definition fespace.cpp:4276
int GetNConformingDofs() const
Definition fespace.cpp:1484
void SetProlongation(const SparseMatrix &p)
Definition fespace.cpp:159
Array< int > edge_min_nghb_order
Minimum order among neighboring elements.
Definition fespace.hpp:301
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3871
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Definition fespace.hpp:1168
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified element.
Definition fespace.cpp:3771
virtual int NumGhostEdges() const
Returns the number of ghost edges (nonzero in ParFiniteElementSpace).
Definition fespace.hpp:464
void Constructor(Mesh *mesh, NURBSExtension *ext, const FiniteElementCollection *fec, int vdim=1, int ordering=Ordering::byNODES)
Help function for constructors + Load().
Definition fespace.cpp:2510
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
Definition fespace.cpp:348
void ConvertFromConformingVDofs(const Array< int > &cdofs, Array< int > &dofs)
For a partially conforming FE space, convert a marker array (nonzero entries are true) on the conform...
Definition fespace.cpp:830
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
Definition fespace.cpp:584
virtual int NumGhostFaces() const
Returns the number of ghost faces (nonzero in ParFiniteElementSpace).
Definition fespace.hpp:467
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition fespace.cpp:294
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering f_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Definition fespace.cpp:1543
int GetCurlDim() const
Return the dimension of the curl of a GridFunction defined on this space.
Definition fespace.cpp:1500
virtual void GhostFaceOrderToEdges(const Array< VarOrderBits > &face_orders, Array< VarOrderBits > &edge_orders) const
Helper function for ParFiniteElementSpace.
Definition fespace.hpp:452
int FindEdgeDof(int edge, int ndof) const
Definition fespace.hpp:479
void GetPatchVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom in vdofs for NURBS patch i.
Definition fespace.cpp:355
void BuildElementToDofTable() const
Definition fespace.cpp:391
std::unique_ptr< SparseMatrix > cR_hp
A version of the conforming restriction matrix for variable-order spaces.
Definition fespace.hpp:337
int FindDofs(const Table &var_dof_table, int row, int ndof) const
Search row of a DOF table for a DOF set of size 'ndof', return first DOF.
Definition fespace.cpp:3339
void VariableOrderMinimumRule(SparseMatrix &deps) const
Definition fespace.cpp:1128
Abstract class for all finite elements.
Definition fe_base.hpp:244
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
Definition fe_base.hpp:325
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:338
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:321
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition fe_base.cpp:119
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition fe_base.hpp:351
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
Definition fe_base.cpp:113
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:331
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition fe_base.cpp:107
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_base.cpp:126
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
Definition fe_base.hpp:328
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
static const int NumGeom
Definition geom.hpp:46
static const int NumVerts[NumGeom]
Definition geom.hpp:53
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:275
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
A standard isoparametric element transformation.
Definition eltrans.hpp:629
void SetPointMat(const DenseMatrix &pm)
Set the underlying point matrix describing the transformation.
Definition eltrans.hpp:668
void SetFE(const FiniteElement *FE)
Set the element that will be used to compute the transformations.
Definition eltrans.hpp:648
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition eltrans.cpp:417
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition eltrans.hpp:671
Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
Operator that extracts Face degrees of freedom for L2 spaces.
Operator that extracts face degrees of freedom for L2 interface spaces.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:346
Class used by MFEM to store pointers to host and/or device memory.
List of mesh geometries stored as Array<Geometry::Type>.
Definition mesh.hpp:1489
Mesh data type.
Definition mesh.hpp:64
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
Definition mesh.cpp:7321
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
Definition mesh.cpp:7254
Operation GetLastOperation() const
Return type of last modification of the mesh.
Definition mesh.hpp:2381
int GetNEdges() const
Return the number of edges.
Definition mesh.hpp:1288
void GetBdrElementFace(int i, int *f, int *o) const
Definition mesh.cpp:7565
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:290
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:298
virtual void GetExteriorFaceMarker(Array< int > &face_marker) const
Populate a marker array identifying exterior faces.
Definition mesh.cpp:1555
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6531
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition mesh.hpp:1512
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
Definition mesh.cpp:1476
Geometry::Type GetElementGeometry(int i) const
Definition mesh.hpp:1434
Geometry::Type GetBdrElementGeometry(int i) const
Definition mesh.hpp:1446
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1508
const FiniteElementSpace * GetNodalFESpace() const
Definition mesh.cpp:6479
Geometry::Type GetTypicalElementGeometry() const
If the local mesh is not empty, return GetElementGeometry(0); otherwise, return a typical Geometry pr...
Definition mesh.cpp:1528
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
Definition mesh.hpp:1588
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition mesh.hpp:1291
long GetSequence() const
Definition mesh.hpp:2387
const CoarseFineTransformations & GetRefinementTransforms() const
Definition mesh.cpp:11407
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
const Element * GetFace(int i) const
Return pointer to the i'th face element object.
Definition mesh.hpp:1366
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1354
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
Definition mesh.cpp:7289
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:7514
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition mesh.hpp:1279
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition mesh.cpp:7351
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i....
Definition mesh.cpp:7584
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition mesh.hpp:1526
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1285
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition mesh.hpp:299
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:7267
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:7243
Geometry::Type GetTypicalFaceGeometry() const
If the local mesh is not empty, return GetFaceGeometry(0); otherwise return a typical face geometry p...
Definition mesh.cpp:1496
Geometry::Type GetElementBaseGeometry(int i) const
Definition mesh.hpp:1455
Operator that extracts face degrees of freedom for L2 nonconforming spaces.
const NCList & GetNCList(int entity)
Return vertex/edge/face list (entity = 0/1/2, respectively).
Definition ncmesh.hpp:341
const CoarseFineTransformations & GetDerefinementTransforms() const
Definition ncmesh.cpp:5204
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Definition ncmesh.hpp:319
Geometry::Type GetFaceGeometry(int index) const
Return face geometry type. index is the Mesh face number.
Definition ncmesh.hpp:446
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges, Array< int > &bdr_faces)
Get a list of vertices (2D/3D), edges (3D) and faces (3D) that coincide with boundary elements with t...
Definition ncmesh.cpp:5776
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition ncmesh.hpp:326
int GetNFaces() const
Return the number of (2D) faces in the NCMesh.
Definition ncmesh.hpp:173
int GetNEdges() const
Return the number of edges in the NCMesh.
Definition ncmesh.hpp:171
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:482
DoF transformation implementation for the Nedelec basis on pyramid elements.
Definition doftrans.hpp:373
DoF transformation implementation for the Nedelec basis on tetrahedra.
Definition doftrans.hpp:350
DoF transformation implementation for the Nedelec basis on wedge elements.
Definition doftrans.hpp:361
NURBSExtension generally contains multiple NURBSPatch objects spanning an entire Mesh....
Definition nurbs.hpp:449
const Vector & GetWeights() const
Access function to the vector of weights weights.
Definition nurbs.hpp:845
const Array< int > & GetSlave() const
Definition nurbs.hpp:710
int GetNKV() const
Return the number of KnotVectors.
Definition nurbs.hpp:758
void LoadBE(int i, const FiniteElement *BE) const
Load boundary element i into BE.
Definition nurbs.cpp:4233
Table * GetElementDofTable()
Definition nurbs.hpp:808
const Array< int > & GetMaster() const
Definition nurbs.hpp:708
void LoadFE(int i, const FiniteElement *FE) const
Load element i into FE.
Definition nurbs.cpp:4212
void GetPatchDofs(const int patch, Array< int > &dofs)
Return the degrees of freedom in dofs on patch patch, in Cartesian order.
Definition nurbs.cpp:3935
Table * GetBdrElementDofTable()
Definition nurbs.hpp:812
int GetNDof() const
Return the number of active DOFs.
Definition nurbs.hpp:776
int GetOrder() const
If all KnotVector orders are identical, return that number. Otherwise, return NURBSFECollection::Vari...
Definition nurbs.hpp:755
NURBSExtension * GetCurlExtension(int component)
Definition nurbs.cpp:4422
const Array< int > & GetOrders() const
Read-only access to the orders of all KnotVectors.
Definition nurbs.hpp:751
NURBSExtension * GetDivExtension(int component)
Definition nurbs.cpp:4407
void ConnectBoundaries()
Set DOF maps for periodic BC.
Definition nurbs.cpp:2584
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition fe_coll.hpp:699
int GetOrder() const
Get the order of the NURBS collection: either a positive number, when using fixed order,...
Definition fe_coll.hpp:733
virtual void SetOrder(int Order) const
Set the order and the name, based on the given Order: either a positive number for fixed order,...
Definition fe_coll.cpp:3514
virtual void SetDim(const int dim)
Definition fe_coll.hpp:728
Arbitrary order H(curl) NURBS finite elements.
Definition fe_coll.hpp:810
Arbitrary order H(div) NURBS finite elements.
Definition fe_coll.hpp:758
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
bool OwnsOperator() const
Return true if the OperatorHandle owns the held Operator.
Definition handle.hpp:117
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition handle.hpp:120
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition handle.hpp:132
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged.
Definition handle.hpp:124
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition handle.hpp:145
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition handle.hpp:108
Operator::Type Type() const
Get the currently set operator type id.
Definition handle.hpp:99
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Type
Enumeration defining IDs for some classes derived from Operator.
Definition operator.hpp:284
@ ANY_TYPE
ID for the base class Operator, i.e. any type.
Definition operator.hpp:285
@ MFEM_SPARSEMAT
ID for class SparseMatrix.
Definition operator.hpp:286
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition fespace.hpp:30
Type
Ordering methods:
Definition fespace.hpp:34
A pair of objects.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Parallel version of NURBSExtension.
Definition nurbs.hpp:924
General product operator: x -> (A*B)(x) = A(B(x)).
Definition operator.hpp:894
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
const IntegrationRule * IntRule
Not owned.
const QuadratureSpace * qspace
Not owned.
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:120
Data type sparse matrix.
Definition sparsemat.hpp:51
int GetRow(const int row, Array< int > &cols, Vector &srow) const override
Extract all column indices and values from a given row.
void Add(const int i, const int j, const real_t val)
void Swap(SparseMatrix &other)
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
int RowSize(const int i) const
Returns the number of elements in row i.
real_t * GetRowEntries(const int row)
Return a pointer to the entries in a row.
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
void Finalize(int skip_zeros=1) override
Finalize the matrix initialization, switching the storage format from LIL to CSR.
void Set(const int i, const int j, const real_t val)
int * GetJ()
Definition table.hpp:122
void AddConnections(int r, const int *c, int nc)
Definition table.cpp:176
int RowSize(int i) const
Definition table.hpp:116
void ShiftUpI()
Definition table.cpp:187
void Clear()
Definition table.cpp:453
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition table.cpp:259
void AddConnection(int r, int c)
Definition table.hpp:88
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition table.cpp:153
int Size() const
Returns the number of TYPE I elements.
Definition table.hpp:100
int Size_of_connections() const
Definition table.hpp:106
void AddColumnsInRow(int r, int ncol)
Definition table.hpp:86
void MakeFromList(int nrows, const Array< Connection > &list)
Definition table.cpp:352
void MakeJ()
Definition table.cpp:163
int * GetI()
Definition table.hpp:121
void AddAColumnInRow(int r)
Definition table.hpp:85
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition operator.hpp:847
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition operator.hpp:958
Vector data type.
Definition vector.hpp:82
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition vector.cpp:830
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:679
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:745
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition vector.cpp:126
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:1164
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
real_t Min() const
Returns the minimal element of the vector.
Definition vector.cpp:1123
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:653
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:322
int dim
Definition ex24.cpp:53
int index(int i, int j, int nx, int ny)
Definition life.cpp:236
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void RemoveBasisAndRestriction(const mfem::FiniteElementSpace *fes)
Remove from ceed_basis_map and ceed_restr_map the entries associated with the given fes.
Definition util.cpp:41
Linear1DFiniteElement SegmentFE
Definition segment.cpp:52
void Swap(Array< T > &, Array< T > &)
Definition array.hpp:672
void mfem_error(const char *msg)
Definition error.cpp:154
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:548
void DofMapHelper(int entity, const Table &var_ent_dofs, const Table &loc_var_ent_dofs, const Array< char > &var_ent_orders, const Array< char > &loc_var_ent_orders, Array< int > &all2local, int &ndof_all, int &ndof_loc)
Definition fespace.cpp:2896
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:486
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
Definition text.hpp:45
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition fespace.hpp:1555
BiLinear2DFiniteElement QuadrilateralFE
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
float real_t
Definition config.hpp:43
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
ElementDofOrdering GetEVectorOrdering(const FiniteElementSpace &fes)
Return LEXICOGRAPHIC if mesh contains only one topology and the elements are tensor elements,...
Definition fespace.cpp:4596
void MarkDofs(const Array< int > &dofs, Array< int > &mark_array)
Definition fespace.cpp:576
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition fespace.hpp:83
@ NATIVE
Native ordering as defined by the FiniteElement.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition fe.cpp:32
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
Definition text.hpp:31
FaceType
Definition mesh.hpp:48
real_t p(const Vector &x, real_t t)
bool operator<(const Data &d1, const Data &d2)
Definition nurbs_ex1.cpp:54
Defines the coarse-fine transformations of all fine elements.
Definition ncmesh.hpp:90
Array< Embedding > embeddings
Fine element positions in their parents.
Definition ncmesh.hpp:92
void MakeCoarseToFineTable(Table &coarse_to_fine, bool want_ghosts=false) const
Definition ncmesh.cpp:5262
DenseTensor point_matrices[Geometry::NumGeom]
Definition ncmesh.hpp:96
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition table.hpp:28
Defines the position of a fine element within a coarse element.
Definition ncmesh.hpp:69
int parent
Coarse Element index in the coarse mesh.
Definition ncmesh.hpp:71
unsigned matrix
Definition ncmesh.hpp:78
int index
Mesh number.
Definition ncmesh.hpp:211
Geometry::Type Geom() const
Definition ncmesh.hpp:216
Lists all edges/faces in the nonconforming mesh.
Definition ncmesh.hpp:251
Array< Slave > slaves
All MeshIds corresponding to slave faces.
Definition ncmesh.hpp:254
Array< Master > masters
All MeshIds corresponding to master faces.
Definition ncmesh.hpp:253
Nonconforming edge/face within a bigger edge/face.
Definition ncmesh.hpp:237
unsigned matrix
index into NCList::point_matrices[geom]
Definition ncmesh.hpp:239