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