MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fespace.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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"
16 #include "../mesh/mesh_headers.hpp"
17 #include "fem.hpp"
18 #include "ceed/util.hpp"
19 
20 #include <cmath>
21 #include <cstdarg>
22 
23 using namespace std;
24 
25 namespace mfem
26 {
27 
28 template <> void Ordering::
29 DofsToVDofs<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 
43 template <> void Ordering::
44 DofsToVDofs<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 
59 FiniteElementSpace::FiniteElementSpace()
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  DoFTrans(0), VDoFTrans(vdim, ordering),
67  cP(NULL), cR(NULL), cR_hp(NULL), cP_is_set(false),
68  Th(Operator::ANY_TYPE),
69  sequence(0), mesh_sequence(0), orders_changed(false), relaxed_hp(false)
70 { }
71 
73  Mesh *mesh_,
74  const FiniteElementCollection *fec_)
75  : VDoFTrans(orig.vdim, orig.ordering)
76 {
77  mesh_ = mesh_ ? mesh_ : orig.mesh;
78  fec_ = fec_ ? fec_ : orig.fec;
79 
80  NURBSExtension *nurbs_ext = NULL;
81  if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
82  {
83 #ifdef MFEM_USE_MPI
84  ParNURBSExtension *pNURBSext =
85  dynamic_cast<ParNURBSExtension *>(orig.NURBSext);
86  if (pNURBSext)
87  {
88  nurbs_ext = new ParNURBSExtension(*pNURBSext);
89  }
90  else
91 #endif
92  {
93  nurbs_ext = new NURBSExtension(*orig.NURBSext);
94  }
95  }
96 
97  Constructor(mesh_, nurbs_ext, fec_, orig.vdim, orig.ordering);
98 }
99 
101  const FiniteElementSpace &fes, const Array<int> *perm)
102 {
103  MFEM_VERIFY(cP == NULL, "");
104  MFEM_VERIFY(cR == NULL, "");
105 
106  SparseMatrix *perm_mat = NULL, *perm_mat_tr = NULL;
107  if (perm)
108  {
109  // Note: although n and fes.GetVSize() are typically equal, in
110  // variable-order spaces they may differ, since nonconforming edges/faces
111  // my have fictitious DOFs.
112  int n = perm->Size();
113  perm_mat = new SparseMatrix(n, fes.GetVSize());
114  for (int i=0; i<n; ++i)
115  {
116  double s;
117  int j = DecodeDof((*perm)[i], s);
118  perm_mat->Set(i, j, s);
119  }
120  perm_mat->Finalize();
121  perm_mat_tr = Transpose(*perm_mat);
122  }
123 
124  if (fes.GetConformingProlongation() != NULL)
125  {
126  if (perm) { cP = Mult(*perm_mat, *fes.GetConformingProlongation()); }
127  else { cP = new SparseMatrix(*fes.GetConformingProlongation()); }
128  cP_is_set = true;
129  }
130  else if (perm != NULL)
131  {
132  cP = perm_mat;
133  cP_is_set = true;
134  perm_mat = NULL;
135  }
136  if (fes.GetConformingRestriction() != NULL)
137  {
138  if (perm) { cR = Mult(*fes.GetConformingRestriction(), *perm_mat_tr); }
139  else { cR = new SparseMatrix(*fes.GetConformingRestriction()); }
140  }
141  else if (perm != NULL)
142  {
143  cR = perm_mat_tr;
144  perm_mat_tr = NULL;
145  }
146 
147  delete perm_mat;
148  delete perm_mat_tr;
149 }
150 
152 {
153  MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),
154  "Space has not been Updated() after a Mesh change.");
155  MFEM_VERIFY(i >= 0 && i < GetNE(), "Invalid element index");
156  MFEM_VERIFY(p >= 0 && p <= MaxVarOrder, "Order out of range");
157  MFEM_ASSERT(!elem_order.Size() || elem_order.Size() == GetNE(),
158  "Internal error");
159 
160  if (elem_order.Size()) // already a variable-order space
161  {
162  if (elem_order[i] != p)
163  {
164  elem_order[i] = p;
165  orders_changed = true;
166  }
167  }
168  else // convert space to variable-order space
169  {
171  elem_order = fec->GetOrder();
172 
173  elem_order[i] = p;
174  orders_changed = true;
175  }
176 }
177 
179 {
180  MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),
181  "Space has not been Updated() after a Mesh change.");
182  MFEM_VERIFY(i >= 0 && i < GetNE(), "Invalid element index");
183  MFEM_ASSERT(!elem_order.Size() || elem_order.Size() == GetNE(),
184  "Internal error");
185 
186  return GetElementOrderImpl(i);
187 }
188 
190 {
191  // (this is an internal version of GetElementOrder without asserts and checks)
192  return elem_order.Size() ? elem_order[i] : fec->GetOrder();
193 }
194 
195 void FiniteElementSpace::GetVDofs(int vd, Array<int>& dofs, int ndofs_) const
196 {
197  if (ndofs_ < 0) { ndofs_ = this->ndofs; }
198 
200  {
201  for (int i = 0; i < dofs.Size(); i++)
202  {
203  dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs_, vdim, i, vd);
204  }
205  }
206  else
207  {
208  for (int i = 0; i < dofs.Size(); i++)
209  {
210  dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs_, vdim, i, vd);
211  }
212  }
213 }
214 
215 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs, int ndofs_) const
216 {
217  if (vdim == 1) { return; }
218  if (ndofs_ < 0) { ndofs_ = this->ndofs; }
219 
221  {
222  Ordering::DofsToVDofs<Ordering::byNODES>(ndofs_, vdim, dofs);
223  }
224  else
225  {
226  Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs_, vdim, dofs);
227  }
228 }
229 
230 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs, int ndofs_) const
231 {
232  if (vdim == 1) { return; }
233  if (ndofs_ < 0) { ndofs_ = this->ndofs; }
234 
236  {
237  for (int i = 0; i < dofs.Size(); i++)
238  {
239  dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs_, vdim, dofs[i], vd);
240  }
241  }
242  else
243  {
244  for (int i = 0; i < dofs.Size(); i++)
245  {
246  dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs_, vdim, dofs[i], vd);
247  }
248  }
249 }
250 
251 int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
252 {
253  if (vdim == 1) { return dof; }
254  if (ndofs_ < 0) { ndofs_ = this->ndofs; }
255 
257  {
258  return Ordering::Map<Ordering::byNODES>(ndofs_, vdim, dof, vd);
259  }
260  else
261  {
262  return Ordering::Map<Ordering::byVDIM>(ndofs_, vdim, dof, vd);
263  }
264 }
265 
266 // static function
268 {
269  int n = vdofs.Size(), *vdof = vdofs;
270  for (int i = 0; i < n; i++)
271  {
272  int j;
273  if ((j = vdof[i]) < 0)
274  {
275  vdof[i] = -1-j;
276  }
277  }
278 }
279 
282 {
283  DofTransformation * doftrans = GetElementDofs(i, vdofs);
284  DofsToVDofs(vdofs);
285  if (vdim == 1 || doftrans == NULL)
286  {
287  return doftrans;
288  }
289  else
290  {
291  VDoFTrans.SetDofTransformation(*doftrans);
292  return &VDoFTrans;
293  }
294 }
295 
298 {
299  DofTransformation * doftrans = GetBdrElementDofs(i, vdofs);
300  DofsToVDofs(vdofs);
301  if (vdim == 1 || doftrans == NULL)
302  {
303  return doftrans;
304  }
305  else
306  {
307  VDoFTrans.SetDofTransformation(*doftrans);
308  return &VDoFTrans;
309  }
310 }
311 
313 {
314  GetFaceDofs(i, vdofs);
315  DofsToVDofs(vdofs);
316 }
317 
319 {
320  GetEdgeDofs(i, vdofs);
321  DofsToVDofs(vdofs);
322 }
323 
325 {
326  GetVertexDofs(i, vdofs);
327  DofsToVDofs(vdofs);
328 }
329 
331 {
332  GetElementInteriorDofs(i, vdofs);
333  DofsToVDofs(vdofs);
334 }
335 
337 {
338  GetEdgeInteriorDofs(i, vdofs);
339  DofsToVDofs(vdofs);
340 }
341 
343 {
344  if (elem_dof) { return; }
345 
346  // TODO: can we call GetElementDofs only once per element?
347  Table *el_dof = new Table;
348  Table *el_fos = (mesh->Dimension() > 2) ? (new Table) : NULL;
349  Array<int> dofs;
350  Array<int> F, Fo;
351  el_dof -> MakeI (mesh -> GetNE());
352  if (el_fos) { el_fos -> MakeI (mesh -> GetNE()); }
353  for (int i = 0; i < mesh -> GetNE(); i++)
354  {
355  GetElementDofs (i, dofs);
356  el_dof -> AddColumnsInRow (i, dofs.Size());
357 
358  if (el_fos)
359  {
360  mesh->GetElementFaces(i, F, Fo);
361  el_fos -> AddColumnsInRow (i, Fo.Size());
362  }
363  }
364  el_dof -> MakeJ();
365  if (el_fos) { el_fos -> MakeJ(); }
366  for (int i = 0; i < mesh -> GetNE(); i++)
367  {
368  GetElementDofs (i, dofs);
369  el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
370 
371  if (el_fos)
372  {
373  mesh->GetElementFaces(i, F, Fo);
374  el_fos -> AddConnections (i, (int *)Fo, Fo.Size());
375  }
376  }
377  el_dof -> ShiftUpI();
378  if (el_fos) { el_fos -> ShiftUpI(); }
379  elem_dof = el_dof;
380  elem_fos = el_fos;
381 }
382 
384 {
385  if (bdr_elem_dof) { return; }
386 
387  Table *bel_dof = new Table;
388  Array<int> dofs;
389  bel_dof->MakeI(mesh->GetNBE());
390  for (int i = 0; i < mesh->GetNBE(); i++)
391  {
392  GetBdrElementDofs(i, dofs);
393  bel_dof->AddColumnsInRow(i, dofs.Size());
394  }
395  bel_dof->MakeJ();
396  for (int i = 0; i < mesh->GetNBE(); i++)
397  {
398  GetBdrElementDofs(i, dofs);
399  bel_dof->AddConnections(i, (int *)dofs, dofs.Size());
400  }
401  bel_dof->ShiftUpI();
402  bdr_elem_dof = bel_dof;
403 }
404 
406 {
407  // Here, "face" == (dim-1)-dimensional mesh entity.
408 
409  if (face_dof) { return; }
410 
411  if (NURBSext) { BuildNURBSFaceToDofTable(); return; }
412 
413  Table *fc_dof = new Table;
414  Array<int> dofs;
415  fc_dof->MakeI(mesh->GetNumFaces());
416  for (int i = 0; i < fc_dof->Size(); i++)
417  {
418  GetFaceDofs(i, dofs, 0);
419  fc_dof->AddColumnsInRow(i, dofs.Size());
420  }
421  fc_dof->MakeJ();
422  for (int i = 0; i < fc_dof->Size(); i++)
423  {
424  GetFaceDofs(i, dofs, 0);
425  fc_dof->AddConnections(i, (int *)dofs, dofs.Size());
426  }
427  fc_dof->ShiftUpI();
428  face_dof = fc_dof;
429 }
430 
432 {
433  delete elem_dof;
434  delete elem_fos;
435  elem_dof = NULL;
436  elem_fos = NULL;
438 }
439 
441 {
442  Array<int> dof_marker(ndofs);
443 
444  dof_marker = -1;
445 
446  int *J = elem_dof->GetJ(), nnz = elem_dof->Size_of_connections();
447  for (int k = 0, dof_counter = 0; k < nnz; k++)
448  {
449  const int sdof = J[k]; // signed dof
450  const int dof = (sdof < 0) ? -1-sdof : sdof;
451  int new_dof = dof_marker[dof];
452  if (new_dof < 0)
453  {
454  dof_marker[dof] = new_dof = dof_counter++;
455  }
456  J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
457  }
458 }
459 
461 {
462  if (dof_elem_array.Size()) { return; }
463 
465 
468  dof_elem_array = -1;
469  for (int i = 0; i < mesh -> GetNE(); i++)
470  {
471  const int *dofs = elem_dof -> GetRow(i);
472  const int n = elem_dof -> RowSize(i);
473  for (int j = 0; j < n; j++)
474  {
475  int dof = DecodeDof(dofs[j]);
476  if (dof_elem_array[dof] < 0)
477  {
478  dof_elem_array[dof] = i;
479  dof_ldof_array[dof] = j;
480  }
481  }
482  }
483 }
484 
485 static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
486 {
487  for (int i = 0; i < dofs.Size(); i++)
488  {
489  int k = dofs[i];
490  if (k < 0) { k = -1 - k; }
491  mark_array[k] = -1;
492  }
493 }
494 
496  Array<int> &ess_vdofs,
497  int component) const
498 {
499  Array<int> vdofs, dofs;
500 
501  ess_vdofs.SetSize(GetVSize());
502  ess_vdofs = 0;
503 
504  for (int i = 0; i < GetNBE(); i++)
505  {
506  if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
507  {
508  if (component < 0)
509  {
510  // Mark all components.
511  GetBdrElementVDofs(i, vdofs);
512  mark_dofs(vdofs, ess_vdofs);
513  }
514  else
515  {
516  GetBdrElementDofs(i, dofs);
517  for (int d = 0; d < dofs.Size(); d++)
518  { dofs[d] = DofToVDof(dofs[d], component); }
519  mark_dofs(dofs, ess_vdofs);
520  }
521  }
522  }
523 
524  // mark possible hidden boundary edges in a non-conforming mesh, also
525  // local DOFs affected by boundary elements on other processors
526  if (Nonconforming())
527  {
528  Array<int> bdr_verts, bdr_edges;
529  mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges);
530 
531  for (int i = 0; i < bdr_verts.Size(); i++)
532  {
533  if (component < 0)
534  {
535  GetVertexVDofs(bdr_verts[i], vdofs);
536  mark_dofs(vdofs, ess_vdofs);
537  }
538  else
539  {
540  GetVertexDofs(bdr_verts[i], dofs);
541  for (int d = 0; d < dofs.Size(); d++)
542  { dofs[d] = DofToVDof(dofs[d], component); }
543  mark_dofs(dofs, ess_vdofs);
544  }
545  }
546  for (int i = 0; i < bdr_edges.Size(); i++)
547  {
548  if (component < 0)
549  {
550  GetEdgeVDofs(bdr_edges[i], vdofs);
551  mark_dofs(vdofs, ess_vdofs);
552  }
553  else
554  {
555  GetEdgeDofs(bdr_edges[i], dofs);
556  for (int d = 0; d < dofs.Size(); d++)
557  { dofs[d] = DofToVDof(dofs[d], component); }
558  mark_dofs(dofs, ess_vdofs);
559  }
560  }
561  }
562 }
563 
565  Array<int> &ess_tdof_list,
566  int component)
567 {
568  Array<int> ess_vdofs, ess_tdofs;
569  GetEssentialVDofs(bdr_attr_is_ess, ess_vdofs, component);
571  if (!R)
572  {
573  ess_tdofs.MakeRef(ess_vdofs);
574  }
575  else
576  {
577  R->BooleanMult(ess_vdofs, ess_tdofs);
578  }
579  MarkerToList(ess_tdofs, ess_tdof_list);
580 }
581 
583  int component)
584 {
585  if (mesh->bdr_attributes.Size())
586  {
587  Array<int> ess_bdr(mesh->bdr_attributes.Max());
588  ess_bdr = 1;
589  GetEssentialTrueDofs(ess_bdr, boundary_dofs, component);
590  }
591  else
592  {
593  boundary_dofs.DeleteAll();
594  }
595 }
596 
597 // static method
599  Array<int> &list)
600 {
601  int num_marked = 0;
602  marker.HostRead(); // make sure we can read the array on host
603  for (int i = 0; i < marker.Size(); i++)
604  {
605  if (marker[i]) { num_marked++; }
606  }
607  list.SetSize(0);
608  list.HostWrite();
609  list.Reserve(num_marked);
610  for (int i = 0; i < marker.Size(); i++)
611  {
612  if (marker[i]) { list.Append(i); }
613  }
614 }
615 
616 // static method
617 void FiniteElementSpace::ListToMarker(const Array<int> &list, int marker_size,
618  Array<int> &marker, int mark_val)
619 {
620  list.HostRead(); // make sure we can read the array on host
621  marker.SetSize(marker_size);
622  marker.HostWrite();
623  marker = 0;
624  for (int i = 0; i < list.Size(); i++)
625  {
626  marker[list[i]] = mark_val;
627  }
628 }
629 
631  Array<int> &cdofs)
632 {
634  if (cP) { cP->BooleanMultTranspose(dofs, cdofs); }
635  else { dofs.Copy(cdofs); }
636 }
637 
639  Array<int> &dofs)
640 {
642  if (cR) { cR->BooleanMultTranspose(cdofs, dofs); }
643  else { cdofs.Copy(dofs); }
644 }
645 
646 SparseMatrix *
648 {
649  int i, j;
650  Array<int> d_vdofs, c_vdofs;
651  SparseMatrix *R;
652 
653  R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
654 
655  for (i = 0; i < mesh -> GetNE(); i++)
656  {
657  this -> GetElementVDofs (i, d_vdofs);
658  cfes -> GetElementVDofs (i, c_vdofs);
659 
660 #ifdef MFEM_DEBUG
661  if (d_vdofs.Size() != c_vdofs.Size())
662  {
663  mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
664  }
665 #endif
666 
667  for (j = 0; j < d_vdofs.Size(); j++)
668  {
669  R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
670  }
671  }
672 
673  R -> Finalize();
674 
675  return R;
676 }
677 
678 SparseMatrix *
680 {
681  int i, j;
682  Array<int> d_dofs, c_dofs;
683  SparseMatrix *R;
684 
685  R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
686 
687  for (i = 0; i < mesh -> GetNE(); i++)
688  {
689  this -> GetElementDofs (i, d_dofs);
690  cfes -> GetElementDofs (i, c_dofs);
691 
692 #ifdef MFEM_DEBUG
693  if (c_dofs.Size() != 1)
694  mfem_error ("FiniteElementSpace::"
695  "D2Const_GlobalRestrictionMatrix (...)");
696 #endif
697 
698  for (j = 0; j < d_dofs.Size(); j++)
699  {
700  R -> Set (c_dofs[0], d_dofs[j], 1.0);
701  }
702  }
703 
704  R -> Finalize();
705 
706  return R;
707 }
708 
709 SparseMatrix *
711 {
712  SparseMatrix *R;
713  DenseMatrix loc_restr;
714  Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
715 
716  int lvdim = lfes->GetVDim();
717  R = new SparseMatrix (lvdim * lfes -> GetNDofs(), lvdim * ndofs);
718 
719  Geometry::Type cached_geom = Geometry::INVALID;
720  const FiniteElement *h_fe = NULL;
721  const FiniteElement *l_fe = NULL;
723 
724  for (int i = 0; i < mesh -> GetNE(); i++)
725  {
726  this -> GetElementDofs (i, h_dofs);
727  lfes -> GetElementDofs (i, l_dofs);
728 
729  // Assuming 'loc_restr' depends only on the Geometry::Type.
730  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
731  if (geom != cached_geom)
732  {
733  h_fe = this -> GetFE (i);
734  l_fe = lfes -> GetFE (i);
736  h_fe->Project(*l_fe, T, loc_restr);
737  cached_geom = geom;
738  }
739 
740  for (int vd = 0; vd < lvdim; vd++)
741  {
742  l_dofs.Copy(l_vdofs);
743  lfes->DofsToVDofs(vd, l_vdofs);
744 
745  h_dofs.Copy(h_vdofs);
746  this->DofsToVDofs(vd, h_vdofs);
747 
748  R -> SetSubMatrix (l_vdofs, h_vdofs, loc_restr, 1);
749  }
750  }
751 
752  R -> Finalize();
753 
754  return R;
755 }
756 
759  Array<int>& slave_dofs, DenseMatrix& I, int skipfirst)
760 {
761  for (int i = skipfirst; i < slave_dofs.Size(); i++)
762  {
763  int sdof = slave_dofs[i];
764  if (!deps.RowSize(sdof)) // not processed yet
765  {
766  for (int j = 0; j < master_dofs.Size(); j++)
767  {
768  double coef = I(i, j);
769  if (std::abs(coef) > 1e-12)
770  {
771  int mdof = master_dofs[j];
772  if (mdof != sdof && mdof != (-1-sdof))
773  {
774  deps.Add(sdof, mdof, coef);
775  }
776  }
777  }
778  }
779  }
780 }
781 
784  const FiniteElement *master_fe,
785  Array<int> &slave_dofs, int slave_face,
786  const DenseMatrix *pm) const
787 {
788  // In variable-order spaces in 3D, we need to only constrain interior face
789  // DOFs (this is done one level up), since edge dependencies can be more
790  // complex and are primarily handled by edge-edge dependencies. The one
791  // exception is edges of slave faces that lie in the interior of the master
792  // face, which are not covered by edge-edge relations. This function finds
793  // such edges and makes them constrained by the master face.
794  // See also https://github.com/mfem/mfem/pull/1423#issuecomment-633916643
795 
796  Array<int> V, E, Eo; // TODO: LocalArray
797  mesh->GetFaceVertices(slave_face, V);
798  mesh->GetFaceEdges(slave_face, E, Eo);
799  MFEM_ASSERT(V.Size() == E.Size(), "");
800 
801  DenseMatrix I;
803  edge_T.SetFE(&SegmentFE);
804 
805  // constrain each edge of the slave face
806  for (int i = 0; i < E.Size(); i++)
807  {
808  int a = i, b = (i+1) % V.Size();
809  if (V[a] > V[b]) { std::swap(a, b); }
810 
811  DenseMatrix &edge_pm = edge_T.GetPointMat();
812  edge_pm.SetSize(2, 2);
813 
814  // copy two points from the face point matrix
815  double mid[2];
816  for (int j = 0; j < 2; j++)
817  {
818  edge_pm(j, 0) = (*pm)(j, a);
819  edge_pm(j, 1) = (*pm)(j, b);
820  mid[j] = 0.5*((*pm)(j, a) + (*pm)(j, b));
821  }
822 
823  // check that the edge does not coincide with the master face's edge
824  const double eps = 1e-14;
825  if (mid[0] > eps && mid[0] < 1-eps &&
826  mid[1] > eps && mid[1] < 1-eps)
827  {
828  int order = GetEdgeDofs(E[i], slave_dofs, 0);
829 
830  const auto *edge_fe = fec->GetFE(Geometry::SEGMENT, order);
831  edge_fe->GetTransferMatrix(*master_fe, edge_T, I);
832 
833  AddDependencies(deps, master_dofs, slave_dofs, I, 0);
834  }
835  }
836 }
837 
838 bool FiniteElementSpace::DofFinalizable(int dof, const Array<bool>& finalized,
839  const SparseMatrix& deps)
840 {
841  const int* dep = deps.GetRowColumns(dof);
842  int ndep = deps.RowSize(dof);
843 
844  // are all constraining DOFs finalized?
845  for (int i = 0; i < ndep; i++)
846  {
847  if (!finalized[dep[i]]) { return false; }
848  }
849  return true;
850 }
851 
853  Geometry::Type master_geom,
854  int variant) const
855 {
856  // In NC meshes with prisms/tets, a special constraint occurs where a
857  // prism/tet edge is slave to another element's face (see illustration
858  // here: https://github.com/mfem/mfem/pull/713#issuecomment-495786362)
859  // Rather than introduce a new edge-face constraint type, we handle such
860  // cases as degenerate face-face constraints, where the point-matrix
861  // rectangle has zero height. This method returns DOFs for the first edge
862  // of the rectangle, duplicated in the orthogonal direction, to resemble
863  // DOFs for a quadrilateral face. The extra DOFs are ignored by
864  // FiniteElementSpace::AddDependencies.
865 
866  Array<int> edof;
867  int order = GetEdgeDofs(-1 - index, edof, variant);
868 
871  int nn = 2*nv + ne;
872 
873  dofs.SetSize(nn*nn);
874  if (!dofs.Size()) { return 0; }
875 
876  dofs = edof[0];
877 
878  // copy first two vertex DOFs
879  for (int i = 0; i < nv; i++)
880  {
881  dofs[i] = edof[i];
882  dofs[nv+i] = edof[nv+i];
883  }
884  // copy first edge DOFs
885  int face_vert = Geometry::NumVerts[master_geom];
886  for (int i = 0; i < ne; i++)
887  {
888  dofs[face_vert*nv + i] = edof[2*nv + i];
889  }
890 
891  return order;
892 }
893 
895 {
896  // return the number of vertex and edge DOFs that precede inner DOFs
897  const int nv = fec->GetNumDof(Geometry::POINT, order);
898  const int ne = fec->GetNumDof(Geometry::SEGMENT, order);
899 
900  return Geometry::NumVerts[geom] * (geom == Geometry::SEGMENT ? nv : (nv + ne));
901 }
902 
904  Geometry::Type master_geom,
905  int variant) const
906 {
907  switch (entity)
908  {
909  case 0:
910  GetVertexDofs(index, dofs);
911  return 0;
912 
913  case 1:
914  return GetEdgeDofs(index, dofs, variant);
915 
916  default:
917  if (index >= 0)
918  {
919  return GetFaceDofs(index, dofs, variant);
920  }
921  else
922  {
923  return GetDegenerateFaceDofs(index, dofs, master_geom, variant);
924  }
925  }
926 }
927 
929 {
930 #ifdef MFEM_USE_MPI
931  MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
932  "This method should not be used with a ParFiniteElementSpace!");
933 #endif
934 
935  if (cP_is_set) { return; }
936  cP_is_set = true;
937 
938  if (FEColl()->GetContType() == FiniteElementCollection::DISCONTINUOUS)
939  {
940  cP = cR = cR_hp = NULL; // will be treated as identities
941  return;
942  }
943 
944  Array<int> master_dofs, slave_dofs, highest_dofs;
945 
947  DenseMatrix I;
948 
949  // For each slave DOF, the dependency matrix will contain a row that
950  // expresses the slave DOF as a linear combination of its immediate master
951  // DOFs. Rows of independent DOFs will remain empty.
952  SparseMatrix deps(ndofs);
953 
954  // Inverse dependencies for the cR_hp matrix in variable order spaces:
955  // For each master edge/face with more DOF sets, the inverse dependency
956  // matrix contains a row that expresses the master true DOF (lowest order)
957  // as a linear combination of the highest order set of DOFs.
958  SparseMatrix inv_deps(ndofs);
959 
960  // collect local face/edge dependencies
961  for (int entity = 2; entity >= 1; entity--)
962  {
963  const NCMesh::NCList &list = mesh->ncmesh->GetNCList(entity);
964  if (!list.masters.Size()) { continue; }
965 
966  // loop through all master edges/faces, constrain their slave edges/faces
967  for (const NCMesh::Master &master : list.masters)
968  {
969  Geometry::Type master_geom = master.Geom();
970 
971  int p = GetEntityDofs(entity, master.index, master_dofs, master_geom);
972  if (!master_dofs.Size()) { continue; }
973 
974  const FiniteElement *master_fe = fec->GetFE(master_geom, p);
975  if (!master_fe) { continue; }
976 
977  switch (master_geom)
978  {
979  case Geometry::SQUARE: T.SetFE(&QuadrilateralFE); break;
980  case Geometry::TRIANGLE: T.SetFE(&TriangleFE); break;
981  case Geometry::SEGMENT: T.SetFE(&SegmentFE); break;
982  default: MFEM_ABORT("unsupported geometry");
983  }
984 
985  for (int si = master.slaves_begin; si < master.slaves_end; si++)
986  {
987  const NCMesh::Slave &slave = list.slaves[si];
988 
989  int q = GetEntityDofs(entity, slave.index, slave_dofs, master_geom);
990  if (!slave_dofs.Size()) { break; }
991 
992  const FiniteElement *slave_fe = fec->GetFE(slave.Geom(), q);
993  list.OrientedPointMatrix(slave, T.GetPointMat());
994  slave_fe->GetTransferMatrix(*master_fe, T, I);
995 
996  // variable order spaces: face edges need to be handled separately
997  int skipfirst = 0;
998  if (IsVariableOrder() && entity == 2 && slave.index >= 0)
999  {
1000  skipfirst = GetNumBorderDofs(master_geom, q);
1001  }
1002 
1003  // make each slave DOF dependent on all master DOFs
1004  AddDependencies(deps, master_dofs, slave_dofs, I, skipfirst);
1005 
1006  if (skipfirst)
1007  {
1008  // constrain internal edge DOFs if they were skipped
1009  const auto *pm = list.point_matrices[master_geom][slave.matrix];
1010  AddEdgeFaceDependencies(deps, master_dofs, master_fe,
1011  slave_dofs, slave.index, pm);
1012  }
1013  }
1014 
1015  // Add inverse dependencies for the cR_hp matrix; if a master has
1016  // more DOF sets, the lowest order set interpolates the highest one.
1017  if (IsVariableOrder())
1018  {
1019  int nvar = GetNVariants(entity, master.index);
1020  if (nvar > 1)
1021  {
1022  int q = GetEntityDofs(entity, master.index, highest_dofs,
1023  master_geom, nvar-1);
1024  const auto *highest_fe = fec->GetFE(master_geom, q);
1025 
1026  T.SetIdentityTransformation(master_geom);
1027  master_fe->GetTransferMatrix(*highest_fe, T, I);
1028 
1029  // add dependencies only for the inner dofs
1030  int skip = GetNumBorderDofs(master_geom, p);
1031  AddDependencies(inv_deps, highest_dofs, master_dofs, I, skip);
1032  }
1033  }
1034  }
1035  }
1036 
1037  // variable order spaces: enforce minimum rule on conforming edges/faces
1038  if (IsVariableOrder())
1039  {
1040  for (int entity = 1; entity < mesh->Dimension(); entity++)
1041  {
1042  const Table &ent_dofs = (entity == 1) ? var_edge_dofs : var_face_dofs;
1043  int num_ent = (entity == 1) ? mesh->GetNEdges() : mesh->GetNFaces();
1044  MFEM_ASSERT(ent_dofs.Size() == num_ent+1, "");
1045 
1046  // add constraints within edges/faces holding multiple DOF sets
1047  Geometry::Type last_geom = Geometry::INVALID;
1048  for (int i = 0; i < num_ent; i++)
1049  {
1050  if (ent_dofs.RowSize(i) <= 1) { continue; }
1051 
1052  Geometry::Type geom =
1053  (entity == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
1054 
1055  if (geom != last_geom)
1056  {
1057  T.SetIdentityTransformation(geom);
1058  last_geom = geom;
1059  }
1060 
1061  // get lowest order variant DOFs and FE
1062  int p = GetEntityDofs(entity, i, master_dofs, geom, 0);
1063  const auto *master_fe = fec->GetFE(geom, p);
1064  if (!master_fe) { break; }
1065 
1066  // constrain all higher order DOFs: interpolate lowest order function
1067  for (int variant = 1; ; variant++)
1068  {
1069  int q = GetEntityDofs(entity, i, slave_dofs, geom, variant);
1070  if (q < 0) { break; }
1071 
1072  const auto *slave_fe = fec->GetFE(geom, q);
1073  slave_fe->GetTransferMatrix(*master_fe, T, I);
1074 
1075  AddDependencies(deps, master_dofs, slave_dofs, I);
1076  }
1077  }
1078  }
1079  }
1080 
1081  deps.Finalize();
1082  inv_deps.Finalize();
1083 
1084  // DOFs that stayed independent are true DOFs
1085  int n_true_dofs = 0;
1086  for (int i = 0; i < ndofs; i++)
1087  {
1088  if (!deps.RowSize(i)) { n_true_dofs++; }
1089  }
1090 
1091  // if all dofs are true dofs leave cP and cR NULL
1092  if (n_true_dofs == ndofs)
1093  {
1094  cP = cR = cR_hp = NULL; // will be treated as identities
1095  return;
1096  }
1097 
1098  // create the conforming prolongation matrix cP
1099  cP = new SparseMatrix(ndofs, n_true_dofs);
1100 
1101  // create the conforming restriction matrix cR
1102  int *cR_J;
1103  {
1104  int *cR_I = Memory<int>(n_true_dofs+1);
1105  double *cR_A = Memory<double>(n_true_dofs);
1106  cR_J = Memory<int>(n_true_dofs);
1107  for (int i = 0; i < n_true_dofs; i++)
1108  {
1109  cR_I[i] = i;
1110  cR_A[i] = 1.0;
1111  }
1112  cR_I[n_true_dofs] = n_true_dofs;
1113  cR = new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, ndofs);
1114  }
1115 
1116  // In var. order spaces, create the restriction matrix cR_hp which is similar
1117  // to cR, but has interpolation in the extra master edge/face DOFs.
1118  cR_hp = IsVariableOrder() ? new SparseMatrix(n_true_dofs, ndofs) : NULL;
1119 
1120  Array<bool> finalized(ndofs);
1121  finalized = false;
1122 
1123  Array<int> cols;
1124  Vector srow;
1125 
1126  // put identity in the prolongation matrix for true DOFs, initialize cR_hp
1127  for (int i = 0, true_dof = 0; i < ndofs; i++)
1128  {
1129  if (!deps.RowSize(i)) // true dof
1130  {
1131  cP->Add(i, true_dof, 1.0);
1132  cR_J[true_dof] = i;
1133  finalized[i] = true;
1134 
1135  if (cR_hp)
1136  {
1137  if (inv_deps.RowSize(i))
1138  {
1139  inv_deps.GetRow(i, cols, srow);
1140  cR_hp->AddRow(true_dof, cols, srow);
1141  }
1142  else
1143  {
1144  cR_hp->Add(true_dof, i, 1.0);
1145  }
1146  }
1147 
1148  true_dof++;
1149  }
1150  }
1151 
1152  // Now calculate cP rows of slave DOFs as combinations of cP rows of their
1153  // master DOFs. It is possible that some slave DOFs depend on DOFs that are
1154  // themselves slaves. Here we resolve such indirect constraints by first
1155  // calculating rows of the cP matrix for DOFs whose master DOF cP rows are
1156  // already known (in the first iteration these are the true DOFs). In the
1157  // second iteration, slaves of slaves can be 'finalized' (given a row in the
1158  // cP matrix), in the third iteration slaves of slaves of slaves, etc.
1159  bool finished;
1160  int n_finalized = n_true_dofs;
1161  do
1162  {
1163  finished = true;
1164  for (int dof = 0; dof < ndofs; dof++)
1165  {
1166  if (!finalized[dof] && DofFinalizable(dof, finalized, deps))
1167  {
1168  const int* dep_col = deps.GetRowColumns(dof);
1169  const double* dep_coef = deps.GetRowEntries(dof);
1170  int n_dep = deps.RowSize(dof);
1171 
1172  for (int j = 0; j < n_dep; j++)
1173  {
1174  cP->GetRow(dep_col[j], cols, srow);
1175  srow *= dep_coef[j];
1176  cP->AddRow(dof, cols, srow);
1177  }
1178 
1179  finalized[dof] = true;
1180  n_finalized++;
1181  finished = false;
1182  }
1183  }
1184  }
1185  while (!finished);
1186 
1187  // If everything is consistent (mesh, face orientations, etc.), we should
1188  // be able to finalize all slave DOFs, otherwise it's a serious error.
1189  MFEM_VERIFY(n_finalized == ndofs,
1190  "Error creating cP matrix: n_finalized = "
1191  << n_finalized << ", ndofs = " << ndofs);
1192 
1193  cP->Finalize();
1194  if (cR_hp) { cR_hp->Finalize(); }
1195 
1196  if (vdim > 1)
1197  {
1198  MakeVDimMatrix(*cP);
1199  MakeVDimMatrix(*cR);
1200  if (cR_hp) { MakeVDimMatrix(*cR_hp); }
1201  }
1202 
1204 }
1205 
1207 {
1208  if (vdim == 1) { return; }
1209 
1210  int height = mat.Height();
1211  int width = mat.Width();
1212 
1213  SparseMatrix *vmat = new SparseMatrix(vdim*height, vdim*width);
1214 
1215  Array<int> dofs, vdofs;
1216  Vector srow;
1217  for (int i = 0; i < height; i++)
1218  {
1219  mat.GetRow(i, dofs, srow);
1220  for (int vd = 0; vd < vdim; vd++)
1221  {
1222  dofs.Copy(vdofs);
1223  DofsToVDofs(vd, vdofs, width);
1224  vmat->SetRow(DofToVDof(i, vd, height), vdofs, srow);
1225  }
1226  }
1227  vmat->Finalize();
1228 
1229  mat.Swap(*vmat);
1230  delete vmat;
1231 }
1232 
1233 
1235 {
1236  if (Conforming()) { return NULL; }
1238  return cP;
1239 }
1240 
1242 {
1243  if (Conforming()) { return NULL; }
1245  return cR;
1246 }
1247 
1249 {
1250  if (Conforming()) { return NULL; }
1252  return IsVariableOrder() ? cR_hp : cR;
1253 }
1254 
1256 {
1258  return P ? (P->Width() / vdim) : ndofs;
1259 }
1260 
1262  ElementDofOrdering e_ordering) const
1263 {
1264  // Check if we have a discontinuous space using the FE collection:
1265  if (IsDGSpace())
1266  {
1267  // TODO: when VDIM is 1, we can return IdentityOperator.
1268  if (L2E_nat.Ptr() == NULL)
1269  {
1270  // The input L-vector layout is:
1271  // * ND x NE x VDIM, for Ordering::byNODES, or
1272  // * VDIM x ND x NE, for Ordering::byVDIM.
1273  // The output E-vector layout is: ND x VDIM x NE.
1274  L2E_nat.Reset(new L2ElementRestriction(*this));
1275  }
1276  return L2E_nat.Ptr();
1277  }
1278  if (e_ordering == ElementDofOrdering::LEXICOGRAPHIC)
1279  {
1280  if (L2E_lex.Ptr() == NULL)
1281  {
1282  L2E_lex.Reset(new ElementRestriction(*this, e_ordering));
1283  }
1284  return L2E_lex.Ptr();
1285  }
1286  // e_ordering == ElementDofOrdering::NATIVE
1287  if (L2E_nat.Ptr() == NULL)
1288  {
1289  L2E_nat.Reset(new ElementRestriction(*this, e_ordering));
1290  }
1291  return L2E_nat.Ptr();
1292 }
1293 
1295  ElementDofOrdering e_ordering, FaceType type, L2FaceValues mul) const
1296 {
1297  const bool is_dg_space = IsDGSpace();
1298  const L2FaceValues m = (is_dg_space && mul==L2FaceValues::DoubleValued) ?
1300  key_face key = std::make_tuple(is_dg_space, e_ordering, type, m);
1301  auto itr = L2F.find(key);
1302  if (itr != L2F.end())
1303  {
1304  return itr->second;
1305  }
1306  else
1307  {
1308  FaceRestriction *res;
1309  if (is_dg_space)
1310  {
1311  if (Conforming())
1312  {
1313  res = new L2FaceRestriction(*this, e_ordering, type, m);
1314  }
1315  else
1316  {
1317  res = new NCL2FaceRestriction(*this, e_ordering, type, m);
1318  }
1319  }
1320  else
1321  {
1322  res = new H1FaceRestriction(*this, e_ordering, type);
1323  }
1324  L2F[key] = res;
1325  return res;
1326  }
1327 }
1328 
1330  const IntegrationRule &ir) const
1331 {
1332  for (int i = 0; i < E2Q_array.Size(); i++)
1333  {
1334  const QuadratureInterpolator *qi = E2Q_array[i];
1335  if (qi->IntRule == &ir) { return qi; }
1336  }
1337 
1338  QuadratureInterpolator *qi = new QuadratureInterpolator(*this, ir);
1339  E2Q_array.Append(qi);
1340  return qi;
1341 }
1342 
1344  const QuadratureSpace &qs) const
1345 {
1346  for (int i = 0; i < E2Q_array.Size(); i++)
1347  {
1348  const QuadratureInterpolator *qi = E2Q_array[i];
1349  if (qi->qspace == &qs) { return qi; }
1350  }
1351 
1352  QuadratureInterpolator *qi = new QuadratureInterpolator(*this, qs);
1353  E2Q_array.Append(qi);
1354  return qi;
1355 }
1356 
1359  const IntegrationRule &ir, FaceType type) const
1360 {
1361  if (type==FaceType::Interior)
1362  {
1363  for (int i = 0; i < E2IFQ_array.Size(); i++)
1364  {
1365  const FaceQuadratureInterpolator *qi = E2IFQ_array[i];
1366  if (qi->IntRule == &ir) { return qi; }
1367  }
1368 
1370  type);
1371  E2IFQ_array.Append(qi);
1372  return qi;
1373  }
1374  else //Boundary
1375  {
1376  for (int i = 0; i < E2BFQ_array.Size(); i++)
1377  {
1378  const FaceQuadratureInterpolator *qi = E2BFQ_array[i];
1379  if (qi->IntRule == &ir) { return qi; }
1380  }
1381 
1383  type);
1384  E2BFQ_array.Append(qi);
1385  return qi;
1386  }
1387 }
1388 
1390  const int coarse_ndofs, const Table &coarse_elem_dof,
1391  const Table *coarse_elem_fos, const DenseTensor localP[]) const
1392 {
1393  /// TODO: Implement DofTransformation support
1394 
1395  MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
1396 
1397  Array<int> dofs, coarse_dofs, coarse_vdofs;
1398  Vector row;
1399 
1400  Mesh::GeometryList elem_geoms(*mesh);
1401 
1402  SparseMatrix *P;
1403  if (elem_geoms.Size() == 1)
1404  {
1405  const int coarse_ldof = localP[elem_geoms[0]].SizeJ();
1406  P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim, coarse_ldof);
1407  }
1408  else
1409  {
1410  P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim);
1411  }
1412 
1413  Array<int> mark(P->Height());
1414  mark = 0;
1415 
1417 
1418  for (int k = 0; k < mesh->GetNE(); k++)
1419  {
1420  const Embedding &emb = rtrans.embeddings[k];
1421  const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
1422  const DenseMatrix &lP = localP[geom](emb.matrix);
1423  const int fine_ldof = localP[geom].SizeI();
1424 
1425  elem_dof->GetRow(k, dofs);
1426  coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
1427 
1428  for (int vd = 0; vd < vdim; vd++)
1429  {
1430  coarse_dofs.Copy(coarse_vdofs);
1431  DofsToVDofs(vd, coarse_vdofs, coarse_ndofs);
1432 
1433  for (int i = 0; i < fine_ldof; i++)
1434  {
1435  int r = DofToVDof(dofs[i], vd);
1436  int m = (r >= 0) ? r : (-1 - r);
1437 
1438  if (!mark[m])
1439  {
1440  lP.GetRow(i, row);
1441  P->SetRow(r, coarse_vdofs, row);
1442  mark[m] = 1;
1443  }
1444  }
1445  }
1446  }
1447 
1448  MFEM_ASSERT(mark.Sum() == P->Height(), "Not all rows of P set.");
1449  if (elem_geoms.Size() != 1) { P->Finalize(); }
1450  return P;
1451 }
1452 
1454  Geometry::Type geom, DenseTensor &localP) const
1455 {
1456  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
1457 
1459  const DenseTensor &pmats = rtrans.point_matrices[geom];
1460 
1461  int nmat = pmats.SizeK();
1462  int ldof = fe->GetDof();
1463 
1465  isotr.SetIdentityTransformation(geom);
1466 
1467  // calculate local interpolation matrices for all refinement types
1468  localP.SetSize(ldof, ldof, nmat);
1469  for (int i = 0; i < nmat; i++)
1470  {
1471  isotr.SetPointMat(pmats(i));
1472  fe->GetLocalInterpolation(isotr, localP(i));
1473  }
1474 }
1475 
1477  const Table* old_elem_dof,
1478  const Table* old_elem_fos)
1479 {
1480  MFEM_VERIFY(GetNE() >= old_elem_dof->Size(),
1481  "Previous mesh is not coarser.");
1482 
1483  Mesh::GeometryList elem_geoms(*mesh);
1484 
1486  for (int i = 0; i < elem_geoms.Size(); i++)
1487  {
1488  GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1489  }
1490 
1491  return RefinementMatrix_main(old_ndofs, *old_elem_dof, old_elem_fos,
1492  localP);
1493 }
1494 
1496 (const FiniteElementSpace* fespace, Table* old_elem_dof, Table* old_elem_fos,
1497  int old_ndofs)
1498  : fespace(fespace)
1499  , old_elem_dof(old_elem_dof)
1500  , old_elem_fos(old_elem_fos)
1501 {
1502  MFEM_VERIFY(fespace->GetNE() >= old_elem_dof->Size(),
1503  "Previous mesh is not coarser.");
1504 
1505  width = old_ndofs * fespace->GetVDim();
1506  height = fespace->GetVSize();
1507 
1508  Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1509 
1510  for (int i = 0; i < elem_geoms.Size(); i++)
1511  {
1512  fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1513  }
1514 
1515  ConstructDoFTrans();
1516 }
1517 
1519  const FiniteElementSpace *fespace, const FiniteElementSpace *coarse_fes)
1520  : Operator(fespace->GetVSize(), coarse_fes->GetVSize()),
1521  fespace(fespace), old_elem_dof(NULL), old_elem_fos(NULL)
1522 {
1523  Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1524 
1525  for (int i = 0; i < elem_geoms.Size(); i++)
1526  {
1527  fespace->GetLocalRefinementMatrices(*coarse_fes, elem_geoms[i],
1528  localP[elem_geoms[i]]);
1529  }
1530 
1531  // Make a copy of the coarse elem_dof Table.
1532  old_elem_dof = new Table(coarse_fes->GetElementToDofTable());
1533 
1534  // Make a copy of the coarse elem_fos Table if it exists.
1535  if (coarse_fes->GetElementToFaceOrientationTable())
1536  {
1537  old_elem_fos = new Table(*coarse_fes->GetElementToFaceOrientationTable());
1538  }
1539 
1540  ConstructDoFTrans();
1541 }
1542 
1544 {
1545  delete old_elem_dof;
1546  delete old_elem_fos;
1547 }
1548 
1549 void FiniteElementSpace::RefinementOperator
1550 ::ConstructDoFTrans()
1551 {
1552  old_DoFTrans.SetSize(Geometry::NUM_GEOMETRIES);
1553  for (int i=0; i<old_DoFTrans.Size(); i++)
1554  {
1555  old_DoFTrans[i] = NULL;
1556  }
1557 
1558  const FiniteElementCollection *fec_ref = fespace->FEColl();
1559  if (dynamic_cast<const ND_FECollection*>(fec_ref))
1560  {
1561  const FiniteElement * nd_tri =
1563  if (nd_tri)
1564  {
1565  old_DoFTrans[Geometry::TRIANGLE] =
1566  new ND_TriDofTransformation(nd_tri->GetOrder());
1567  }
1568 
1569  const FiniteElement * nd_tet =
1571  if (nd_tet)
1572  {
1573  old_DoFTrans[Geometry::TETRAHEDRON] =
1574  new ND_TetDofTransformation(nd_tet->GetOrder());
1575  }
1576 
1577  const FiniteElement * nd_pri =
1579  if (nd_pri)
1580  {
1581  old_DoFTrans[Geometry::PRISM] =
1582  new ND_WedgeDofTransformation(nd_pri->GetOrder());
1583  }
1584  }
1585 }
1586 
1588 ::Mult(const Vector &x, Vector &y) const
1589 {
1590  Mesh* mesh_ref = fespace->GetMesh();
1591  const CoarseFineTransformations &trans_ref =
1592  mesh_ref->GetRefinementTransforms();
1593 
1594  Array<int> dofs, vdofs, old_dofs, old_vdofs, old_Fo;
1595 
1596  int rvdim = fespace->GetVDim();
1597  int old_ndofs = width / rvdim;
1598 
1599  Vector subY, subX;
1600 
1601  for (int k = 0; k < mesh_ref->GetNE(); k++)
1602  {
1603  const Embedding &emb = trans_ref.embeddings[k];
1604  const Geometry::Type geom = mesh_ref->GetElementBaseGeometry(k);
1605  const DenseMatrix &lP = localP[geom](emb.matrix);
1606 
1607  subY.SetSize(lP.Height());
1608 
1609  DofTransformation *doftrans = fespace->GetElementDofs(k, dofs);
1610  old_elem_dof->GetRow(emb.parent, old_dofs);
1611 
1612  if (!doftrans)
1613  {
1614  for (int vd = 0; vd < rvdim; vd++)
1615  {
1616  dofs.Copy(vdofs);
1617  fespace->DofsToVDofs(vd, vdofs);
1618  old_dofs.Copy(old_vdofs);
1619  fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
1620  x.GetSubVector(old_vdofs, subX);
1621  lP.Mult(subX, subY);
1622  y.SetSubVector(vdofs, subY);
1623  }
1624  }
1625  else
1626  {
1627  old_elem_fos->GetRow(emb.parent, old_Fo);
1628  old_DoFTrans[geom]->SetFaceOrientations(old_Fo);
1629 
1630  DofTransformation *new_doftrans = NULL;
1631  VDofTransformation *vdoftrans =
1632  dynamic_cast<VDofTransformation*>(doftrans);
1633  if (vdoftrans)
1634  {
1635  new_doftrans = doftrans;
1636  doftrans = vdoftrans->GetDofTransformation();
1637  }
1638 
1639  for (int vd = 0; vd < rvdim; vd++)
1640  {
1641  dofs.Copy(vdofs);
1642  fespace->DofsToVDofs(vd, vdofs);
1643  old_dofs.Copy(old_vdofs);
1644  fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
1645  x.GetSubVector(old_vdofs, subX);
1646  old_DoFTrans[geom]->InvTransformPrimal(subX);
1647  lP.Mult(subX, subY);
1648  doftrans->TransformPrimal(subY);
1649  y.SetSubVector(vdofs, subY);
1650  }
1651 
1652  if (vdoftrans)
1653  {
1654  doftrans = new_doftrans;
1655  }
1656  }
1657  }
1658 }
1659 
1661 ::MultTranspose(const Vector &x, Vector &y) const
1662 {
1663  y = 0.0;
1664 
1665  Mesh* mesh_ref = fespace->GetMesh();
1666  const CoarseFineTransformations &trans_ref =
1667  mesh_ref->GetRefinementTransforms();
1668 
1669  Array<char> processed(fespace->GetVSize());
1670  processed = 0;
1671 
1672  Array<int> f_dofs, c_dofs, f_vdofs, c_vdofs, old_Fo;
1673 
1674  int rvdim = fespace->GetVDim();
1675  int old_ndofs = width / rvdim;
1676 
1677  Vector subY, subX, subYt;
1678 
1679  for (int k = 0; k < mesh_ref->GetNE(); k++)
1680  {
1681  const Embedding &emb = trans_ref.embeddings[k];
1682  const Geometry::Type geom = mesh_ref->GetElementBaseGeometry(k);
1683  const DenseMatrix &lP = localP[geom](emb.matrix);
1684 
1685  DofTransformation * doftrans = fespace->GetElementDofs(k, f_dofs);
1686  old_elem_dof->GetRow(emb.parent, c_dofs);
1687 
1688  if (!doftrans)
1689  {
1690  subY.SetSize(lP.Width());
1691 
1692  for (int vd = 0; vd < rvdim; vd++)
1693  {
1694  f_dofs.Copy(f_vdofs);
1695  fespace->DofsToVDofs(vd, f_vdofs);
1696  c_dofs.Copy(c_vdofs);
1697  fespace->DofsToVDofs(vd, c_vdofs, old_ndofs);
1698 
1699  x.GetSubVector(f_vdofs, subX);
1700 
1701  for (int p = 0; p < f_dofs.Size(); ++p)
1702  {
1703  if (processed[DecodeDof(f_dofs[p])])
1704  {
1705  subX[p] = 0.0;
1706  }
1707  }
1708 
1709  lP.MultTranspose(subX, subY);
1710  y.AddElementVector(c_vdofs, subY);
1711  }
1712  }
1713  else
1714  {
1715  subYt.SetSize(lP.Width());
1716 
1717  old_elem_fos->GetRow(emb.parent, old_Fo);
1718  old_DoFTrans[geom]->SetFaceOrientations(old_Fo);
1719 
1720  DofTransformation *new_doftrans = NULL;
1721  VDofTransformation *vdoftrans =
1722  dynamic_cast<VDofTransformation*>(doftrans);
1723  if (vdoftrans)
1724  {
1725  new_doftrans = doftrans;
1726  doftrans = vdoftrans->GetDofTransformation();
1727  }
1728 
1729  for (int vd = 0; vd < rvdim; vd++)
1730  {
1731  f_dofs.Copy(f_vdofs);
1732  fespace->DofsToVDofs(vd, f_vdofs);
1733  c_dofs.Copy(c_vdofs);
1734  fespace->DofsToVDofs(vd, c_vdofs, old_ndofs);
1735 
1736  x.GetSubVector(f_vdofs, subX);
1737  doftrans->InvTransformDual(subX);
1738  for (int p = 0; p < f_dofs.Size(); ++p)
1739  {
1740  if (processed[DecodeDof(f_dofs[p])])
1741  {
1742  subX[p] = 0.0;
1743  }
1744  }
1745 
1746  lP.MultTranspose(subX, subYt);
1747  old_DoFTrans[geom]->TransformDual(subYt);
1748  y.AddElementVector(c_vdofs, subYt);
1749  }
1750 
1751  if (vdoftrans)
1752  {
1753  doftrans = new_doftrans;
1754  }
1755  }
1756 
1757  for (int p = 0; p < f_dofs.Size(); ++p)
1758  {
1759  processed[DecodeDof(f_dofs[p])] = 1;
1760  }
1761  }
1762 }
1763 
1764 namespace internal
1765 {
1766 
1767 // Used in GetCoarseToFineMap() below.
1768 struct RefType
1769 {
1770  Geometry::Type geom;
1771  int num_children;
1772  const Pair<int,int> *children;
1773 
1774  RefType(Geometry::Type g, int n, const Pair<int,int> *c)
1775  : geom(g), num_children(n), children(c) { }
1776 
1777  bool operator<(const RefType &other) const
1778  {
1779  if (geom < other.geom) { return true; }
1780  if (geom > other.geom) { return false; }
1781  if (num_children < other.num_children) { return true; }
1782  if (num_children > other.num_children) { return false; }
1783  for (int i = 0; i < num_children; i++)
1784  {
1785  if (children[i].one < other.children[i].one) { return true; }
1786  if (children[i].one > other.children[i].one) { return false; }
1787  }
1788  return false; // everything is equal
1789  }
1790 };
1791 
1792 void GetCoarseToFineMap(const CoarseFineTransformations &cft,
1793  const mfem::Mesh &fine_mesh,
1794  Table &coarse_to_fine,
1795  Array<int> &coarse_to_ref_type,
1796  Table &ref_type_to_matrix,
1797  Array<Geometry::Type> &ref_type_to_geom)
1798 {
1799  const int fine_ne = cft.embeddings.Size();
1800  int coarse_ne = -1;
1801  for (int i = 0; i < fine_ne; i++)
1802  {
1803  coarse_ne = std::max(coarse_ne, cft.embeddings[i].parent);
1804  }
1805  coarse_ne++;
1806 
1807  coarse_to_ref_type.SetSize(coarse_ne);
1808  coarse_to_fine.SetDims(coarse_ne, fine_ne);
1809 
1810  Array<int> cf_i(coarse_to_fine.GetI(), coarse_ne+1);
1811  Array<Pair<int,int> > cf_j(fine_ne);
1812  cf_i = 0;
1813  for (int i = 0; i < fine_ne; i++)
1814  {
1815  cf_i[cft.embeddings[i].parent+1]++;
1816  }
1817  cf_i.PartialSum();
1818  MFEM_ASSERT(cf_i.Last() == cf_j.Size(), "internal error");
1819  for (int i = 0; i < fine_ne; i++)
1820  {
1821  const Embedding &e = cft.embeddings[i];
1822  cf_j[cf_i[e.parent]].one = e.matrix; // used as sort key below
1823  cf_j[cf_i[e.parent]].two = i;
1824  cf_i[e.parent]++;
1825  }
1826  std::copy_backward(cf_i.begin(), cf_i.end()-1, cf_i.end());
1827  cf_i[0] = 0;
1828  for (int i = 0; i < coarse_ne; i++)
1829  {
1830  std::sort(&cf_j[cf_i[i]], cf_j.GetData() + cf_i[i+1]);
1831  }
1832  for (int i = 0; i < fine_ne; i++)
1833  {
1834  coarse_to_fine.GetJ()[i] = cf_j[i].two;
1835  }
1836 
1837  using std::map;
1838  using std::pair;
1839 
1840  map<RefType,int> ref_type_map;
1841  for (int i = 0; i < coarse_ne; i++)
1842  {
1843  const int num_children = cf_i[i+1]-cf_i[i];
1844  MFEM_ASSERT(num_children > 0, "");
1845  const int fine_el = cf_j[cf_i[i]].two;
1846  // Assuming the coarse and the fine elements have the same geometry:
1847  const Geometry::Type geom = fine_mesh.GetElementBaseGeometry(fine_el);
1848  const RefType ref_type(geom, num_children, &cf_j[cf_i[i]]);
1849  pair<map<RefType,int>::iterator,bool> res =
1850  ref_type_map.insert(
1851  pair<const RefType,int>(ref_type, (int)ref_type_map.size()));
1852  coarse_to_ref_type[i] = res.first->second;
1853  }
1854 
1855  ref_type_to_matrix.MakeI((int)ref_type_map.size());
1856  ref_type_to_geom.SetSize((int)ref_type_map.size());
1857  for (map<RefType,int>::iterator it = ref_type_map.begin();
1858  it != ref_type_map.end(); ++it)
1859  {
1860  ref_type_to_matrix.AddColumnsInRow(it->second, it->first.num_children);
1861  ref_type_to_geom[it->second] = it->first.geom;
1862  }
1863 
1864  ref_type_to_matrix.MakeJ();
1865  for (map<RefType,int>::iterator it = ref_type_map.begin();
1866  it != ref_type_map.end(); ++it)
1867  {
1868  const RefType &rt = it->first;
1869  for (int j = 0; j < rt.num_children; j++)
1870  {
1871  ref_type_to_matrix.AddConnection(it->second, rt.children[j].one);
1872  }
1873  }
1874  ref_type_to_matrix.ShiftUpI();
1875 }
1876 
1877 } // namespace internal
1878 
1879 
1880 /// TODO: Implement DofTransformation support
1882  const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes,
1883  BilinearFormIntegrator *mass_integ)
1884  : Operator(c_fes->GetVSize(), f_fes->GetVSize()),
1885  fine_fes(f_fes)
1886 {
1887  MFEM_VERIFY(c_fes->GetOrdering() == f_fes->GetOrdering() &&
1888  c_fes->GetVDim() == f_fes->GetVDim(),
1889  "incompatible coarse and fine FE spaces");
1890 
1892  Mesh *f_mesh = f_fes->GetMesh();
1893  const CoarseFineTransformations &rtrans = f_mesh->GetRefinementTransforms();
1894 
1895  Mesh::GeometryList elem_geoms(*f_mesh);
1897  for (int gi = 0; gi < elem_geoms.Size(); gi++)
1898  {
1899  const Geometry::Type geom = elem_geoms[gi];
1900  DenseTensor &lP = localP[geom], &lM = localM[geom];
1901  const FiniteElement *fine_fe =
1902  f_fes->fec->FiniteElementForGeometry(geom);
1903  const FiniteElement *coarse_fe =
1904  c_fes->fec->FiniteElementForGeometry(geom);
1905  const DenseTensor &pmats = rtrans.point_matrices[geom];
1906 
1907  lP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), pmats.SizeK());
1908  lM.SetSize(fine_fe->GetDof(), fine_fe->GetDof(), pmats.SizeK());
1909  emb_tr.SetIdentityTransformation(geom);
1910  for (int i = 0; i < pmats.SizeK(); i++)
1911  {
1912  emb_tr.SetPointMat(pmats(i));
1913  // Get the local interpolation matrix for this refinement type
1914  fine_fe->GetTransferMatrix(*coarse_fe, emb_tr, lP(i));
1915  // Get the local mass matrix for this refinement type
1916  mass_integ->AssembleElementMatrix(*fine_fe, emb_tr, lM(i));
1917  }
1918  }
1919 
1920  Table ref_type_to_matrix;
1921  internal::GetCoarseToFineMap(rtrans, *f_mesh, coarse_to_fine,
1922  coarse_to_ref_type, ref_type_to_matrix,
1923  ref_type_to_geom);
1924  MFEM_ASSERT(coarse_to_fine.Size() == c_fes->GetNE(), "");
1925 
1926  const int total_ref_types = ref_type_to_geom.Size();
1927  int num_ref_types[Geometry::NumGeom], num_fine_elems[Geometry::NumGeom];
1928  Array<int> ref_type_to_coarse_elem_offset(total_ref_types);
1929  ref_type_to_fine_elem_offset.SetSize(total_ref_types);
1930  std::fill(num_ref_types, num_ref_types+Geometry::NumGeom, 0);
1931  std::fill(num_fine_elems, num_fine_elems+Geometry::NumGeom, 0);
1932  for (int i = 0; i < total_ref_types; i++)
1933  {
1934  Geometry::Type g = ref_type_to_geom[i];
1935  ref_type_to_coarse_elem_offset[i] = num_ref_types[g];
1936  ref_type_to_fine_elem_offset[i] = num_fine_elems[g];
1937  num_ref_types[g]++;
1938  num_fine_elems[g] += ref_type_to_matrix.RowSize(i);
1939  }
1940  DenseTensor localPtMP[Geometry::NumGeom];
1941  for (int g = 0; g < Geometry::NumGeom; g++)
1942  {
1943  if (num_ref_types[g] == 0) { continue; }
1944  const int fine_dofs = localP[g].SizeI();
1945  const int coarse_dofs = localP[g].SizeJ();
1946  localPtMP[g].SetSize(coarse_dofs, coarse_dofs, num_ref_types[g]);
1947  localR[g].SetSize(coarse_dofs, fine_dofs, num_fine_elems[g]);
1948  }
1949  for (int i = 0; i < total_ref_types; i++)
1950  {
1951  Geometry::Type g = ref_type_to_geom[i];
1952  DenseMatrix &lPtMP = localPtMP[g](ref_type_to_coarse_elem_offset[i]);
1953  int lR_offset = ref_type_to_fine_elem_offset[i]; // offset in localR[g]
1954  const int *mi = ref_type_to_matrix.GetRow(i);
1955  const int nm = ref_type_to_matrix.RowSize(i);
1956  lPtMP = 0.0;
1957  for (int s = 0; s < nm; s++)
1958  {
1959  DenseMatrix &lP = localP[g](mi[s]);
1960  DenseMatrix &lM = localM[g](mi[s]);
1961  DenseMatrix &lR = localR[g](lR_offset+s);
1962  MultAtB(lP, lM, lR); // lR = lP^T lM
1963  AddMult(lR, lP, lPtMP); // lPtMP += lP^T lM lP
1964  }
1965  DenseMatrixInverse lPtMP_inv(lPtMP);
1966  for (int s = 0; s < nm; s++)
1967  {
1968  DenseMatrix &lR = localR[g](lR_offset+s);
1969  lPtMP_inv.Mult(lR); // lR <- (P^T M P)^{-1} P^T M
1970  }
1971  }
1972 
1973  // Make a copy of the coarse element-to-dof Table.
1974  coarse_elem_dof = new Table(c_fes->GetElementToDofTable());
1975 }
1976 
1978 {
1979  delete coarse_elem_dof;
1980 }
1981 
1983 ::Mult(const Vector &x, Vector &y) const
1984 {
1985  Array<int> c_vdofs, f_vdofs;
1986  Vector loc_x, loc_y;
1987  DenseMatrix loc_x_mat, loc_y_mat;
1988  const int fine_vdim = fine_fes->GetVDim();
1989  const int coarse_ndofs = height/fine_vdim;
1990  for (int coarse_el = 0; coarse_el < coarse_to_fine.Size(); coarse_el++)
1991  {
1992  coarse_elem_dof->GetRow(coarse_el, c_vdofs);
1993  fine_fes->DofsToVDofs(c_vdofs, coarse_ndofs);
1994  loc_y.SetSize(c_vdofs.Size());
1995  loc_y = 0.0;
1996  loc_y_mat.UseExternalData(loc_y.GetData(), c_vdofs.Size()/fine_vdim,
1997  fine_vdim);
1998  const int ref_type = coarse_to_ref_type[coarse_el];
1999  const Geometry::Type geom = ref_type_to_geom[ref_type];
2000  const int *fine_elems = coarse_to_fine.GetRow(coarse_el);
2001  const int num_fine_elems = coarse_to_fine.RowSize(coarse_el);
2002  const int lR_offset = ref_type_to_fine_elem_offset[ref_type];
2003  for (int s = 0; s < num_fine_elems; s++)
2004  {
2005  const DenseMatrix &lR = localR[geom](lR_offset+s);
2006  fine_fes->GetElementVDofs(fine_elems[s], f_vdofs);
2007  x.GetSubVector(f_vdofs, loc_x);
2008  loc_x_mat.UseExternalData(loc_x.GetData(), f_vdofs.Size()/fine_vdim,
2009  fine_vdim);
2010  AddMult(lR, loc_x_mat, loc_y_mat);
2011  }
2012  y.SetSubVector(c_vdofs, loc_y);
2013  }
2014 }
2015 
2017  DenseTensor &localR) const
2018 {
2019  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
2020 
2021  const CoarseFineTransformations &dtrans =
2023  const DenseTensor &pmats = dtrans.point_matrices[geom];
2024 
2025  const int nmat = pmats.SizeK();
2026  const int ldof = fe->GetDof();
2027 
2029  isotr.SetIdentityTransformation(geom);
2030 
2031  // calculate local restriction matrices for all refinement types
2032  localR.SetSize(ldof, ldof, nmat);
2033  for (int i = 0; i < nmat; i++)
2034  {
2035  isotr.SetPointMat(pmats(i));
2036  fe->GetLocalRestriction(isotr, localR(i));
2037  }
2038 }
2039 
2041  const Table* old_elem_dof,
2042  const Table* old_elem_fos)
2043 {
2044  /// TODO: Implement DofTransformation support
2045 
2046  MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
2047  MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
2048  MFEM_VERIFY(ndofs <= old_ndofs, "Previous space is not finer.");
2049 
2050  Array<int> dofs, old_dofs, old_vdofs;
2051  Vector row;
2052 
2053  Mesh::GeometryList elem_geoms(*mesh);
2054 
2056  for (int i = 0; i < elem_geoms.Size(); i++)
2057  {
2058  GetLocalDerefinementMatrices(elem_geoms[i], localR[elem_geoms[i]]);
2059  }
2060 
2061  SparseMatrix *R = (elem_geoms.Size() != 1)
2062  ? new SparseMatrix(ndofs*vdim, old_ndofs*vdim) // variable row size
2063  : new SparseMatrix(ndofs*vdim, old_ndofs*vdim,
2064  localR[elem_geoms[0]].SizeI());
2065 
2066  Array<int> mark(R->Height());
2067  mark = 0;
2068 
2069  const CoarseFineTransformations &dtrans =
2071 
2072  MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
2073 
2074  int num_marked = 0;
2075  for (int k = 0; k < dtrans.embeddings.Size(); k++)
2076  {
2077  const Embedding &emb = dtrans.embeddings[k];
2079  DenseMatrix &lR = localR[geom](emb.matrix);
2080 
2081  elem_dof->GetRow(emb.parent, dofs);
2082  old_elem_dof->GetRow(k, old_dofs);
2083 
2084  for (int vd = 0; vd < vdim; vd++)
2085  {
2086  old_dofs.Copy(old_vdofs);
2087  DofsToVDofs(vd, old_vdofs, old_ndofs);
2088 
2089  for (int i = 0; i < lR.Height(); i++)
2090  {
2091  if (!std::isfinite(lR(i, 0))) { continue; }
2092 
2093  int r = DofToVDof(dofs[i], vd);
2094  int m = (r >= 0) ? r : (-1 - r);
2095 
2096  if (!mark[m])
2097  {
2098  lR.GetRow(i, row);
2099  R->SetRow(r, old_vdofs, row);
2100  mark[m] = 1;
2101  num_marked++;
2102  }
2103  }
2104  }
2105  }
2106 
2107  MFEM_VERIFY(num_marked == R->Height(),
2108  "internal error: not all rows of R were set.");
2109 
2110  R->Finalize(); // no-op if fixed width
2111  return R;
2112 }
2113 
2115  const FiniteElementSpace &coarse_fes, Geometry::Type geom,
2116  DenseTensor &localP) const
2117 {
2118  // Assumptions: see the declaration of the method.
2119 
2120  const FiniteElement *fine_fe = fec->FiniteElementForGeometry(geom);
2121  const FiniteElement *coarse_fe =
2122  coarse_fes.fec->FiniteElementForGeometry(geom);
2123 
2125  const DenseTensor &pmats = rtrans.point_matrices[geom];
2126 
2127  int nmat = pmats.SizeK();
2128 
2130  isotr.SetIdentityTransformation(geom);
2131 
2132  // Calculate the local interpolation matrices for all refinement types
2133  localP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), nmat);
2134  for (int i = 0; i < nmat; i++)
2135  {
2136  isotr.SetPointMat(pmats(i));
2137  fine_fe->GetTransferMatrix(*coarse_fe, isotr, localP(i));
2138  }
2139 }
2140 
2142  const FiniteElementCollection *fec_,
2143  int vdim_, int ordering_)
2144 {
2145  mesh = mesh_;
2146  fec = fec_;
2147  vdim = vdim_;
2148  ordering = (Ordering::Type) ordering_;
2149 
2150  elem_dof = NULL;
2151  elem_fos = NULL;
2152  face_dof = NULL;
2153 
2154  sequence = 0;
2155  orders_changed = false;
2156  relaxed_hp = false;
2157 
2159 
2160  const NURBSFECollection *nurbs_fec =
2161  dynamic_cast<const NURBSFECollection *>(fec_);
2162  if (nurbs_fec)
2163  {
2164  MFEM_VERIFY(mesh_->NURBSext, "NURBS FE space requires a NURBS mesh.");
2165 
2166  if (NURBSext_ == NULL)
2167  {
2168  NURBSext = mesh_->NURBSext;
2169  own_ext = 0;
2170  }
2171  else
2172  {
2173  NURBSext = NURBSext_;
2174  own_ext = 1;
2175  }
2176  UpdateNURBS();
2177  cP = cR = cR_hp = NULL;
2178  cP_is_set = false;
2179 
2181  }
2182  else
2183  {
2184  NURBSext = NULL;
2185  own_ext = 0;
2186  Construct();
2187  }
2188 
2190 }
2191 
2193 {
2194  DestroyDoFTrans();
2195 
2198  for (int i=0; i<DoFTrans.Size(); i++)
2199  {
2200  DoFTrans[i] = NULL;
2201  }
2202  if (mesh->Dimension() < 3) { return; }
2203  if (dynamic_cast<const ND_FECollection*>(fec))
2204  {
2205  const FiniteElement * nd_tri =
2207  if (nd_tri)
2208  {
2210  new ND_TriDofTransformation(nd_tri->GetOrder());
2211  }
2212 
2213  const FiniteElement * nd_tet =
2215  if (nd_tet)
2216  {
2218  new ND_TetDofTransformation(nd_tet->GetOrder());
2219  }
2220 
2221  const FiniteElement * nd_pri =
2223  if (nd_pri)
2224  {
2226  new ND_WedgeDofTransformation(nd_pri->GetOrder());
2227  }
2228  }
2229 }
2230 
2232 {
2233  if (NURBSext && !own_ext)
2234  {
2235  mfem_error("FiniteElementSpace::StealNURBSext");
2236  }
2237  own_ext = 0;
2238 
2239  return NURBSext;
2240 }
2241 
2243 {
2244  MFEM_VERIFY(NURBSext, "NURBSExt not defined.");
2245 
2246  nvdofs = 0;
2247  nedofs = 0;
2248  nfdofs = 0;
2249  nbdofs = 0;
2250  bdofs = NULL;
2251 
2252  delete face_dof;
2253  face_dof = NULL;
2255 
2256  dynamic_cast<const NURBSFECollection *>(fec)->Reset();
2257 
2258  ndofs = NURBSext->GetNDof();
2261 
2263  sequence++;
2264 }
2265 
2267 {
2268  if (face_dof) { return; }
2269 
2270  const int dim = mesh->Dimension();
2271 
2272  // Find bdr to face mapping
2274  face_to_be = -1;
2275  for (int b = 0; b < GetNBE(); b++)
2276  {
2277  int f = mesh->GetBdrElementEdgeIndex(b);
2278  face_to_be[f] = b;
2279  }
2280 
2281  // Loop over faces in correct order, to prevent a sort
2282  // Sort will destroy orientation info in ordering of dofs
2283  Array<Connection> face_dof_list;
2284  Array<int> row;
2285  for (int f = 0; f < GetNF(); f++)
2286  {
2287  int b = face_to_be[f];
2288  if (b == -1) { continue; }
2289  // FIXME: this assumes that the boundary element and the face element have
2290  // the same orientation.
2291  if (dim > 1)
2292  {
2293  const Element *fe = mesh->GetFace(f);
2294  const Element *be = mesh->GetBdrElement(b);
2295  const int nv = be->GetNVertices();
2296  const int *fv = fe->GetVertices();
2297  const int *bv = be->GetVertices();
2298  for (int i = 0; i < nv; i++)
2299  {
2300  MFEM_VERIFY(fv[i] == bv[i],
2301  "non-matching face and boundary elements detected!");
2302  }
2303  }
2304  GetBdrElementDofs(b, row);
2305  Connection conn(f,0);
2306  for (int i = 0; i < row.Size(); i++)
2307  {
2308  conn.to = row[i];
2309  face_dof_list.Append(conn);
2310  }
2311  }
2312  face_dof = new Table(GetNF(), face_dof_list);
2313 }
2314 
2316 {
2317  // This method should be used only for non-NURBS spaces.
2318  MFEM_VERIFY(!NURBSext, "internal error");
2319 
2320  // Variable order space needs a nontrivial P matrix + also ghost elements
2321  // in parallel, we thus require the mesh to be NC.
2322  MFEM_VERIFY(!IsVariableOrder() || Nonconforming(),
2323  "Variable order space requires a nonconforming mesh.");
2324 
2325  elem_dof = NULL;
2326  elem_fos = NULL;
2327  bdr_elem_dof = NULL;
2328  bdr_elem_fos = NULL;
2329  face_dof = NULL;
2330 
2331  ndofs = 0;
2332  nvdofs = nedofs = nfdofs = nbdofs = 0;
2333  bdofs = NULL;
2334 
2335  cP = NULL;
2336  cR = NULL;
2337  cR_hp = NULL;
2338  cP_is_set = false;
2339  // 'Th' is initialized/destroyed before this method is called.
2340 
2341  int dim = mesh->Dimension();
2342  int order = fec->GetOrder();
2343 
2344  MFEM_VERIFY((mesh->GetNumGeometries(dim) > 0) || (mesh->GetNE() == 0),
2345  "Mesh was not correctly finalized.");
2346 
2347  bool mixed_elements = (mesh->GetNumGeometries(dim) > 1);
2348  bool mixed_faces = (dim > 2 && mesh->GetNumGeometries(2) > 1);
2349 
2350  Array<VarOrderBits> edge_orders, face_orders;
2351  if (IsVariableOrder())
2352  {
2353  // for variable order spaces, calculate orders of edges and faces
2354  CalcEdgeFaceVarOrders(edge_orders, face_orders);
2355  }
2356  else if (mixed_faces)
2357  {
2358  // for mixed faces we also create the var_face_dofs table, see below
2359  face_orders.SetSize(mesh->GetNFaces());
2360  face_orders = (VarOrderBits(1) << order);
2361  }
2362 
2363  // assign vertex DOFs
2364  if (mesh->GetNV())
2365  {
2366  nvdofs = mesh->GetNV() * fec->GetNumDof(Geometry::POINT, order);
2367  }
2368 
2369  // assign edge DOFs
2370  if (mesh->GetNEdges())
2371  {
2372  if (IsVariableOrder())
2373  {
2374  nedofs = MakeDofTable(1, edge_orders, var_edge_dofs, &var_edge_orders);
2375  }
2376  else
2377  {
2378  // the simple case: all edges are of the same order
2380  }
2381  }
2382 
2383  // assign face DOFs
2384  if (mesh->GetNFaces())
2385  {
2386  if (IsVariableOrder() || mixed_faces)
2387  {
2388  // NOTE: for simplicity, we also use Table var_face_dofs for mixed faces
2389  nfdofs = MakeDofTable(2, face_orders, var_face_dofs,
2390  IsVariableOrder() ? &var_face_orders : NULL);
2391  uni_fdof = -1;
2392  }
2393  else
2394  {
2395  // the simple case: all faces are of the same geometry and order
2396  uni_fdof = fec->GetNumDof(mesh->GetFaceGeometry(0), order);
2397  nfdofs = mesh->GetNFaces() * uni_fdof;
2398  }
2399  }
2400 
2401  // assign internal ("bubble") DOFs
2402  if (mesh->GetNE() && dim > 0)
2403  {
2404  if (IsVariableOrder() || mixed_elements)
2405  {
2406  bdofs = new int[mesh->GetNE()+1];
2407  bdofs[0] = 0;
2408  for (int i = 0; i < mesh->GetNE(); i++)
2409  {
2410  int p = GetElementOrderImpl(i);
2412  bdofs[i+1] = nbdofs;
2413  }
2414  }
2415  else
2416  {
2417  // the simple case: all elements are the same
2418  bdofs = NULL;
2420  nbdofs = mesh->GetNE() * fec->GetNumDof(geom, order);
2421  }
2422  }
2423 
2424  ndofs = nvdofs + nedofs + nfdofs + nbdofs;
2425 
2427 
2428  // record the current mesh sequence number to detect refinement etc.
2430 
2431  // increment our sequence number to let GridFunctions know they need updating
2432  sequence++;
2433 
2434  // DOFs are now assigned according to current element orders
2435  orders_changed = false;
2436 
2437  // Do not build elem_dof Table here: in parallel it has to be constructed
2438  // later.
2439 }
2440 
2442 {
2443  MFEM_ASSERT(bits != 0, "invalid bit mask");
2444  for (int order = 0; bits != 0; order++, bits >>= 1)
2445  {
2446  if (bits & 1) { return order; }
2447  }
2448  return 0;
2449 }
2450 
2451 void FiniteElementSpace
2453  Array<VarOrderBits> &face_orders) const
2454 {
2455  MFEM_ASSERT(IsVariableOrder(), "");
2456  MFEM_ASSERT(Nonconforming(), "");
2457  MFEM_ASSERT(elem_order.Size() == mesh->GetNE(), "");
2458 
2459  edge_orders.SetSize(mesh->GetNEdges()); edge_orders = 0;
2460  face_orders.SetSize(mesh->GetNFaces()); face_orders = 0;
2461 
2462  // Calculate initial edge/face orders, as required by incident elements.
2463  // For each edge/face we accumulate in a bit-mask the orders of elements
2464  // sharing the edge/face.
2465  Array<int> E, F, ori;
2466  for (int i = 0; i < mesh->GetNE(); i++)
2467  {
2468  int order = elem_order[i];
2469  MFEM_ASSERT(order <= MaxVarOrder, "");
2470  VarOrderBits mask = (VarOrderBits(1) << order);
2471 
2472  mesh->GetElementEdges(i, E, ori);
2473  for (int j = 0; j < E.Size(); j++)
2474  {
2475  edge_orders[E[j]] |= mask;
2476  }
2477 
2478  if (mesh->Dimension() > 2)
2479  {
2480  mesh->GetElementFaces(i, F, ori);
2481  for (int j = 0; j < F.Size(); j++)
2482  {
2483  face_orders[F[j]] |= mask;
2484  }
2485  }
2486  }
2487 
2488  if (relaxed_hp)
2489  {
2490  // for relaxed conformity we don't need the masters to match the minimum
2491  // orders of the slaves, we can stop now
2492  return;
2493  }
2494 
2495  // Iterate while minimum orders propagate by master/slave relations
2496  // (and new orders also propagate from faces to incident edges).
2497  // See https://github.com/mfem/mfem/pull/1423#issuecomment-638930559
2498  // for an illustration of why this is necessary in hp meshes.
2499  bool done;
2500  do
2501  {
2502  done = true;
2503 
2504  // propagate from slave edges to master edges
2505  const NCMesh::NCList &edge_list = mesh->ncmesh->GetEdgeList();
2506  for (const NCMesh::Master &master : edge_list.masters)
2507  {
2508  VarOrderBits slave_orders = 0;
2509  for (int i = master.slaves_begin; i < master.slaves_end; i++)
2510  {
2511  slave_orders |= edge_orders[edge_list.slaves[i].index];
2512  }
2513 
2514  int min_order = MinOrder(slave_orders);
2515  if (min_order < MinOrder(edge_orders[master.index]))
2516  {
2517  edge_orders[master.index] |= (VarOrderBits(1) << min_order);
2518  done = false;
2519  }
2520  }
2521 
2522  // propagate from slave faces(+edges) to master faces(+edges)
2523  const NCMesh::NCList &face_list = mesh->ncmesh->GetFaceList();
2524  for (const NCMesh::Master &master : face_list.masters)
2525  {
2526  VarOrderBits slave_orders = 0;
2527  for (int i = master.slaves_begin; i < master.slaves_end; i++)
2528  {
2529  const NCMesh::Slave &slave = face_list.slaves[i];
2530  if (slave.index >= 0)
2531  {
2532  slave_orders |= face_orders[slave.index];
2533 
2534  mesh->GetFaceEdges(slave.index, E, ori);
2535  for (int j = 0; j < E.Size(); j++)
2536  {
2537  slave_orders |= edge_orders[E[j]];
2538  }
2539  }
2540  else
2541  {
2542  // degenerate face (i.e., edge-face constraint)
2543  slave_orders |= edge_orders[-1 - slave.index];
2544  }
2545  }
2546 
2547  int min_order = MinOrder(slave_orders);
2548  if (min_order < MinOrder(face_orders[master.index]))
2549  {
2550  face_orders[master.index] |= (VarOrderBits(1) << min_order);
2551  done = false;
2552  }
2553  }
2554 
2555  // make sure edges support (new) orders required by incident faces
2556  for (int i = 0; i < mesh->GetNFaces(); i++)
2557  {
2558  mesh->GetFaceEdges(i, E, ori);
2559  for (int j = 0; j < E.Size(); j++)
2560  {
2561  edge_orders[E[j]] |= face_orders[i];
2562  }
2563  }
2564  }
2565  while (!done);
2566 }
2567 
2569  const Array<int> &entity_orders,
2570  Table &entity_dofs,
2571  Array<char> *var_ent_order)
2572 {
2573  // The tables var_edge_dofs and var_face_dofs hold DOF assignments for edges
2574  // and faces of a variable order space, in which each edge/face may host
2575  // several DOF sets, called DOF set variants. Example: an edge 'i' shared by
2576  // 4 hexes of orders 2, 3, 4, 5 will hold four DOF sets, each starting at
2577  // indices e.g. 100, 101, 103, 106, respectively. These numbers are stored
2578  // in row 'i' of var_edge_dofs. Variant zero is always the lowest order DOF
2579  // set, followed by consecutive ranges of higher order DOFs. Variable order
2580  // faces are handled similarly by var_face_dofs, which holds at most two DOF
2581  // set variants per face. The tables are empty for constant-order spaces.
2582 
2583  int num_ent = entity_orders.Size();
2584  int total_dofs = 0;
2585 
2586  Array<Connection> list;
2587  list.Reserve(2*num_ent);
2588 
2589  if (var_ent_order)
2590  {
2591  var_ent_order->SetSize(0);
2592  var_ent_order->Reserve(num_ent);
2593  }
2594 
2595  // assign DOFs according to order bit masks
2596  for (int i = 0; i < num_ent; i++)
2597  {
2598  auto geom = (ent_dim == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
2599 
2600  VarOrderBits bits = entity_orders[i];
2601  for (int order = 0; bits != 0; order++, bits >>= 1)
2602  {
2603  if (bits & 1)
2604  {
2605  int dofs = fec->GetNumDof(geom, order);
2606  list.Append(Connection(i, total_dofs));
2607  total_dofs += dofs;
2608 
2609  if (var_ent_order) { var_ent_order->Append(order); }
2610  }
2611  }
2612  }
2613 
2614  // append a dummy row as terminator
2615  list.Append(Connection(num_ent, total_dofs));
2616 
2617  // build the table
2618  entity_dofs.MakeFromList(num_ent+1, list);
2619 
2620  return total_dofs;
2621 }
2622 
2623 int FiniteElementSpace::FindDofs(const Table &var_dof_table,
2624  int row, int ndof) const
2625 {
2626  const int *beg = var_dof_table.GetRow(row);
2627  const int *end = var_dof_table.GetRow(row + 1); // terminator, see above
2628 
2629  while (beg < end)
2630  {
2631  // return the appropriate range of DOFs
2632  if ((beg[1] - beg[0]) == ndof) { return beg[0]; }
2633  beg++;
2634  }
2635 
2636  MFEM_ABORT("DOFs not found for ndof = " << ndof);
2637  return 0;
2638 }
2639 
2640 int FiniteElementSpace::GetEdgeOrder(int edge, int variant) const
2641 {
2642  if (!IsVariableOrder()) { return fec->GetOrder(); }
2643 
2644  const int* beg = var_edge_dofs.GetRow(edge);
2645  const int* end = var_edge_dofs.GetRow(edge + 1);
2646  if (variant >= end - beg) { return -1; } // past last variant
2647 
2648  return var_edge_orders[var_edge_dofs.GetI()[edge] + variant];
2649 }
2650 
2651 int FiniteElementSpace::GetFaceOrder(int face, int variant) const
2652 {
2653  if (!IsVariableOrder())
2654  {
2655  // face order can be different from fec->GetOrder()
2656  Geometry::Type geom = mesh->GetFaceGeometry(face);
2657  return fec->FiniteElementForGeometry(geom)->GetOrder();
2658  }
2659 
2660  const int* beg = var_face_dofs.GetRow(face);
2661  const int* end = var_face_dofs.GetRow(face + 1);
2662  if (variant >= end - beg) { return -1; } // past last variant
2663 
2664  return var_face_orders[var_face_dofs.GetI()[face] + variant];
2665 }
2666 
2667 int FiniteElementSpace::GetNVariants(int entity, int index) const
2668 {
2669  MFEM_ASSERT(IsVariableOrder(), "");
2670  const Table &dof_table = (entity == 1) ? var_edge_dofs : var_face_dofs;
2671 
2672  MFEM_ASSERT(index >= 0 && index < dof_table.Size(), "");
2673  return dof_table.GetRow(index + 1) - dof_table.GetRow(index);
2674 }
2675 
2676 static const char* msg_orders_changed =
2677  "Element orders changed, you need to Update() the space first.";
2678 
2681 {
2682  MFEM_VERIFY(!orders_changed, msg_orders_changed);
2683 
2684  if (elem_dof)
2685  {
2686  elem_dof->GetRow(elem, dofs);
2687 
2688  if (DoFTrans[mesh->GetElementBaseGeometry(elem)])
2689  {
2690  Array<int> Fo;
2691  elem_fos -> GetRow (elem, Fo);
2692  DoFTrans[mesh->GetElementBaseGeometry(elem)]->SetFaceOrientations(Fo);
2693  }
2694  return DoFTrans[mesh->GetElementBaseGeometry(elem)];
2695  }
2696 
2697  Array<int> V, E, Eo, F, Fo; // TODO: LocalArray
2698 
2699  int dim = mesh->Dimension();
2700  auto geom = mesh->GetElementGeometry(elem);
2701  int order = GetElementOrderImpl(elem);
2702 
2703  int nv = fec->GetNumDof(Geometry::POINT, order);
2704  int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
2705  int nb = (dim > 0) ? fec->GetNumDof(geom, order) : 0;
2706 
2707  if (nv) { mesh->GetElementVertices(elem, V); }
2708  if (ne) { mesh->GetElementEdges(elem, E, Eo); }
2709 
2710  int nfd = 0;
2711  if (dim > 2 && fec->HasFaceDofs(geom, order))
2712  {
2713  mesh->GetElementFaces(elem, F, Fo);
2714  for (int i = 0; i < F.Size(); i++)
2715  {
2716  nfd += fec->GetNumDof(mesh->GetFaceGeometry(F[i]), order);
2717  }
2718  if (DoFTrans[mesh->GetElementBaseGeometry(elem)])
2719  {
2721  -> SetFaceOrientations(Fo);
2722  }
2723  }
2724 
2725  dofs.SetSize(0);
2726  dofs.Reserve(nv*V.Size() + ne*E.Size() + nfd + nb);
2727 
2728  if (nv) // vertex DOFs
2729  {
2730  for (int i = 0; i < V.Size(); i++)
2731  {
2732  for (int j = 0; j < nv; j++)
2733  {
2734  dofs.Append(V[i]*nv + j);
2735  }
2736  }
2737  }
2738 
2739  if (ne) // edge DOFs
2740  {
2741  for (int i = 0; i < E.Size(); i++)
2742  {
2743  int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
2744  const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
2745 
2746  for (int j = 0; j < ne; j++)
2747  {
2748  dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
2749  }
2750  }
2751  }
2752 
2753  if (nfd) // face DOFs
2754  {
2755  for (int i = 0; i < F.Size(); i++)
2756  {
2757  auto fgeom = mesh->GetFaceGeometry(F[i]);
2758  int nf = fec->GetNumDof(fgeom, order);
2759 
2760  int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F[i], nf) : F[i]*nf;
2761  const int *ind = fec->GetDofOrdering(fgeom, order, Fo[i]);
2762 
2763  for (int j = 0; j < nf; j++)
2764  {
2765  dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
2766  }
2767  }
2768  }
2769 
2770  if (nb) // interior ("bubble") DOFs
2771  {
2772  int bbase = bdofs ? bdofs[elem] : elem*nb;
2773  bbase += nvdofs + nedofs + nfdofs;
2774 
2775  for (int j = 0; j < nb; j++)
2776  {
2777  dofs.Append(bbase + j);
2778  }
2779  }
2780  return DoFTrans[mesh->GetElementBaseGeometry(elem)];
2781 }
2782 
2784 {
2785  if (i < 0 || !mesh->GetNE()) { return NULL; }
2786  MFEM_VERIFY(i < mesh->GetNE(),
2787  "Invalid element id " << i << ", maximum allowed " << mesh->GetNE()-1);
2788 
2789  const FiniteElement *FE =
2791 
2792  if (NURBSext)
2793  {
2794  NURBSext->LoadFE(i, FE);
2795  }
2796  else
2797  {
2798 #ifdef MFEM_DEBUG
2799  // consistency check: fec->GetOrder() and FE->GetOrder() should return
2800  // the same value (for standard, constant-order spaces)
2801  if (!IsVariableOrder() && FE->GetDim() > 0)
2802  {
2803  MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
2804  "internal error: " <<
2805  FE->GetOrder() << " != " << fec->GetOrder());
2806  }
2807 #endif
2808  }
2809 
2810  return FE;
2811 }
2812 
2815 {
2816  MFEM_VERIFY(!orders_changed, msg_orders_changed);
2817 
2818  if (bdr_elem_dof)
2819  {
2820  bdr_elem_dof->GetRow(bel, dofs);
2821 
2823  {
2824  Array<int> Fo;
2825  bdr_elem_fos -> GetRow (bel, Fo);
2827  SetFaceOrientations(Fo);
2828  }
2829  return DoFTrans[mesh->GetBdrElementBaseGeometry(bel)];
2830  }
2831 
2832  Array<int> V, E, Eo, Fo; // TODO: LocalArray
2833  int F, oF;
2834 
2835  int dim = mesh->Dimension();
2836  auto geom = mesh->GetBdrElementGeometry(bel);
2837  int order = fec->GetOrder();
2838 
2839  if (IsVariableOrder()) // determine order from adjacent element
2840  {
2841  int elem, info;
2842  mesh->GetBdrElementAdjacentElement(bel, elem, info);
2843  order = elem_order[elem];
2844  }
2845 
2846  int nv = fec->GetNumDof(Geometry::POINT, order);
2847  int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
2848  int nf = (dim > 2) ? fec->GetNumDof(geom, order) : 0;
2849 
2850  if (nv) { mesh->GetBdrElementVertices(bel, V); }
2851  if (ne) { mesh->GetBdrElementEdges(bel, E, Eo); }
2852  if (nf)
2853  {
2854  mesh->GetBdrElementFace(bel, &F, &oF);
2855 
2857  {
2858  Fo.Append(oF);
2860  SetFaceOrientations(Fo);
2861  }
2862  }
2863 
2864  dofs.SetSize(0);
2865  dofs.Reserve(nv*V.Size() + ne*E.Size() + nf);
2866 
2867  if (nv) // vertex DOFs
2868  {
2869  for (int i = 0; i < V.Size(); i++)
2870  {
2871  for (int j = 0; j < nv; j++)
2872  {
2873  dofs.Append(V[i]*nv + j);
2874  }
2875  }
2876  }
2877 
2878  if (ne) // edge DOFs
2879  {
2880  for (int i = 0; i < E.Size(); i++)
2881  {
2882  int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
2883  const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
2884 
2885  for (int j = 0; j < ne; j++)
2886  {
2887  dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
2888  }
2889  }
2890  }
2891 
2892  if (nf) // face DOFs
2893  {
2894  int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F, nf) : F*nf;
2895  const int *ind = fec->GetDofOrdering(geom, order, oF);
2896 
2897  for (int j = 0; j < nf; j++)
2898  {
2899  dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
2900  }
2901  }
2902 
2903  return DoFTrans[mesh->GetBdrElementBaseGeometry(bel)];
2904 }
2905 
2907  int variant) const
2908 {
2909  MFEM_VERIFY(!orders_changed, msg_orders_changed);
2910 
2911  // If face_dof is already built, use it.
2912  // If it is not and we have a NURBS space, build the face_dof and use it.
2913  if ((face_dof && variant == 0) ||
2914  (NURBSext && (BuildNURBSFaceToDofTable(), true)))
2915  {
2916  face_dof->GetRow(face, dofs);
2917  return fec->GetOrder();
2918  }
2919 
2920  int order, nf, fbase;
2921  int dim = mesh->Dimension();
2922  auto fgeom = (dim > 2) ? mesh->GetFaceGeometry(face) : Geometry::INVALID;
2923 
2924  if (var_face_dofs.Size() > 0) // variable orders or *mixed* faces
2925  {
2926  const int* beg = var_face_dofs.GetRow(face);
2927  const int* end = var_face_dofs.GetRow(face + 1);
2928  if (variant >= end - beg) { return -1; } // past last face DOFs
2929 
2930  fbase = beg[variant];
2931  nf = beg[variant+1] - fbase;
2932 
2933  order = !IsVariableOrder() ? fec->GetOrder() :
2934  var_face_orders[var_face_dofs.GetI()[face] + variant];
2935  MFEM_ASSERT(fec->GetNumDof(fgeom, order) == nf, "");
2936  }
2937  else
2938  {
2939  if (variant > 0) { return -1; }
2940  order = fec->GetOrder();
2941  nf = (dim > 2) ? fec->GetNumDof(fgeom, order) : 0;
2942  fbase = face*nf;
2943  }
2944 
2945  // for 1D, 2D and 3D faces
2946  int nv = fec->GetNumDof(Geometry::POINT, order);
2947  int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
2948 
2949  Array<int> V, E, Eo;
2950  if (nv) { mesh->GetFaceVertices(face, V); }
2951  if (ne) { mesh->GetFaceEdges(face, E, Eo); }
2952 
2953  dofs.SetSize(0);
2954  dofs.Reserve(V.Size() * nv + E.Size() * ne + nf);
2955 
2956  if (nv) // vertex DOFs
2957  {
2958  for (int i = 0; i < V.Size(); i++)
2959  {
2960  for (int j = 0; j < nv; j++)
2961  {
2962  dofs.Append(V[i]*nv + j);
2963  }
2964  }
2965  }
2966  if (ne) // edge DOFs
2967  {
2968  for (int i = 0; i < E.Size(); i++)
2969  {
2970  int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
2971  const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
2972 
2973  for (int j = 0; j < ne; j++)
2974  {
2975  dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
2976  }
2977  }
2978  }
2979  for (int j = 0; j < nf; j++)
2980  {
2981  dofs.Append(nvdofs + nedofs + fbase + j);
2982  }
2983 
2984  return order;
2985 }
2986 
2988  int variant) const
2989 {
2990  MFEM_VERIFY(!orders_changed, msg_orders_changed);
2991 
2992  int order, ne, base;
2993  if (IsVariableOrder())
2994  {
2995  const int* beg = var_edge_dofs.GetRow(edge);
2996  const int* end = var_edge_dofs.GetRow(edge + 1);
2997  if (variant >= end - beg) { return -1; } // past last edge DOFs
2998 
2999  base = beg[variant];
3000  ne = beg[variant+1] - base;
3001 
3002  order = var_edge_orders[var_edge_dofs.GetI()[edge] + variant];
3003  MFEM_ASSERT(fec->GetNumDof(Geometry::SEGMENT, order) == ne, "");
3004  }
3005  else
3006  {
3007  if (variant > 0) { return -1; }
3008  order = fec->GetOrder();
3009  ne = fec->GetNumDof(Geometry::SEGMENT, order);
3010  base = edge*ne;
3011  }
3012 
3013  Array<int> V; // TODO: LocalArray
3014  int nv = fec->GetNumDof(Geometry::POINT, order);
3015  if (nv) { mesh->GetEdgeVertices(edge, V); }
3016 
3017  dofs.SetSize(0);
3018  dofs.Reserve(2*nv + ne);
3019 
3020  for (int i = 0; i < 2; i++)
3021  {
3022  for (int j = 0; j < nv; j++)
3023  {
3024  dofs.Append(V[i]*nv + j);
3025  }
3026  }
3027  for (int j = 0; j < ne; j++)
3028  {
3029  dofs.Append(nvdofs + base + j);
3030  }
3031 
3032  return order;
3033 }
3034 
3036 {
3037  int nv = fec->DofForGeometry(Geometry::POINT);
3038  dofs.SetSize(nv);
3039  for (int j = 0; j < nv; j++)
3040  {
3041  dofs[j] = i*nv+j;
3042  }
3043 }
3044 
3046 {
3047  MFEM_VERIFY(!orders_changed, msg_orders_changed);
3048 
3050  int base = bdofs ? bdofs[i] : i*nb;
3051 
3052  dofs.SetSize(nb);
3053  base += nvdofs + nedofs + nfdofs;
3054  for (int j = 0; j < nb; j++)
3055  {
3056  dofs[j] = base + j;
3057  }
3058 }
3059 
3061 {
3062  return fec->GetNumDof(mesh->GetElementGeometry(i),
3063  GetElementOrderImpl(i));
3064 }
3065 
3067 {
3068  MFEM_VERIFY(!IsVariableOrder(), "not implemented");
3069 
3071  dofs.SetSize (ne);
3072  for (int j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
3073  {
3074  dofs[j] = k;
3075  }
3076 }
3077 
3079 {
3080  MFEM_VERIFY(!IsVariableOrder(), "not implemented");
3081 
3082  int nf, base;
3083  if (var_face_dofs.Size() > 0) // mixed faces
3084  {
3085  base = var_face_dofs.GetRow(i)[0];
3086  nf = var_face_dofs.GetRow(i)[1] - base;
3087  }
3088  else
3089  {
3090  auto geom = mesh->GetFaceGeometry(0);
3091  nf = fec->GetNumDof(geom, fec->GetOrder());
3092  base = i*nf;
3093  }
3094 
3095  dofs.SetSize(nf);
3096  for (int j = 0; j < nf; j++)
3097  {
3098  dofs[j] = nvdofs + nedofs + base + j;
3099  }
3100 }
3101 
3103 {
3104  int order = fec->GetOrder();
3105 
3106  if (IsVariableOrder()) // determine order from adjacent element
3107  {
3108  int elem, info;
3109  mesh->GetBdrElementAdjacentElement(i, elem, info);
3110  order = elem_order[elem];
3111  }
3112 
3113  const FiniteElement *BE;
3114  switch (mesh->Dimension())
3115  {
3116  case 1:
3117  BE = fec->GetFE(Geometry::POINT, order);
3118  break;
3119  case 2:
3120  BE = fec->GetFE(Geometry::SEGMENT, order);
3121  break;
3122  case 3:
3123  default:
3124  BE = fec->GetFE(mesh->GetBdrElementBaseGeometry(i), order);
3125  }
3126 
3127  if (NURBSext)
3128  {
3129  NURBSext->LoadBE(i, BE);
3130  }
3131 
3132  return BE;
3133 }
3134 
3136 {
3137  MFEM_VERIFY(!IsVariableOrder(), "not implemented");
3138 
3139  const FiniteElement *fe;
3140  switch (mesh->Dimension())
3141  {
3142  case 1:
3144  break;
3145  case 2:
3147  break;
3148  case 3:
3149  default:
3151  }
3152 
3153  if (NURBSext)
3154  {
3155  // Ensure 'face_to_be' is built:
3156  if (!face_dof) { BuildNURBSFaceToDofTable(); }
3157  MFEM_ASSERT(face_to_be[i] >= 0,
3158  "NURBS mesh: only boundary faces are supported!");
3159  NURBSext->LoadBE(face_to_be[i], fe);
3160  }
3161 
3162  return fe;
3163 }
3164 
3166  int variant) const
3167 {
3168  MFEM_ASSERT(mesh->Dimension() > 1, "No edges with mesh dimension < 2");
3169 
3170  int eo = IsVariableOrder() ? GetEdgeOrder(i, variant) : fec->GetOrder();
3171  return fec->GetFE(Geometry::SEGMENT, eo);
3172 }
3173 
3175 ::GetTraceElement(int i, Geometry::Type geom_type) const
3176 {
3177  return fec->TraceFiniteElementForGeometry(geom_type);
3178 }
3179 
3181 {
3182  Destroy();
3183 }
3184 
3186 {
3187  delete cR;
3188  delete cR_hp;
3189  delete cP;
3190  Th.Clear();
3191  L2E_nat.Clear();
3192  L2E_lex.Clear();
3193  for (int i = 0; i < E2Q_array.Size(); i++)
3194  {
3195  delete E2Q_array[i];
3196  }
3197  E2Q_array.SetSize(0);
3198  for (auto &x : L2F)
3199  {
3200  delete x.second;
3201  }
3202  for (int i = 0; i < E2IFQ_array.Size(); i++)
3203  {
3204  delete E2IFQ_array[i];
3205  }
3206  E2IFQ_array.SetSize(0);
3207  for (int i = 0; i < E2BFQ_array.Size(); i++)
3208  {
3209  delete E2BFQ_array[i];
3210  }
3211  E2BFQ_array.SetSize(0);
3212 
3213  DestroyDoFTrans();
3214 
3217 
3218  if (NURBSext)
3219  {
3220  if (own_ext) { delete NURBSext; }
3221  delete face_dof;
3223  }
3224  else
3225  {
3226  delete elem_dof;
3227  delete elem_fos;
3228  delete bdr_elem_dof;
3229  delete bdr_elem_fos;
3230  delete face_dof;
3231 
3232  delete [] bdofs;
3233  }
3234  ceed::RemoveBasisAndRestriction(this);
3235 }
3236 
3238 {
3239  for (int i = 0; i < DoFTrans.Size(); i++)
3240  {
3241  delete DoFTrans[i];
3242  }
3243  DoFTrans.SetSize(0);
3244 }
3245 
3247  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
3248 {
3249  // Assumptions: see the declaration of the method.
3250 
3251  if (T.Type() == Operator::MFEM_SPARSEMAT)
3252  {
3253  Mesh::GeometryList elem_geoms(*mesh);
3254 
3256  for (int i = 0; i < elem_geoms.Size(); i++)
3257  {
3258  GetLocalRefinementMatrices(coarse_fes, elem_geoms[i],
3259  localP[elem_geoms[i]]);
3260  }
3261  T.Reset(RefinementMatrix_main(coarse_fes.GetNDofs(),
3262  coarse_fes.GetElementToDofTable(),
3263  coarse_fes.
3265  localP));
3266  }
3267  else
3268  {
3269  T.Reset(new RefinementOperator(this, &coarse_fes));
3270  }
3271 }
3272 
3274  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
3275 {
3276  const SparseMatrix *coarse_P = coarse_fes.GetConformingProlongation();
3277 
3278  Operator::Type req_type = T.Type();
3279  GetTransferOperator(coarse_fes, T);
3280 
3281  if (req_type == Operator::MFEM_SPARSEMAT)
3282  {
3284  {
3285  T.Reset(mfem::Mult(*cR, *T.As<SparseMatrix>()));
3286  }
3287  if (coarse_P)
3288  {
3289  T.Reset(mfem::Mult(*T.As<SparseMatrix>(), *coarse_P));
3290  }
3291  }
3292  else
3293  {
3294  const int RP_case = bool(GetConformingRestriction()) + 2*bool(coarse_P);
3295  if (RP_case == 0) { return; }
3296  const bool owner = T.OwnsOperator();
3297  T.SetOperatorOwner(false);
3298  switch (RP_case)
3299  {
3300  case 1:
3301  T.Reset(new ProductOperator(cR, T.Ptr(), false, owner));
3302  break;
3303  case 2:
3304  T.Reset(new ProductOperator(T.Ptr(), coarse_P, owner, false));
3305  break;
3306  case 3:
3308  cR, T.Ptr(), coarse_P, false, owner, false));
3309  break;
3310  }
3311  }
3312 }
3313 
3315 {
3317 
3318  Array<char> new_order(mesh->GetNE());
3319  switch (mesh->GetLastOperation())
3320  {
3321  case Mesh::REFINE:
3322  {
3323  for (int i = 0; i < mesh->GetNE(); i++)
3324  {
3325  new_order[i] = elem_order[cf_tr.embeddings[i].parent];
3326  }
3327  break;
3328  }
3329  default:
3330  MFEM_ABORT("not implemented yet");
3331  }
3332 
3333  mfem::Swap(elem_order, new_order);
3334 }
3335 
3336 void FiniteElementSpace::Update(bool want_transform)
3337 {
3338  if (!orders_changed)
3339  {
3340  if (mesh->GetSequence() == mesh_sequence)
3341  {
3342  return; // mesh and space are in sync, no-op
3343  }
3344  if (want_transform && mesh->GetSequence() != mesh_sequence + 1)
3345  {
3346  MFEM_ABORT("Error in update sequence. Space needs to be updated after "
3347  "each mesh modification.");
3348  }
3349  }
3350  else
3351  {
3352  if (mesh->GetSequence() != mesh_sequence)
3353  {
3354  MFEM_ABORT("Updating space after both mesh change and element order "
3355  "change is not supported. Please update separately after "
3356  "each change.");
3357  }
3358  }
3359 
3360  if (NURBSext)
3361  {
3362  UpdateNURBS();
3363  return;
3364  }
3365 
3366  Table* old_elem_dof = NULL;
3367  Table* old_elem_fos = NULL;
3368  int old_ndofs;
3369  bool old_orders_changed = orders_changed;
3370 
3371  // save old DOF table
3372  if (want_transform)
3373  {
3374  old_elem_dof = elem_dof;
3375  old_elem_fos = elem_fos;
3376  elem_dof = NULL;
3377  elem_fos = NULL;
3378  old_ndofs = ndofs;
3379  }
3380 
3381  // update the 'elem_order' array if the mesh has changed
3383  {
3385  }
3386 
3387  Destroy(); // calls Th.Clear()
3388  Construct();
3390 
3391  if (want_transform)
3392  {
3393  MFEM_VERIFY(!old_orders_changed, "Interpolation for element order change "
3394  "is not implemented yet, sorry.");
3395 
3396  // calculate appropriate GridFunction transformation
3397  switch (mesh->GetLastOperation())
3398  {
3399  case Mesh::REFINE:
3400  {
3401  if (Th.Type() != Operator::MFEM_SPARSEMAT)
3402  {
3403  Th.Reset(new RefinementOperator(this, old_elem_dof,
3404  old_elem_fos, old_ndofs));
3405  // The RefinementOperator takes ownership of 'old_elem_dof', so
3406  // we no longer own it:
3407  old_elem_dof = NULL;
3408  old_elem_fos = NULL;
3409  }
3410  else
3411  {
3412  // calculate fully assembled matrix
3413  Th.Reset(RefinementMatrix(old_ndofs, old_elem_dof,
3414  old_elem_fos));
3415  }
3416  break;
3417  }
3418 
3419  case Mesh::DEREFINE:
3420  {
3422  Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof, old_elem_fos));
3423  if (cP && cR)
3424  {
3425  Th.SetOperatorOwner(false);
3427  false, false, true));
3428  }
3429  break;
3430  }
3431 
3432  default:
3433  break;
3434  }
3435 
3436  delete old_elem_dof;
3437  delete old_elem_fos;
3438  }
3439 }
3440 
3442 {
3443  mesh = new_mesh;
3444 }
3445 
3446 void FiniteElementSpace::Save(std::ostream &os) const
3447 {
3448  int fes_format = 90; // the original format, v0.9
3449  bool nurbs_unit_weights = false;
3450 
3451  // Determine the format that should be used.
3452  if (!NURBSext)
3453  {
3454  // TODO: if this is a variable-order FE space, use fes_format = 100.
3455  }
3456  else
3457  {
3458  const NURBSFECollection *nurbs_fec =
3459  dynamic_cast<const NURBSFECollection *>(fec);
3460  MFEM_VERIFY(nurbs_fec, "invalid FE collection");
3461  nurbs_fec->SetOrder(NURBSext->GetOrder());
3462  const double eps = 5e-14;
3463  nurbs_unit_weights = (NURBSext->GetWeights().Min() >= 1.0-eps &&
3464  NURBSext->GetWeights().Max() <= 1.0+eps);
3466  (NURBSext != mesh->NURBSext && !nurbs_unit_weights) ||
3467  (NURBSext->GetMaster().Size() != 0 ))
3468  {
3469  fes_format = 100; // v1.0 format
3470  }
3471  }
3472 
3473  os << (fes_format == 90 ?
3474  "FiniteElementSpace\n" : "MFEM FiniteElementSpace v1.0\n")
3475  << "FiniteElementCollection: " << fec->Name() << '\n'
3476  << "VDim: " << vdim << '\n'
3477  << "Ordering: " << ordering << '\n';
3478 
3479  if (fes_format == 100) // v1.0
3480  {
3481  if (!NURBSext)
3482  {
3483  // TODO: this is a variable-order FE space --> write 'element_orders'.
3484  }
3485  else if (NURBSext != mesh->NURBSext)
3486  {
3488  {
3489  os << "NURBS_order\n" << NURBSext->GetOrder() << '\n';
3490  }
3491  else
3492  {
3493  os << "NURBS_orders\n";
3494  // 1 = do not write the size, just the entries:
3495  NURBSext->GetOrders().Save(os, 1);
3496  }
3497  // If periodic BCs are given, write connectivity
3498  if (NURBSext->GetMaster().Size() != 0 )
3499  {
3500  os <<"NURBS_periodic\n";
3501  NURBSext->GetMaster().Save(os);
3502  NURBSext->GetSlave().Save(os);
3503  }
3504  // If the weights are not unit, write them to the output:
3505  if (!nurbs_unit_weights)
3506  {
3507  os << "NURBS_weights\n";
3508  NURBSext->GetWeights().Print(os, 1);
3509  }
3510  }
3511  os << "End: MFEM FiniteElementSpace v1.0\n";
3512  }
3513 }
3514 
3516 {
3517  string buff;
3518  int fes_format = 0, ord;
3519  FiniteElementCollection *r_fec;
3520 
3521  Destroy();
3522 
3523  input >> std::ws;
3524  getline(input, buff); // 'FiniteElementSpace'
3525  filter_dos(buff);
3526  if (buff == "FiniteElementSpace") { fes_format = 90; /* v0.9 */ }
3527  else if (buff == "MFEM FiniteElementSpace v1.0") { fes_format = 100; }
3528  else { MFEM_ABORT("input stream is not a FiniteElementSpace!"); }
3529  getline(input, buff, ' '); // 'FiniteElementCollection:'
3530  input >> std::ws;
3531  getline(input, buff);
3532  filter_dos(buff);
3533  r_fec = FiniteElementCollection::New(buff.c_str());
3534  getline(input, buff, ' '); // 'VDim:'
3535  input >> vdim;
3536  getline(input, buff, ' '); // 'Ordering:'
3537  input >> ord;
3538 
3539  NURBSFECollection *nurbs_fec = dynamic_cast<NURBSFECollection*>(r_fec);
3540  NURBSExtension *nurbs_ext = NULL;
3541  if (fes_format == 90) // original format, v0.9
3542  {
3543  if (nurbs_fec)
3544  {
3545  MFEM_VERIFY(m->NURBSext, "NURBS FE collection requires a NURBS mesh!");
3546  const int order = nurbs_fec->GetOrder();
3547  if (order != m->NURBSext->GetOrder() &&
3549  {
3550  nurbs_ext = new NURBSExtension(m->NURBSext, order);
3551  }
3552  }
3553  }
3554  else if (fes_format == 100) // v1.0
3555  {
3556  while (1)
3557  {
3558  skip_comment_lines(input, '#');
3559  MFEM_VERIFY(input.good(), "error reading FiniteElementSpace v1.0");
3560  getline(input, buff);
3561  filter_dos(buff);
3562  if (buff == "NURBS_order" || buff == "NURBS_orders")
3563  {
3564  MFEM_VERIFY(nurbs_fec,
3565  buff << ": NURBS FE collection is required!");
3566  MFEM_VERIFY(m->NURBSext, buff << ": NURBS mesh is required!");
3567  MFEM_VERIFY(!nurbs_ext, buff << ": order redefinition!");
3568  if (buff == "NURBS_order")
3569  {
3570  int order;
3571  input >> order;
3572  nurbs_ext = new NURBSExtension(m->NURBSext, order);
3573  }
3574  else
3575  {
3576  Array<int> orders;
3577  orders.Load(m->NURBSext->GetNKV(), input);
3578  nurbs_ext = new NURBSExtension(m->NURBSext, orders);
3579  }
3580  }
3581  else if (buff == "NURBS_periodic")
3582  {
3583  Array<int> master, slave;
3584  master.Load(input);
3585  slave.Load(input);
3586  nurbs_ext->ConnectBoundaries(master,slave);
3587  }
3588  else if (buff == "NURBS_weights")
3589  {
3590  MFEM_VERIFY(nurbs_ext, "NURBS_weights: NURBS_orders have to be "
3591  "specified before NURBS_weights!");
3592  nurbs_ext->GetWeights().Load(input, nurbs_ext->GetNDof());
3593  }
3594  else if (buff == "element_orders")
3595  {
3596  MFEM_VERIFY(!nurbs_fec, "section element_orders cannot be used "
3597  "with a NURBS FE collection");
3598  MFEM_ABORT("element_orders: not implemented yet!");
3599  }
3600  else if (buff == "End: MFEM FiniteElementSpace v1.0")
3601  {
3602  break;
3603  }
3604  else
3605  {
3606  MFEM_ABORT("unknown section: " << buff);
3607  }
3608  }
3609  }
3610 
3611  Constructor(m, nurbs_ext, r_fec, vdim, ord);
3612 
3613  return r_fec;
3614 }
3615 
3616 
3618 {
3619  // protected method
3620  int offset = 0;
3621  const int num_elem = mesh->GetNE();
3622  element_offsets = new int[num_elem + 1];
3623  for (int g = 0; g < Geometry::NumGeom; g++)
3624  {
3625  int_rule[g] = NULL;
3626  }
3627  for (int i = 0; i < num_elem; i++)
3628  {
3629  element_offsets[i] = offset;
3630  int geom = mesh->GetElementBaseGeometry(i);
3631  if (int_rule[geom] == NULL)
3632  {
3633  int_rule[geom] = &IntRules.Get(geom, order);
3634  }
3635  offset += int_rule[geom]->GetNPoints();
3636  }
3637  element_offsets[num_elem] = size = offset;
3638 }
3639 
3640 QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
3641  : mesh(mesh_)
3642 {
3643  const char *msg = "invalid input stream";
3644  string ident;
3645 
3646  in >> ident; MFEM_VERIFY(ident == "QuadratureSpace", msg);
3647  in >> ident; MFEM_VERIFY(ident == "Type:", msg);
3648  in >> ident;
3649  if (ident == "default_quadrature")
3650  {
3651  in >> ident; MFEM_VERIFY(ident == "Order:", msg);
3652  in >> order;
3653  }
3654  else
3655  {
3656  MFEM_ABORT("unknown QuadratureSpace type: " << ident);
3657  return;
3658  }
3659 
3660  Construct();
3661 }
3662 
3663 void QuadratureSpace::Save(std::ostream &os) const
3664 {
3665  os << "QuadratureSpace\n"
3666  << "Type: default_quadrature\n"
3667  << "Order: " << order << '\n';
3668 }
3669 
3670 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof, const Table *old_elem_fos)
Definition: fespace.cpp:1476
Abstract class for all finite elements.
Definition: fe_base.hpp:235
int GetOrder() const
Get the order of the NURBS collection: either a positive number, when using fixed order...
Definition: fe_coll.hpp:634
VDofTransformation VDoFTrans
Definition: fespace.hpp:144
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:602
Field is discontinuous across element interfaces.
Definition: fe_coll.hpp:48
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:575
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:336
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:560
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:563
void GetVertexVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:324
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
const Vector & GetWeights() const
Definition: nurbs.hpp:386
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
Definition: fespace.cpp:1389
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
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:560
const FiniteElement * GetFE(Geometry::Type geom, int p) const
Variable order version of FiniteElementForGeometry().
Definition: fe_coll.hpp:161
int * GetJ()
Definition: table.hpp:114
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
Definition: fespace.hpp:111
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
int * bdofs
internal DOFs of elements if mixed/var-order; NULL otherwise
Definition: fespace.hpp:119
const Array< int > & GetMaster() const
Definition: nurbs.hpp:329
void GetElementInteriorVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:330
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:251
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:5923
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1234
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I, int skipfirst=0)
Definition: fespace.cpp:758
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:4399
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
void SetSize(int dim, int connections_per_row)
Set the size and the number of connections for the table.
Definition: table.cpp:124
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:455
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:2934
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:5985
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:2040
Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
void BuildElementToDofTable() const
Definition: fespace.cpp:342
void AddColumnsInRow(int r, int ncol)
Definition: table.hpp:78
int MakeDofTable(int ent_dim, const Array< int > &entity_orders, Table &entity_dofs, Array< char > *var_ent_order)
Definition: fespace.cpp:2568
void BuildNURBSFaceToDofTable() const
Generates partial face_dof table for a NURBS space.
Definition: fespace.cpp:2266
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th edge in the ...
Definition: fespace.cpp:3165
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition: table.cpp:81
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:125
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
DerefinementOperator(const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes, BilinearFormIntegrator *mass_integ)
TODO: Implement DofTransformation support.
Definition: fespace.cpp:1881
static const int NumGeom
Definition: geom.hpp:42
Array< Slave > slaves
Definition: ncmesh.hpp:231
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element &#39;i&#39; to the given &#39;geom_type&#39;.
Definition: fespace.cpp:3175
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:515
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:6140
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:5877
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:931
Ordering::Type ordering
Definition: fespace.hpp:108
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element &#39;bel&#39;.
Definition: fespace.cpp:2814
void SetElementOrder(int i, int p)
Sets the order of the i&#39;th finite element.
Definition: fespace.cpp:151
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:383
static int DecodeDof(int dof)
Helpers to remove encoded sign from a DOF.
Definition: fespace.hpp:254
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:308
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:933
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:227
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int SizeK() const
Definition: densemat.hpp:834
static int MinOrder(VarOrderBits bits)
Return the minimum order (least significant bit set) in the bit mask.
Definition: fespace.cpp:2441
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:99
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
static constexpr int MaxVarOrder
Definition: fespace.hpp:213
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
unsigned matrix
index into NCList::point_matrices[geom]
Definition: ncmesh.hpp:217
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
int GetNKV() const
Definition: nurbs.hpp:355
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
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:495
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition: mesh.hpp:1120
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int vdim
Vector dimension (number of unknowns per degree of freedom).
Definition: fespace.hpp:103
Geometry::Type Geom() const
Definition: ncmesh.hpp:194
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
OperatorHandle L2E_lex
Definition: fespace.hpp:160
int GetNumElementInteriorDofs(int i) const
Definition: fespace.cpp:3060
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
int GetElementOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:178
int GetBdrAttribute(int i) const
Definition: fespace.hpp:637
Operator that extracts face degrees of freedom for L2 nonconforming spaces.
int GetBdrElementEdgeIndex(int i) const
Definition: mesh.cpp:6163
const SparseMatrix * GetHpConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1248
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
const Element * GetFace(int i) const
Definition: mesh.hpp:1041
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:2987
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:106
T Sum()
Return the sum of all the array entries using the &#39;+&#39;&#39; operator for class &#39;T&#39;.
Definition: array.cpp:115
MFEM_DEPRECATED void RebuildElementToDofTable()
(
Definition: fespace.cpp:431
Geometry::Type GetFaceBaseGeometry(int i) const
Definition: mesh.hpp:1059
bool operator<(const Pair< A, B > &p, const Pair< A, B > &q)
Comparison operator for class Pair, based on the first element only.
Definition: sort_pairs.hpp:36
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:1661
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th face in the ...
Definition: fespace.cpp:3135
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, defined on a coarse mesh, to this FE space, defined on a refined mesh.
Definition: fespace.cpp:3273
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:783
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
Definition: fespace.cpp:3515
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2884
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:397
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes, const Array< int > *perm)
Copies the prolongation and restriction matrices from fes.
Definition: fespace.cpp:100
const FiniteElementCollection * fec
Associated FE collection (not owned).
Definition: fespace.hpp:100
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
Definition: fespace.hpp:157
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:2906
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1062
int Size_of_connections() const
Definition: table.hpp:98
void GetLocalDerefinementMatrices(Geometry::Type geom, DenseTensor &localR) const
Definition: fespace.cpp:2016
Geometry::Type GetBdrElementGeometry(int i) const
Definition: mesh.hpp:1053
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
void DeleteAll()
Delete the whole array.
Definition: array.hpp:846
int GetNDof() const
Definition: nurbs.hpp:365
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:963
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:195
void AddConnections(int r, const int *c, int nc)
Definition: table.cpp:104
std::uint64_t VarOrderBits
Bit-mask representing a set of orders needed by an edge/face.
Definition: fespace.hpp:212
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Definition: ncmesh.hpp:253
Array< int > dof_elem_array
Definition: fespace.hpp:137
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:598
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:297
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:3078
virtual int DofForGeometry(Geometry::Type GeomType) const =0
virtual void UpdateMeshPointer(Mesh *new_mesh)
Definition: fespace.cpp:3441
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:638
const IntegrationRule * IntRule
Not owned.
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1261
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5345
int uni_fdof
of single face DOFs if all faces uniform; -1 otherwise
Definition: fespace.hpp:118
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1241
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
ID for class SparseMatrix.
Definition: operator.hpp:258
const Array< int > & GetSlave() const
Definition: nurbs.hpp:331
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: fespace.cpp:1983
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:57
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:590
void Set(const int i, const int j, const double val)
Definition: sparsemat.cpp:2630
long GetSequence() const
Definition: mesh.hpp:1639
void GetVertexDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:3035
void SetPointMat(const DenseMatrix &pm)
Set the underlying point matrix describing the transformation.
Definition: eltrans.hpp:401
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:188
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:743
FaceType
Definition: mesh.hpp:45
void BuildBdrElementToDofTable() const
Definition: fespace.cpp:383
int HasFaceDofs(Geometry::Type geom, int p) const
Definition: fe_coll.cpp:25
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:257
int FindFaceDof(int face, int ndof) const
Similar to FindEdgeDof, but used for mixed meshes too.
Definition: fespace.hpp:240
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
int GetFaceOrder(int face, int variant=0) const
Returns the polynomial degree of the i&#39;th face finite element.
Definition: fespace.cpp:2651
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Extract all column indices and values from a given row.
Definition: sparsemat.cpp:2797
Array< DofTransformation * > DoFTrans
Definition: fespace.hpp:143
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: fespace.cpp:1588
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2042
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:599
void MakeFromList(int nrows, const Array< Connection > &list)
Definition: table.cpp:276
bool orders_changed
True if at least one element order changed (variable-order space only).
Definition: fespace.hpp:189
DenseTensor point_matrices[Geometry::NumGeom]
Definition: ncmesh.hpp:77
Data type sparse matrix.
Definition: sparsemat.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int FindEdgeDof(int edge, int ndof) const
Definition: fespace.hpp:236
unsigned matrix
Definition: ncmesh.hpp:58
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
void EnsureMultTranspose() const
Ensures that the matrix is capable of performing MultTranspose(), AddMultTranspose(), and AbsMultTranspose().
Definition: sparsemat.cpp:881
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
Array< char > elem_order
Definition: fespace.hpp:115
Operator that extracts Face degrees of freedom for L2 spaces.
void GetVDofs(int vd, Array< int > &dofs, int ndofs=-1) const
Returns the indices of all of the VDofs for the specified dimension &#39;vd&#39;.
Definition: fespace.cpp:195
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
void SetDofTransformation(DofTransformation &doftrans)
Set or change the nested DofTransformation object.
Definition: doftrans.hpp:175
double b
Definition: lissajous.cpp:42
Type
Ordering methods:
Definition: fespace.hpp:33
SparseMatrix * cR
Conforming restriction matrix such that cR.cP=I.
Definition: fespace.hpp:151
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1254
const IntegrationRule * IntRule
Not owned.
SparseMatrix * D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Generate the global restriction matrix from a discontinuous FE space to the piecewise constant FE spa...
Definition: fespace.cpp:679
static const int NumVerts[NumGeom]
Definition: geom.hpp:49
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:647
A pair of objects.
Definition: sort_pairs.hpp:23
void GetBoundaryTrueDofs(Array< int > &boundary_dofs, int component=-1)
Get a list of all boundary true dofs, boundary_dofs. For spaces with &#39;vdim&#39; &gt; 1, the &#39;component&#39; param...
Definition: fespace.cpp:582
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void BuildConformingInterpolation() const
Calculate the cP and cR matrices for a nonconforming mesh.
Definition: fespace.cpp:928
Geometry::Type GetFaceGeometry(int i) const
Definition: mesh.hpp:1043
NURBSExtension * StealNURBSext()
Definition: fespace.cpp:2231
Array< DenseMatrix * > point_matrices[Geometry::NumGeom]
List of unique point matrices for each slave geometry.
Definition: ncmesh.hpp:234
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:156
Array< Embedding > embeddings
Fine element positions in their parents.
Definition: ncmesh.hpp:73
virtual const FiniteElement * TraceFiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.hpp:67
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:6175
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Definition: fespace.hpp:901
Array< int > dof_ldof_array
Definition: fespace.hpp:137
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:838
const Array< int > & GetOrders() const
Read-only access to the orders of all knot vectors.
Definition: nurbs.hpp:350
bool Conforming() const
Definition: fespace.hpp:439
SparseMatrix * cP
Definition: fespace.hpp:149
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:557
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
const int * GetDofOrdering(Geometry::Type geom, int p, int ori) const
Variable order version of DofOrderForOrientation().
Definition: fe_coll.hpp:185
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition: fespace.hpp:596
int Dimension() const
Definition: mesh.hpp:999
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:250
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Definition: fespace.cpp:1294
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:4964
Array< QuadratureInterpolator * > E2Q_array
Definition: fespace.hpp:176
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:626
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
int GetNumBorderDofs(Geometry::Type geom, int order) const
Definition: fespace.cpp:894
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
void Add(const int i, const int j, const double val)
Definition: sparsemat.cpp:2649
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:903
OperatorHandle L2E_nat
The element restriction operators, see GetElementRestriction().
Definition: fespace.hpp:160
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2680
int GetNumDof(Geometry::Type geom, int p) const
Variable order version of DofForGeometry().
Definition: fe_coll.hpp:173
int slaves_end
slave faces
Definition: ncmesh.hpp:205
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:3246
int SizeI() const
Definition: densemat.hpp:832
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:3045
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
Definition: fespace.cpp:1453
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:353
List of mesh geometries stored as Array&lt;Geometry::Type&gt;.
Definition: mesh.hpp:1085
DoF transformation implementation for the Nedelec basis on triangles.
Definition: doftrans.hpp:232
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1182
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:255
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:2913
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:3363
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:383
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:27
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:711
Nonconforming edge/face within a bigger edge/face.
Definition: ncmesh.hpp:214
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
virtual ~FiniteElementSpace()
Definition: fespace.cpp:3180
bool OwnsOperator() const
Return true if the OperatorHandle owns the held Operator.
Definition: handle.hpp:117
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Given a particular Finite Element computes the element matrix elmat.
Definition: bilininteg.cpp:123
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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:710
virtual void TransformPrimal(double *v) const =0
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:410
RefinementOperator(const FiniteElementSpace *fespace, Table *old_elem_dof, Table *old_elem_fos, int old_ndofs)
Definition: fespace.cpp:1496
void SetVDim(int vdim)
Set or change the vdim parameter.
Definition: doftrans.hpp:162
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Definition: fespace.cpp:3663
Array< char > var_face_orders
Definition: fespace.hpp:128
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:124
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition: table.hpp:27
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3252
Geometry::Type GetElementGeometry(int i) const
Definition: mesh.hpp:1048
Table * GetElementDofTable()
Definition: nurbs.hpp:376
void MakeVDimMatrix(SparseMatrix &mat) const
Replicate &#39;mat&#39; in the vector dimension, according to vdim ordering mode.
Definition: fespace.cpp:1206
Table * GetBdrElementDofTable()
Definition: nurbs.hpp:377
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:29
void GetEdgeInteriorVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:336
void ShiftUpI()
Definition: table.cpp:115
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
int GetNVariants(int entity, int index) const
Return number of possible DOF variants for edge/face (var. order spaces).
Definition: fespace.cpp:2667
const CoarseFineTransformations & GetRefinementTransforms()
Definition: mesh.cpp:9849
Operation GetLastOperation() const
Return type of last modification of the mesh.
Definition: mesh.hpp:1633
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:51
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
Definition: array.hpp:316
bool Nonconforming() const
Definition: fespace.hpp:440
int SizeJ() const
Definition: densemat.hpp:833
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:925
int GetDegenerateFaceDofs(int index, Array< int > &dofs, Geometry::Type master_geom, int variant) const
Definition: fespace.cpp:852
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:630
virtual const char * Name() const
Definition: fe_coll.hpp:61
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
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:5901
Mesh * mesh
The mesh that FE space lives on (not owned).
Definition: fespace.hpp:97
void ReorderElementToDofTable()
Reorder the scalar DOFs based on the element ordering.
Definition: fespace.cpp:440
GridFunction interpolation operator applicable after mesh refinement.
Definition: fespace.hpp:290
void CalcEdgeFaceVarOrders(Array< VarOrderBits > &edge_orders, Array< VarOrderBits > &face_orders) const
Definition: fespace.cpp:2452
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:882
const NCList & GetNCList(int entity)
Return vertex/edge/face list (entity = 0/1/2, respectively).
Definition: ncmesh.hpp:275
Array< FaceQuadratureInterpolator * > E2BFQ_array
Definition: fespace.hpp:178
const Table * GetElementToFaceOrientationTable() const
Definition: fespace.hpp:722
int GetNConformingDofs() const
Definition: fespace.cpp:1255
void ConnectBoundaries()
Definition: nurbs.cpp:1793
Table var_face_dofs
NOTE: also used for spaces with mixed faces.
Definition: fespace.hpp:124
Array< Master > masters
Definition: ncmesh.hpp:230
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
Definition: mesh.cpp:5955
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:3066
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition: fespace.cpp:1329
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:66
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:807
void MakeJ()
Definition: table.cpp:91
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:112
int dim
Definition: ex24.cpp:53
Array< FaceQuadratureInterpolator * > E2IFQ_array
Definition: fespace.hpp:177
void BuildFaceToDofTable() const
Definition: fespace.cpp:405
void UpdateElementOrders()
Resize the elem_order array on mesh change.
Definition: fespace.cpp:3314
void SetFE(const FiniteElement *FE)
Set the element that will be used to compute the transformations.
Definition: eltrans.hpp:381
void OrientedPointMatrix(const Slave &slave, DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
Definition: ncmesh.cpp:3224
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
int index(int i, int j, int nx, int ny)
Definition: life.cpp:237
Array< char > var_edge_orders
Definition: fespace.hpp:128
int GetEdgeOrder(int edge, int variant=0) const
Definition: fespace.cpp:2640
Lexicographic ordering for tensor-product FiniteElements.
virtual int GetNVertices() const =0
int parent
Coarse Element index in the coarse mesh.
Definition: ncmesh.hpp:52
FiniteElementSpace()
Default constructor: the object is invalid until initialized using the method Load().
Definition: fespace.cpp:59
QuadratureSpace(Mesh *mesh_, int order_)
Create a QuadratureSpace based on the global rules from IntRules.
Definition: fespace.hpp:955
void Constructor(Mesh *mesh, NURBSExtension *ext, const FiniteElementCollection *fec, int vdim=1, int ordering=Ordering::byNODES)
Help function for constructors + Load().
Definition: fespace.cpp:2141
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:273
Geometry::Type GetBdrElementBaseGeometry(int i) const
Definition: mesh.hpp:1065
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th face element (2D and 3D).
Definition: fespace.cpp:312
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2837
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th edge.
Definition: fespace.cpp:318
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:617
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:934
void filter_dos(std::string &line)
Check for, and remove, a trailing &#39;\r&#39; from and std::string.
Definition: text.hpp:45
MFEM_HOST_DEVICE void Set(const int height, const int width, const double alpha, const TA *Adata, TB *Bdata)
Compute B = alpha*A, where the matrices A and B are of size height x width with data Adata and Bdata...
Definition: kernels.hpp:326
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:864
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:937
BiLinear2DFiniteElement QuadrilateralFE
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2633
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:577
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:160
std::tuple< bool, ElementDofOrdering, FaceType, L2FaceValues > key_face
The face restriction operators, see GetFaceRestriction().
Definition: fespace.hpp:162
Array< int > face_to_be
Definition: fespace.hpp:141
int * GetI()
Definition: table.hpp:113
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:726
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:70
RefCoord s[3]
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:118
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition: ncmesh.hpp:260
SparseMatrix * cR_hp
A version of the conforming restriction matrix for variable-order spaces.
Definition: fespace.hpp:153
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Definition: fespace.cpp:838
Class representing the storage layout of a QuadratureFunction.
Definition: fespace.hpp:935
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
int RowSize(int i) const
Definition: table.hpp:108
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3102
NURBSExtension * NURBSext
Definition: fespace.hpp:139
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Base class for operators that extracts Face degrees of freedom.
Abstract operator.
Definition: operator.hpp:24
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:404
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:4044
int index
Mesh number.
Definition: ncmesh.hpp:189
DoF transformation implementation for the Nedelec basis on tetrahedra.
Definition: doftrans.hpp:252
int GetElementOrderImpl(int i) const
Return element order: internal version of GetElementOrder without checks.
Definition: fespace.cpp:189
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:786
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
void Save(std::ostream &out) const
Save finite element space to output stream out.
Definition: fespace.cpp:3446
Abstract data type element.
Definition: element.hpp:28
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
Definition: fespace.cpp:460
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:1106
Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:84
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1102
DofTransformation * GetDofTransformation() const
Return the nested DofTransformation object.
Definition: doftrans.hpp:182
int FindDofs(const Table &var_dof_table, int row, int ndof) const
Search row of a DOF table for a DOF set of size &#39;ndof&#39;, return first DOF.
Definition: fespace.cpp:2623
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:267
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:132
DoF transformation implementation for the Nedelec basis on wedge elements.
Definition: doftrans.hpp:272
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1037
Defines the position of a fine element within a coarse element.
Definition: ncmesh.hpp:49
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
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:6117
const QuadratureSpace * qspace
Not owned.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:215
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:1358