MFEM  v4.1.0 Finite element discretization library
fespace.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 // Implementation of FiniteElementSpace
13
14 #include "../general/text.hpp"
15 #include "../general/forall.hpp"
17 #include "fem.hpp"
18
19 #include <cmath>
20 #include <cstdarg>
21 #include <limits>
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  fdofs(NULL), bdofs(NULL),
63  elem_dof(NULL), bdrElem_dof(NULL),
64  NURBSext(NULL), own_ext(false),
65  cP(NULL), cR(NULL), cP_is_set(false),
66  Th(Operator::ANY_TYPE),
67  sequence(0)
68 { }
69
71  Mesh *mesh,
72  const FiniteElementCollection *fec)
73 {
74  mesh = mesh ? mesh : orig.mesh;
75  fec = fec ? fec : orig.fec;
76  NURBSExtension *NURBSext = NULL;
77  if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
78  {
79 #ifdef MFEM_USE_MPI
80  ParNURBSExtension *pNURBSext =
81  dynamic_cast<ParNURBSExtension *>(orig.NURBSext);
82  if (pNURBSext)
83  {
84  NURBSext = new ParNURBSExtension(*pNURBSext);
85  }
86  else
87 #endif
88  {
89  NURBSext = new NURBSExtension(*orig.NURBSext);
90  }
91  }
92  Constructor(mesh, NURBSext, fec, orig.vdim, orig.ordering);
93 }
94
96 {
98  return fec->FiniteElementForGeometry(GeomType)->GetOrder();
99 }
100
102 {
103  Geometry::Type GeomType = mesh->GetFaceBaseGeometry(i);
104  return fec->FiniteElementForGeometry(GeomType)->GetOrder();
105 }
106
107 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs, int ndofs) const
108 {
109  if (vdim == 1) { return; }
110  if (ndofs < 0) { ndofs = this->ndofs; }
111
113  {
114  Ordering::DofsToVDofs<Ordering::byNODES>(ndofs, vdim, dofs);
115  }
116  else
117  {
118  Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs, vdim, dofs);
119  }
120 }
121
122 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs, int ndofs) const
123 {
124  if (vdim == 1) { return; }
125  if (ndofs < 0) { ndofs = this->ndofs; }
126
128  {
129  for (int i = 0; i < dofs.Size(); i++)
130  {
131  dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs, vdim, dofs[i], vd);
132  }
133  }
134  else
135  {
136  for (int i = 0; i < dofs.Size(); i++)
137  {
138  dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dofs[i], vd);
139  }
140  }
141 }
142
143 int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs) const
144 {
145  if (vdim == 1) { return dof; }
146  if (ndofs < 0) { ndofs = this->ndofs; }
147
149  {
150  return Ordering::Map<Ordering::byNODES>(ndofs, vdim, dof, vd);
151  }
152  else
153  {
154  return Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dof, vd);
155  }
156 }
157
158 // static function
160 {
161  int n = vdofs.Size(), *vdof = vdofs;
162  for (int i = 0; i < n; i++)
163  {
164  int j;
165  if ((j = vdof[i]) < 0)
166  {
167  vdof[i] = -1-j;
168  }
169  }
170 }
171
173 {
174  GetElementDofs(i, vdofs);
175  DofsToVDofs(vdofs);
176 }
177
179 {
180  GetBdrElementDofs(i, vdofs);
181  DofsToVDofs(vdofs);
182 }
183
185 {
186  GetFaceDofs(i, vdofs);
187  DofsToVDofs(vdofs);
188 }
189
191 {
192  GetEdgeDofs(i, vdofs);
193  DofsToVDofs(vdofs);
194 }
195
197 {
198  GetVertexDofs(i, vdofs);
199  DofsToVDofs(vdofs);
200 }
201
203 {
204  GetElementInteriorDofs(i, vdofs);
205  DofsToVDofs(vdofs);
206 }
207
209 {
210  GetEdgeInteriorDofs(i, vdofs);
211  DofsToVDofs(vdofs);
212 }
213
215 {
216  if (elem_dof) { return; }
217
218  Table *el_dof = new Table;
219  Array<int> dofs;
220  el_dof -> MakeI (mesh -> GetNE());
221  for (int i = 0; i < mesh -> GetNE(); i++)
222  {
223  GetElementDofs (i, dofs);
224  el_dof -> AddColumnsInRow (i, dofs.Size());
225  }
226  el_dof -> MakeJ();
227  for (int i = 0; i < mesh -> GetNE(); i++)
228  {
229  GetElementDofs (i, dofs);
230  el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
231  }
232  el_dof -> ShiftUpI();
233  elem_dof = el_dof;
234 }
235
237 {
238  delete elem_dof;
239  elem_dof = NULL;
241 }
242
244 {
245  Array<int> dof_marker(ndofs);
246
247  dof_marker = -1;
248
249  int *J = elem_dof->GetJ(), nnz = elem_dof->Size_of_connections();
250  for (int k = 0, dof_counter = 0; k < nnz; k++)
251  {
252  const int sdof = J[k]; // signed dof
253  const int dof = (sdof < 0) ? -1-sdof : sdof;
254  int new_dof = dof_marker[dof];
255  if (new_dof < 0)
256  {
257  dof_marker[dof] = new_dof = dof_counter++;
258  }
259  J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
260  }
261 }
262
264 {
265  if (dof_elem_array.Size()) { return; }
266
268
271  dof_elem_array = -1;
272  for (int i = 0; i < mesh -> GetNE(); i++)
273  {
274  const int *dofs = elem_dof -> GetRow(i);
275  const int n = elem_dof -> RowSize(i);
276  for (int j = 0; j < n; j++)
277  {
278  if (dof_elem_array[dofs[j]] < 0)
279  {
280  dof_elem_array[dofs[j]] = i;
281  dof_ldof_array[dofs[j]] = j;
282  }
283  }
284  }
285 }
286
287 static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
288 {
289  for (int i = 0; i < dofs.Size(); i++)
290  {
291  int k = dofs[i];
292  if (k < 0) { k = -1 - k; }
293  mark_array[k] = -1;
294  }
295 }
296
298  Array<int> &ess_vdofs,
299  int component) const
300 {
301  Array<int> vdofs, dofs;
302
303  ess_vdofs.SetSize(GetVSize());
304  ess_vdofs = 0;
305
306  for (int i = 0; i < GetNBE(); i++)
307  {
308  if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
309  {
310  if (component < 0)
311  {
312  // Mark all components.
313  GetBdrElementVDofs(i, vdofs);
314  mark_dofs(vdofs, ess_vdofs);
315  }
316  else
317  {
318  GetBdrElementDofs(i, dofs);
319  for (int d = 0; d < dofs.Size(); d++)
320  { dofs[d] = DofToVDof(dofs[d], component); }
321  mark_dofs(dofs, ess_vdofs);
322  }
323  }
324  }
325
326  // mark possible hidden boundary edges in a non-conforming mesh, also
327  // local DOFs affected by boundary elements on other processors
328  if (Nonconforming())
329  {
330  Array<int> bdr_verts, bdr_edges;
331  mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges);
332
333  for (int i = 0; i < bdr_verts.Size(); i++)
334  {
335  if (component < 0)
336  {
337  GetVertexVDofs(bdr_verts[i], vdofs);
338  mark_dofs(vdofs, ess_vdofs);
339  }
340  else
341  {
342  GetVertexDofs(bdr_verts[i], dofs);
343  for (int d = 0; d < dofs.Size(); d++)
344  { dofs[d] = DofToVDof(dofs[d], component); }
345  mark_dofs(dofs, ess_vdofs);
346  }
347  }
348  for (int i = 0; i < bdr_edges.Size(); i++)
349  {
350  if (component < 0)
351  {
352  GetEdgeVDofs(bdr_edges[i], vdofs);
353  mark_dofs(vdofs, ess_vdofs);
354  }
355  else
356  {
357  GetEdgeDofs(bdr_edges[i], dofs);
358  for (int d = 0; d < dofs.Size(); d++)
359  { dofs[d] = DofToVDof(dofs[d], component); }
360  mark_dofs(dofs, ess_vdofs);
361  }
362  }
363  }
364 }
365
367  Array<int> &ess_tdof_list,
368  int component)
369 {
370  Array<int> ess_vdofs, ess_tdofs;
371  GetEssentialVDofs(bdr_attr_is_ess, ess_vdofs, component);
373  if (!R)
374  {
375  ess_tdofs.MakeRef(ess_vdofs);
376  }
377  else
378  {
379  R->BooleanMult(ess_vdofs, ess_tdofs);
380  }
381  MarkerToList(ess_tdofs, ess_tdof_list);
382 }
383
384 // static method
386  Array<int> &list)
387 {
388  int num_marked = 0;
389  marker.HostRead(); // make sure we can read the array on host
390  for (int i = 0; i < marker.Size(); i++)
391  {
392  if (marker[i]) { num_marked++; }
393  }
394  list.SetSize(0);
395  list.Reserve(num_marked);
396  for (int i = 0; i < marker.Size(); i++)
397  {
398  if (marker[i]) { list.Append(i); }
399  }
400 }
401
402 // static method
403 void FiniteElementSpace::ListToMarker(const Array<int> &list, int marker_size,
404  Array<int> &marker, int mark_val)
405 {
406  marker.SetSize(marker_size);
407  marker = 0;
408  for (int i = 0; i < list.Size(); i++)
409  {
410  marker[list[i]] = mark_val;
411  }
412 }
413
415  Array<int> &cdofs)
416 {
418  if (cP) { cP->BooleanMultTranspose(dofs, cdofs); }
419  else { dofs.Copy(cdofs); }
420 }
421
423  Array<int> &dofs)
424 {
426  if (cR) { cR->BooleanMultTranspose(cdofs, dofs); }
427  else { cdofs.Copy(dofs); }
428 }
429
430 SparseMatrix *
432 {
433  int i, j;
434  Array<int> d_vdofs, c_vdofs;
435  SparseMatrix *R;
436
437  R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
438
439  for (i = 0; i < mesh -> GetNE(); i++)
440  {
441  this -> GetElementVDofs (i, d_vdofs);
442  cfes -> GetElementVDofs (i, c_vdofs);
443
444 #ifdef MFEM_DEBUG
445  if (d_vdofs.Size() != c_vdofs.Size())
446  {
447  mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
448  }
449 #endif
450
451  for (j = 0; j < d_vdofs.Size(); j++)
452  {
453  R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
454  }
455  }
456
457  R -> Finalize();
458
459  return R;
460 }
461
462 SparseMatrix *
464 {
465  int i, j;
466  Array<int> d_dofs, c_dofs;
467  SparseMatrix *R;
468
469  R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
470
471  for (i = 0; i < mesh -> GetNE(); i++)
472  {
473  this -> GetElementDofs (i, d_dofs);
474  cfes -> GetElementDofs (i, c_dofs);
475
476 #ifdef MFEM_DEBUG
477  if (c_dofs.Size() != 1)
478  mfem_error ("FiniteElementSpace::"
479  "D2Const_GlobalRestrictionMatrix (...)");
480 #endif
481
482  for (j = 0; j < d_dofs.Size(); j++)
483  {
484  R -> Set (c_dofs[0], d_dofs[j], 1.0);
485  }
486  }
487
488  R -> Finalize();
489
490  return R;
491 }
492
493 SparseMatrix *
495 {
496  SparseMatrix *R;
497  DenseMatrix loc_restr;
498  Array<int> l_dofs, h_dofs;
499
500  R = new SparseMatrix (lfes -> GetNDofs(), ndofs);
501
502  Geometry::Type cached_geom = Geometry::INVALID;
503  const FiniteElement *h_fe = NULL;
504  const FiniteElement *l_fe = NULL;
506
507  for (int i = 0; i < mesh -> GetNE(); i++)
508  {
509  this -> GetElementDofs (i, h_dofs);
510  lfes -> GetElementDofs (i, l_dofs);
511
512  // Assuming 'loc_restr' depends only on the Geometry::Type.
514  if (geom != cached_geom)
515  {
516  h_fe = this -> GetFE (i);
517  l_fe = lfes -> GetFE (i);
519  h_fe->Project(*l_fe, T, loc_restr);
520  cached_geom = geom;
521  }
522
523  R -> SetSubMatrix (l_dofs, h_dofs, loc_restr, 1);
524  }
525
526  R -> Finalize();
527
528  return R;
529 }
530
531 void
533  Array<int>& slave_dofs, DenseMatrix& I)
534 {
535  for (int i = 0; i < slave_dofs.Size(); i++)
536  {
537  int sdof = slave_dofs[i];
538  if (!deps.RowSize(sdof)) // not processed yet?
539  {
540  for (int j = 0; j < master_dofs.Size(); j++)
541  {
542  double coef = I(i, j);
543  if (std::abs(coef) > 1e-12)
544  {
545  int mdof = master_dofs[j];
546  if (mdof != sdof && mdof != (-1-sdof))
547  {
549  }
550  }
551  }
552  }
553  }
554 }
555
556 bool FiniteElementSpace::DofFinalizable(int dof, const Array<bool>& finalized,
557  const SparseMatrix& deps)
558 {
559  const int* dep = deps.GetRowColumns(dof);
560  int ndep = deps.RowSize(dof);
561
562  // are all constraining DOFs finalized?
563  for (int i = 0; i < ndep; i++)
564  {
565  if (!finalized[dep[i]]) { return false; }
566  }
567  return true;
568 }
569
571  Geometry::Type master_geom) const
572 {
573  // In NC meshes with prisms/tets, a special constraint occurs where a
574  // prism/tet edge is slave to another element's face. Rather than introduce a
575  // new edge-face constraint type, we handle such cases as degenerate
576  // face-face constraints, where the point-matrix rectangle has zero height.
577  // This method returns DOFs for the first edge of the rectangle, duplicated
578  // in the orthogonal direction, to resemble DOFs for a quadrilateral face.
579  // The extra DOFs are ignored by FiniteElementSpace::AddDependencies.
580
581  Array<int> edof;
582  GetEdgeDofs(-1 - index, edof);
583
586  int nn = 2*nv + ne;
587
588  dofs.SetSize(nn*nn);
589  if (!dofs.Size()) { return; }
590
591  dofs = edof[0];
592
593  // copy first two vertex DOFs
594  for (int i = 0; i < nv; i++)
595  {
596  dofs[i] = edof[i];
597  dofs[nv+i] = edof[nv+i];
598  }
599  // copy first edge DOFs
600  int face_vert = Geometry::NumVerts[master_geom];
601  for (int i = 0; i < ne; i++)
602  {
603  dofs[face_vert*nv + i] = edof[2*nv + i];
604  }
605 }
606
607 void
609  Geometry::Type master_geom) const
610 {
611  switch (entity)
612  {
613  case 0: GetVertexDofs(index, dofs); break;
614  case 1: GetEdgeDofs(index, dofs); break;
615  case 2: (index >= 0) ? GetFaceDofs(index, dofs)
616  /* */ : GetDegenerateFaceDofs(index, dofs, master_geom);
617  }
618 }
619
620
622 {
623 #ifdef MFEM_USE_MPI
624  MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
625  "This method should not be used with a ParFiniteElementSpace!");
626 #endif
627
628  if (cP_is_set) { return; }
629  cP_is_set = true;
630
631  // For each slave DOF, the dependency matrix will contain a row that
632  // expresses the slave DOF as a linear combination of its immediate master
633  // DOFs. Rows of independent DOFs will remain empty.
634  SparseMatrix deps(ndofs);
635
636  // collect local edge/face dependencies
637  for (int entity = 1; entity <= 2; entity++)
638  {
639  const NCMesh::NCList &list = mesh->ncmesh->GetNCList(entity);
640  if (!list.masters.size()) { continue; }
641
642  Array<int> master_dofs, slave_dofs;
643
645  DenseMatrix I;
646
647  // loop through all master edges/faces, constrain their slave edges/faces
648  for (unsigned mi = 0; mi < list.masters.size(); mi++)
649  {
650  const NCMesh::Master &master = list.masters[mi];
651
652  GetEntityDofs(entity, master.index, master_dofs);
653  if (!master_dofs.Size()) { continue; }
654
655  const FiniteElement* fe = fec->FiniteElementForGeometry(master.Geom());
656  if (!fe) { continue; }
657
658  switch (master.geom)
659  {
661  case Geometry::TRIANGLE: T.SetFE(&TriangleFE); break;
662  case Geometry::SEGMENT: T.SetFE(&SegmentFE); break;
663  default: MFEM_ABORT("unsupported geometry");
664  }
665
666  for (int si = master.slaves_begin; si < master.slaves_end; si++)
667  {
668  const NCMesh::Slave &slave = list.slaves[si];
669  GetEntityDofs(entity, slave.index, slave_dofs, master.Geom());
670  if (!slave_dofs.Size()) { continue; }
671
672  slave.OrientedPointMatrix(T.GetPointMat());
674  fe->GetLocalInterpolation(T, I);
675
676  // make each slave DOF dependent on all master DOFs
678  }
679  }
680  }
681
682  deps.Finalize();
683
684  // DOFs that stayed independent are true DOFs
685  int n_true_dofs = 0;
686  for (int i = 0; i < ndofs; i++)
687  {
688  if (!deps.RowSize(i)) { n_true_dofs++; }
689  }
690
691  // if all dofs are true dofs leave cP and cR NULL
692  if (n_true_dofs == ndofs)
693  {
694  cP = cR = NULL; // will be treated as identities
695  return;
696  }
697
698  // create the conforming restriction matrix cR
699  int *cR_J;
700  {
701  int *cR_I = Memory<int>(n_true_dofs+1);
702  double *cR_A = Memory<double>(n_true_dofs);
703  cR_J = Memory<int>(n_true_dofs);
704  for (int i = 0; i < n_true_dofs; i++)
705  {
706  cR_I[i] = i;
707  cR_A[i] = 1.0;
708  }
709  cR_I[n_true_dofs] = n_true_dofs;
710  cR = new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, ndofs);
711  }
712
713  // create the conforming prolongation matrix cP
714  cP = new SparseMatrix(ndofs, n_true_dofs);
715
716  Array<bool> finalized(ndofs);
717  finalized = false;
718
719  // put identity in the restriction and prolongation matrices for true DOFs
720  for (int i = 0, true_dof = 0; i < ndofs; i++)
721  {
722  if (!deps.RowSize(i))
723  {
724  cR_J[true_dof] = i;
726  finalized[i] = true;
727  }
728  }
729
730  // Now calculate cP rows of slave DOFs as combinations of cP rows of their
731  // master DOFs. It is possible that some slave DOFs depend on DOFs that are
732  // themselves slaves. Here we resolve such indirect constraints by first
733  // calculating rows of the cP matrix for DOFs whose master DOF cP rows are
734  // already known (in the first iteration these are the true DOFs). In the
735  // second iteration, slaves of slaves can be 'finalized' (given a row in the
736  // cP matrix), in the third iteration slaves of slaves of slaves, etc.
737  bool finished;
738  int n_finalized = n_true_dofs;
739  Array<int> cols;
740  Vector srow;
741  do
742  {
743  finished = true;
744  for (int dof = 0; dof < ndofs; dof++)
745  {
746  if (!finalized[dof] && DofFinalizable(dof, finalized, deps))
747  {
748  const int* dep_col = deps.GetRowColumns(dof);
749  const double* dep_coef = deps.GetRowEntries(dof);
750  int n_dep = deps.RowSize(dof);
751
752  for (int j = 0; j < n_dep; j++)
753  {
754  cP->GetRow(dep_col[j], cols, srow);
755  srow *= dep_coef[j];
757  }
758
759  finalized[dof] = true;
760  n_finalized++;
761  finished = false;
762  }
763  }
764  }
765  while (!finished);
766
767  // if everything is consistent (mesh, face orientations, etc.), we should
768  // be able to finalize all slave DOFs, otherwise it's a serious error
769  if (n_finalized != ndofs)
770  {
771  MFEM_ABORT("Error creating cP matrix.");
772  }
773
774  cP->Finalize();
775
776  if (vdim > 1)
777  {
778  MakeVDimMatrix(*cP);
779  MakeVDimMatrix(*cR);
780  }
781
782  if (Device::IsEnabled()) { cP->BuildTranspose(); }
783 }
784
786 {
787  if (vdim == 1) { return; }
788
789  int height = mat.Height();
790  int width = mat.Width();
791
792  SparseMatrix *vmat = new SparseMatrix(vdim*height, vdim*width);
793
794  Array<int> dofs, vdofs;
795  Vector srow;
796  for (int i = 0; i < height; i++)
797  {
798  mat.GetRow(i, dofs, srow);
799  for (int vd = 0; vd < vdim; vd++)
800  {
801  dofs.Copy(vdofs);
802  DofsToVDofs(vd, vdofs, width);
803  vmat->SetRow(DofToVDof(i, vd, height), vdofs, srow);
804  }
805  }
806  vmat->Finalize();
807
808  mat.Swap(*vmat);
809  delete vmat;
810 }
811
812
814 {
815  if (Conforming()) { return NULL; }
817  return cP;
818 }
819
821 {
822  if (Conforming()) { return NULL; }
824  return cR;
825 }
826
828 {
830  return P ? (P->Width() / vdim) : ndofs;
831 }
832
834  ElementDofOrdering e_ordering) const
835 {
836  // Check if we have a discontinuous space using the FE collection:
837  const L2_FECollection *dg_space = dynamic_cast<const L2_FECollection*>(fec);
838  if (dg_space)
839  {
840  if (L2E_nat.Ptr() == NULL)
841  {
842  L2E_nat.Reset(new L2ElementRestriction(*this));
843  }
844  return L2E_nat.Ptr();
845  }
846  if (e_ordering == ElementDofOrdering::LEXICOGRAPHIC)
847  {
848  if (L2E_lex.Ptr() == NULL)
849  {
850  L2E_lex.Reset(new ElementRestriction(*this, e_ordering));
851  }
852  return L2E_lex.Ptr();
853  }
854  // e_ordering == ElementDofOrdering::NATIVE
855  if (L2E_nat.Ptr() == NULL)
856  {
857  L2E_nat.Reset(new ElementRestriction(*this, e_ordering));
858  }
859  return L2E_nat.Ptr();
860 }
861
863  ElementDofOrdering e_ordering, FaceType type, L2FaceValues mul) const
864 {
865  const bool is_dg_space = dynamic_cast<const L2_FECollection*>(fec)!=nullptr;
866  const L2FaceValues m = (is_dg_space && mul==L2FaceValues::DoubleValued) ?
868  key_face key = std::make_tuple(is_dg_space, e_ordering, type, m);
869  auto itr = L2F.find(key);
870  if (itr != L2F.end())
871  {
872  return itr->second;
873  }
874  else
875  {
876  Operator* res;
877  if (is_dg_space)
878  {
879  res = new L2FaceRestriction(*this, e_ordering, type, m);
880  }
881  else
882  {
883  res = new H1FaceRestriction(*this, e_ordering, type);
884  }
885  L2F[key] = res;
886  return res;
887  }
888 }
889
891  const IntegrationRule &ir) const
892 {
893  for (int i = 0; i < E2Q_array.Size(); i++)
894  {
895  const QuadratureInterpolator *qi = E2Q_array[i];
896  if (qi->IntRule == &ir) { return qi; }
897  }
898
900  E2Q_array.Append(qi);
901  return qi;
902 }
903
906 {
907  for (int i = 0; i < E2Q_array.Size(); i++)
908  {
909  const QuadratureInterpolator *qi = E2Q_array[i];
910  if (qi->qspace == &qs) { return qi; }
911  }
912
914  E2Q_array.Append(qi);
915  return qi;
916 }
917
920  const IntegrationRule &ir, FaceType type) const
921 {
922  if (type==FaceType::Interior)
923  {
924  for (int i = 0; i < E2IFQ_array.Size(); i++)
925  {
927  if (qi->IntRule == &ir) { return qi; }
928  }
929
931  type);
932  E2IFQ_array.Append(qi);
933  return qi;
934  }
935  else //Boundary
936  {
937  for (int i = 0; i < E2BFQ_array.Size(); i++)
938  {
940  if (qi->IntRule == &ir) { return qi; }
941  }
942
944  type);
945  E2BFQ_array.Append(qi);
946  return qi;
947  }
948 }
949
951  const int coarse_ndofs, const Table &coarse_elem_dof,
952  const DenseTensor localP[]) const
953 {
954  MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
955
956  Array<int> dofs, coarse_dofs, coarse_vdofs;
957  Vector row;
958
959  Mesh::GeometryList elem_geoms(*mesh);
960
961  SparseMatrix *P;
962  if (elem_geoms.Size() == 1)
963  {
964  const int coarse_ldof = localP[elem_geoms[0]].SizeJ();
965  P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim, coarse_ldof);
966  }
967  else
968  {
969  P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim);
970  }
971
972  Array<int> mark(P->Height());
973  mark = 0;
974
976
977  for (int k = 0; k < mesh->GetNE(); k++)
978  {
979  const Embedding &emb = rtrans.embeddings[k];
981  const DenseMatrix &lP = localP[geom](emb.matrix);
982  const int fine_ldof = localP[geom].SizeI();
983
984  elem_dof->GetRow(k, dofs);
985  coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
986
987  for (int vd = 0; vd < vdim; vd++)
988  {
989  coarse_dofs.Copy(coarse_vdofs);
990  DofsToVDofs(vd, coarse_vdofs, coarse_ndofs);
991
992  for (int i = 0; i < fine_ldof; i++)
993  {
994  int r = DofToVDof(dofs[i], vd);
995  int m = (r >= 0) ? r : (-1 - r);
996
997  if (!mark[m])
998  {
999  lP.GetRow(i, row);
1000  P->SetRow(r, coarse_vdofs, row);
1001  mark[m] = 1;
1002  }
1003  }
1004  }
1005  }
1006
1007  MFEM_ASSERT(mark.Sum() == P->Height(), "Not all rows of P set.");
1008  if (elem_geoms.Size() != 1) { P->Finalize(); }
1009  return P;
1010 }
1011
1013  Geometry::Type geom, DenseTensor &localP) const
1014 {
1015  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
1016
1018  const DenseTensor &pmats = rtrans.point_matrices[geom];
1019
1020  int nmat = pmats.SizeK();
1021  int ldof = fe->GetDof();
1022
1024  isotr.SetIdentityTransformation(geom);
1025
1026  // calculate local interpolation matrices for all refinement types
1027  localP.SetSize(ldof, ldof, nmat);
1028  for (int i = 0; i < nmat; i++)
1029  {
1030  isotr.GetPointMat() = pmats(i);
1031  isotr.FinalizeTransformation();
1032  fe->GetLocalInterpolation(isotr, localP(i));
1033  }
1034 }
1035
1037  const Table* old_elem_dof)
1038 {
1039  MFEM_VERIFY(GetNE() >= old_elem_dof->Size(),
1040  "Previous mesh is not coarser.");
1041
1042  Mesh::GeometryList elem_geoms(*mesh);
1043
1045  for (int i = 0; i < elem_geoms.Size(); i++)
1046  {
1047  GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1048  }
1049
1050  return RefinementMatrix_main(old_ndofs, *old_elem_dof, localP);
1051 }
1052
1054 (const FiniteElementSpace* fespace, Table* old_elem_dof, int old_ndofs)
1055  : fespace(fespace)
1056  , old_elem_dof(old_elem_dof)
1057 {
1058  MFEM_VERIFY(fespace->GetNE() >= old_elem_dof->Size(),
1059  "Previous mesh is not coarser.");
1060
1061  width = old_ndofs * fespace->GetVDim();
1062  height = fespace->GetVSize();
1063
1064  Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1065
1066  for (int i = 0; i < elem_geoms.Size(); i++)
1067  {
1068  fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1069  }
1070 }
1071
1073  const FiniteElementSpace *fespace, const FiniteElementSpace *coarse_fes)
1074  : Operator(fespace->GetVSize(), coarse_fes->GetVSize()),
1075  fespace(fespace), old_elem_dof(NULL)
1076 {
1077  Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1078
1079  for (int i = 0; i < elem_geoms.Size(); i++)
1080  {
1081  fespace->GetLocalRefinementMatrices(*coarse_fes, elem_geoms[i],
1082  localP[elem_geoms[i]]);
1083  }
1084
1085  // Make a copy of the coarse elem_dof Table.
1086  old_elem_dof = new Table(coarse_fes->GetElementToDofTable());
1087 }
1088
1090 {
1091  delete old_elem_dof;
1092 }
1093
1095 ::Mult(const Vector &x, Vector &y) const
1096 {
1097  Mesh* mesh = fespace->GetMesh();
1098  const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
1099
1100  Array<int> dofs, old_dofs, old_vdofs;
1101
1102  Array<char> processed(fespace->GetVSize());
1103  processed = 0;
1104
1105  int vdim = fespace->GetVDim();
1106  int old_ndofs = width / vdim;
1107
1108  for (int k = 0; k < mesh->GetNE(); k++)
1109  {
1110  const Embedding &emb = rtrans.embeddings[k];
1111  const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
1112  const DenseMatrix &lP = localP[geom](emb.matrix);
1113
1114  fespace->GetElementDofs(k, dofs);
1115  old_elem_dof->GetRow(emb.parent, old_dofs);
1116
1117  for (int vd = 0; vd < vdim; vd++)
1118  {
1119  old_dofs.Copy(old_vdofs);
1120  fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
1121
1122  for (int i = 0; i < dofs.Size(); i++)
1123  {
1124  double rsign, osign;
1125  int r = fespace->DofToVDof(dofs[i], vd);
1126  r = DecodeDof(r, rsign);
1127
1128  if (!processed[r])
1129  {
1130  double value = 0.0;
1131  for (int j = 0; j < old_vdofs.Size(); j++)
1132  {
1133  int o = DecodeDof(old_vdofs[j], osign);
1134  value += x[o] * lP(i, j) * osign;
1135  }
1136  y[r] = value * rsign;
1137  processed[r] = 1;
1138  }
1139  }
1140  }
1141  }
1142 }
1143
1145  const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes,
1146  BilinearFormIntegrator *mass_integ)
1147  : Operator(c_fes->GetVSize(), f_fes->GetVSize()),
1148  fine_fes(f_fes)
1149 {
1150  MFEM_VERIFY(c_fes->GetOrdering() == f_fes->GetOrdering() &&
1151  c_fes->GetVDim() == f_fes->GetVDim(),
1152  "incompatible coarse and fine FE spaces");
1153
1155  Mesh *f_mesh = f_fes->GetMesh();
1156  const CoarseFineTransformations &rtrans = f_mesh->GetRefinementTransforms();
1157
1158  Mesh::GeometryList elem_geoms(*f_mesh);
1160  for (int gi = 0; gi < elem_geoms.Size(); gi++)
1161  {
1162  const Geometry::Type geom = elem_geoms[gi];
1163  DenseTensor &lP = localP[geom], &lM = localM[geom];
1164  const FiniteElement *fine_fe =
1165  f_fes->fec->FiniteElementForGeometry(geom);
1166  const FiniteElement *coarse_fe =
1167  c_fes->fec->FiniteElementForGeometry(geom);
1168  const DenseTensor &pmats = rtrans.point_matrices[geom];
1169
1170  lP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), pmats.SizeK());
1171  lM.SetSize(fine_fe->GetDof(), fine_fe->GetDof(), pmats.SizeK());
1172  emb_tr.SetIdentityTransformation(geom);
1173  for (int i = 0; i < pmats.SizeK(); i++)
1174  {
1175  emb_tr.GetPointMat() = pmats(i);
1176  emb_tr.FinalizeTransformation();
1177  // Get the local interpolation matrix for this refinement type
1178  fine_fe->GetTransferMatrix(*coarse_fe, emb_tr, lP(i));
1179  // Get the local mass matrix for this refinement type
1180  mass_integ->AssembleElementMatrix(*fine_fe, emb_tr, lM(i));
1181  }
1182  }
1183
1184  Table ref_type_to_matrix;
1185  rtrans.GetCoarseToFineMap(*f_mesh, coarse_to_fine, coarse_to_ref_type,
1186  ref_type_to_matrix, ref_type_to_geom);
1187  MFEM_ASSERT(coarse_to_fine.Size() == c_fes->GetNE(), "");
1188
1189  const int total_ref_types = ref_type_to_geom.Size();
1190  int num_ref_types[Geometry::NumGeom], num_fine_elems[Geometry::NumGeom];
1191  Array<int> ref_type_to_coarse_elem_offset(total_ref_types);
1192  ref_type_to_fine_elem_offset.SetSize(total_ref_types);
1193  std::fill(num_ref_types, num_ref_types+Geometry::NumGeom, 0);
1194  std::fill(num_fine_elems, num_fine_elems+Geometry::NumGeom, 0);
1195  for (int i = 0; i < total_ref_types; i++)
1196  {
1197  Geometry::Type g = ref_type_to_geom[i];
1198  ref_type_to_coarse_elem_offset[i] = num_ref_types[g];
1199  ref_type_to_fine_elem_offset[i] = num_fine_elems[g];
1200  num_ref_types[g]++;
1201  num_fine_elems[g] += ref_type_to_matrix.RowSize(i);
1202  }
1203  DenseTensor localPtMP[Geometry::NumGeom];
1204  for (int g = 0; g < Geometry::NumGeom; g++)
1205  {
1206  if (num_ref_types[g] == 0) { continue; }
1207  const int fine_dofs = localP[g].SizeI();
1208  const int coarse_dofs = localP[g].SizeJ();
1209  localPtMP[g].SetSize(coarse_dofs, coarse_dofs, num_ref_types[g]);
1210  localR[g].SetSize(coarse_dofs, fine_dofs, num_fine_elems[g]);
1211  }
1212  for (int i = 0; i < total_ref_types; i++)
1213  {
1214  Geometry::Type g = ref_type_to_geom[i];
1215  DenseMatrix &lPtMP = localPtMP[g](ref_type_to_coarse_elem_offset[i]);
1216  int lR_offset = ref_type_to_fine_elem_offset[i]; // offset in localR[g]
1217  const int *mi = ref_type_to_matrix.GetRow(i);
1218  const int nm = ref_type_to_matrix.RowSize(i);
1219  lPtMP = 0.0;
1220  for (int s = 0; s < nm; s++)
1221  {
1222  DenseMatrix &lP = localP[g](mi[s]);
1223  DenseMatrix &lM = localM[g](mi[s]);
1224  DenseMatrix &lR = localR[g](lR_offset+s);
1225  MultAtB(lP, lM, lR); // lR = lP^T lM
1226  AddMult(lR, lP, lPtMP); // lPtMP += lP^T lM lP
1227  }
1228  DenseMatrixInverse lPtMP_inv(lPtMP);
1229  for (int s = 0; s < nm; s++)
1230  {
1231  DenseMatrix &lR = localR[g](lR_offset+s);
1232  lPtMP_inv.Mult(lR); // lR <- (P^T M P)^{-1} P^T M
1233  }
1234  }
1235
1236  // Make a copy of the coarse element-to-dof Table.
1237  coarse_elem_dof = new Table(c_fes->GetElementToDofTable());
1238 }
1239
1241 {
1242  delete coarse_elem_dof;
1243 }
1244
1246 ::Mult(const Vector &x, Vector &y) const
1247 {
1248  Array<int> c_vdofs, f_vdofs;
1249  Vector loc_x, loc_y;
1250  DenseMatrix loc_x_mat, loc_y_mat;
1251  const int vdim = fine_fes->GetVDim();
1252  const int coarse_ndofs = height/vdim;
1253  for (int coarse_el = 0; coarse_el < coarse_to_fine.Size(); coarse_el++)
1254  {
1255  coarse_elem_dof->GetRow(coarse_el, c_vdofs);
1256  fine_fes->DofsToVDofs(c_vdofs, coarse_ndofs);
1257  loc_y.SetSize(c_vdofs.Size());
1258  loc_y = 0.0;
1259  loc_y_mat.UseExternalData(loc_y.GetData(), c_vdofs.Size()/vdim, vdim);
1260  const int ref_type = coarse_to_ref_type[coarse_el];
1261  const Geometry::Type geom = ref_type_to_geom[ref_type];
1262  const int *fine_elems = coarse_to_fine.GetRow(coarse_el);
1263  const int num_fine_elems = coarse_to_fine.RowSize(coarse_el);
1264  const int lR_offset = ref_type_to_fine_elem_offset[ref_type];
1265  for (int s = 0; s < num_fine_elems; s++)
1266  {
1267  const DenseMatrix &lR = localR[geom](lR_offset+s);
1268  fine_fes->GetElementVDofs(fine_elems[s], f_vdofs);
1269  x.GetSubVector(f_vdofs, loc_x);
1270  loc_x_mat.UseExternalData(loc_x.GetData(), f_vdofs.Size()/vdim, vdim);
1272  }
1273  y.SetSubVector(c_vdofs, loc_y);
1274  }
1275 }
1276
1278  DenseTensor &localR) const
1279 {
1280  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
1281
1282  const CoarseFineTransformations &dtrans =
1284  const DenseTensor &pmats = dtrans.point_matrices[geom];
1285
1286  const int nmat = pmats.SizeK();
1287  const int ldof = fe->GetDof();
1288
1290  isotr.SetIdentityTransformation(geom);
1291
1292  // calculate local restriction matrices for all refinement types
1293  localR.SetSize(ldof, ldof, nmat);
1294  for (int i = 0; i < nmat; i++)
1295  {
1296  isotr.GetPointMat() = pmats(i);
1297  isotr.FinalizeTransformation();
1298
1299  fe->GetLocalRestriction(isotr, localR(i));
1300  }
1301 }
1302
1304  const Table* old_elem_dof)
1305 {
1306  MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
1307  MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
1308  MFEM_VERIFY(ndofs <= old_ndofs, "Previous space is not finer.");
1309
1310  Array<int> dofs, old_dofs, old_vdofs;
1311  Vector row;
1312
1313  Mesh::GeometryList elem_geoms(*mesh);
1314
1316  for (int i = 0; i < elem_geoms.Size(); i++)
1317  {
1318  GetLocalDerefinementMatrices(elem_geoms[i], localR[elem_geoms[i]]);
1319  }
1320
1321  SparseMatrix *R = (elem_geoms.Size() != 1)
1322  ? new SparseMatrix(ndofs*vdim, old_ndofs*vdim) // variable row size
1323  : new SparseMatrix(ndofs*vdim, old_ndofs*vdim,
1324  localR[elem_geoms[0]].SizeI());
1325
1326  Array<int> mark(R->Height());
1327  mark = 0;
1328
1329  const CoarseFineTransformations &dtrans =
1331
1332  MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
1333
1334  int num_marked = 0;
1335  for (int k = 0; k < dtrans.embeddings.Size(); k++)
1336  {
1337  const Embedding &emb = dtrans.embeddings[k];
1339  DenseMatrix &lR = localR[geom](emb.matrix);
1340
1341  elem_dof->GetRow(emb.parent, dofs);
1342  old_elem_dof->GetRow(k, old_dofs);
1343
1344  for (int vd = 0; vd < vdim; vd++)
1345  {
1346  old_dofs.Copy(old_vdofs);
1347  DofsToVDofs(vd, old_vdofs, old_ndofs);
1348
1349  for (int i = 0; i < lR.Height(); i++)
1350  {
1351  if (!std::isfinite(lR(i, 0))) { continue; }
1352
1353  int r = DofToVDof(dofs[i], vd);
1354  int m = (r >= 0) ? r : (-1 - r);
1355
1356  if (!mark[m])
1357  {
1358  lR.GetRow(i, row);
1359  R->SetRow(r, old_vdofs, row);
1360  mark[m] = 1;
1361  num_marked++;
1362  }
1363  }
1364  }
1365  }
1366
1367  MFEM_VERIFY(num_marked == R->Height(),
1368  "internal error: not all rows of R were set.");
1369
1370  R->Finalize(); // no-op if fixed width
1371  return R;
1372 }
1373
1375  const FiniteElementSpace &coarse_fes, Geometry::Type geom,
1376  DenseTensor &localP) const
1377 {
1378  // Assumptions: see the declaration of the method.
1379
1380  const FiniteElement *fine_fe = fec->FiniteElementForGeometry(geom);
1381  const FiniteElement *coarse_fe =
1382  coarse_fes.fec->FiniteElementForGeometry(geom);
1383
1385  const DenseTensor &pmats = rtrans.point_matrices[geom];
1386
1387  int nmat = pmats.SizeK();
1388
1390  isotr.SetIdentityTransformation(geom);
1391
1392  // Calculate the local interpolation matrices for all refinement types
1393  localP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), nmat);
1394  for (int i = 0; i < nmat; i++)
1395  {
1396  isotr.GetPointMat() = pmats(i);
1397  isotr.FinalizeTransformation();
1398  fine_fe->GetTransferMatrix(*coarse_fe, isotr, localP(i));
1399  }
1400 }
1401
1404  int vdim, int ordering)
1405 {
1406  this->mesh = mesh;
1407  this->fec = fec;
1408  this->vdim = vdim;
1409  this->ordering = (Ordering::Type) ordering;
1410
1411  elem_dof = NULL;
1412  sequence = mesh->GetSequence();
1414
1415  const NURBSFECollection *nurbs_fec =
1416  dynamic_cast<const NURBSFECollection *>(fec);
1417  if (nurbs_fec)
1418  {
1419  if (!mesh->NURBSext)
1420  {
1421  mfem_error("FiniteElementSpace::FiniteElementSpace :\n"
1422  " NURBS FE space requires NURBS mesh.");
1423  }
1424
1425  if (NURBSext == NULL)
1426  {
1427  this->NURBSext = mesh->NURBSext;
1428  own_ext = 0;
1429  }
1430  else
1431  {
1432  this->NURBSext = NURBSext;
1433  own_ext = 1;
1434  }
1435  UpdateNURBS();
1436  cP = cR = NULL;
1437  cP_is_set = false;
1438  }
1439  else
1440  {
1441  this->NURBSext = NULL;
1442  own_ext = 0;
1443  Construct();
1444  }
1446 }
1447
1449 {
1450  if (NURBSext && !own_ext)
1451  {
1452  mfem_error("FiniteElementSpace::StealNURBSext");
1453  }
1454  own_ext = 0;
1455
1456  return NURBSext;
1457 }
1458
1460 {
1461  nvdofs = 0;
1462  nedofs = 0;
1463  nfdofs = 0;
1464  nbdofs = 0;
1465  fdofs = NULL;
1466  bdofs = NULL;
1467
1468  dynamic_cast<const NURBSFECollection *>(fec)->Reset();
1469
1470  ndofs = NURBSext->GetNDof();
1473 }
1474
1476 {
1477  // This method should be used only for non-NURBS spaces.
1478  MFEM_VERIFY(!NURBSext, "internal error");
1479
1480  elem_dof = NULL;
1481  bdrElem_dof = NULL;
1482
1483  ndofs = 0;
1484  nedofs = nfdofs = nbdofs = 0;
1485  bdofs = NULL;
1486  fdofs = NULL;
1487  cP = NULL;
1488  cR = NULL;
1489  cP_is_set = false;
1490  // 'Th' is initialized/destroyed before this method is called.
1491
1493
1494  if (mesh->Dimension() > 1)
1495  {
1497  }
1498
1499  if (mesh->GetNFaces() > 0)
1500  {
1501  bool have_face_dofs = false;
1502  for (int g = Geometry::DimStart[2]; g < Geometry::DimStart[3]; g++)
1503  {
1504  if (mesh->HasGeometry(Geometry::Type(g)) &&
1506  {
1507  have_face_dofs = true;
1508  break;
1509  }
1510  }
1511  if (have_face_dofs)
1512  {
1513  fdofs = new int[mesh->GetNFaces()+1];
1514  fdofs[0] = 0;
1515  for (int i = 0; i < mesh->GetNFaces(); i++)
1516  {
1518  fdofs[i+1] = nfdofs;
1519  }
1520  }
1521  }
1522
1523  if (mesh->Dimension() > 0)
1524  {
1525  bdofs = new int[mesh->GetNE()+1];
1526  bdofs[0] = 0;
1527  for (int i = 0; i < mesh->GetNE(); i++)
1528  {
1530  bdofs[i+1] = nbdofs;
1531  }
1532  }
1533
1534  ndofs = nvdofs + nedofs + nfdofs + nbdofs;
1535
1536  // Do not build elem_dof Table here: in parallel it has to be constructed
1537  // later.
1538 }
1539
1541 {
1542  if (elem_dof)
1543  {
1544  elem_dof -> GetRow (i, dofs);
1545  }
1546  else
1547  {
1548  Array<int> V, E, Eo, F, Fo;
1549  int k, j, nv, ne, nf, nb, nfd, nd, dim;
1550  const int *ind;
1551
1552  dim = mesh->Dimension();
1554  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
1555  nb = (dim > 0) ? fec->DofForGeometry(mesh->GetElementBaseGeometry(i)) : 0;
1556  if (nv > 0)
1557  {
1558  mesh->GetElementVertices(i, V);
1559  }
1560  if (ne > 0)
1561  {
1562  mesh->GetElementEdges(i, E, Eo);
1563  }
1564  nfd = 0;
1565  if (dim == 3)
1566  {
1568  {
1569  mesh->GetElementFaces(i, F, Fo);
1570  for (k = 0; k < F.Size(); k++)
1571  {
1572  nfd += fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
1573  }
1574  }
1575  }
1576  nd = V.Size() * nv + E.Size() * ne + nfd + nb;
1577  dofs.SetSize(nd);
1578  if (nv > 0)
1579  {
1580  for (k = 0; k < V.Size(); k++)
1581  {
1582  for (j = 0; j < nv; j++)
1583  {
1584  dofs[k*nv+j] = V[k]*nv+j;
1585  }
1586  }
1587  nv *= V.Size();
1588  }
1589  if (ne > 0)
1590  {
1591  // if (dim > 1)
1592  for (k = 0; k < E.Size(); k++)
1593  {
1595  for (j = 0; j < ne; j++)
1596  {
1597  if (ind[j] < 0)
1598  {
1599  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1600  }
1601  else
1602  {
1603  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1604  }
1605  }
1606  }
1607  }
1608  ne = nv + ne * E.Size();
1609  if (nfd > 0)
1610  // if (dim == 3)
1611  {
1612  for (k = 0; k < F.Size(); k++)
1613  {
1615  Fo[k]);
1617  for (j = 0; j < nf; j++)
1618  {
1619  if (ind[j] < 0)
1620  {
1621  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[F[k]]+(-1-ind[j]) );
1622  }
1623  else
1624  {
1625  dofs[ne+j] = nvdofs+nedofs+fdofs[F[k]]+ind[j];
1626  }
1627  }
1628  ne += nf;
1629  }
1630  }
1631  if (nb > 0)
1632  {
1633  k = nvdofs + nedofs + nfdofs + bdofs[i];
1634  for (j = 0; j < nb; j++)
1635  {
1636  dofs[ne+j] = k + j;
1637  }
1638  }
1639  }
1640 }
1641
1643 {
1644  if (i < 0 || !mesh->GetNE()) { return NULL; }
1645  MFEM_VERIFY(i < mesh->GetNE(),
1646  "Invalid element id " << i << ", maximum allowed " << mesh->GetNE()-1);
1647
1648  const FiniteElement *FE =
1650
1651  if (NURBSext)
1652  {
1654  }
1655
1656  return FE;
1657 }
1658
1660 {
1661  if (bdrElem_dof)
1662  {
1663  bdrElem_dof->GetRow(i, dofs);
1664  }
1665  else
1666  {
1667  Array<int> V, E, Eo;
1668  int k, j, nv, ne, nf, nd, iF, oF, dim;
1669  const int *ind;
1670
1671  dim = mesh->Dimension();
1673  if (nv > 0)
1674  {
1675  mesh->GetBdrElementVertices(i, V);
1676  }
1677  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
1678  if (ne > 0)
1679  {
1680  mesh->GetBdrElementEdges(i, E, Eo);
1681  }
1682  nd = V.Size() * nv + E.Size() * ne;
1683  nf = (dim == 3) ? (fec->DofForGeometry(
1684  mesh->GetBdrElementBaseGeometry(i))) : (0);
1685  if (nf > 0)
1686  {
1687  nd += nf;
1688  mesh->GetBdrElementFace(i, &iF, &oF);
1689  }
1690  dofs.SetSize(nd);
1691  if (nv > 0)
1692  {
1693  for (k = 0; k < V.Size(); k++)
1694  {
1695  for (j = 0; j < nv; j++)
1696  {
1697  dofs[k*nv+j] = V[k]*nv+j;
1698  }
1699  }
1700  nv *= V.Size();
1701  }
1702  if (ne > 0)
1703  {
1704  // if (dim > 1)
1705  for (k = 0; k < E.Size(); k++)
1706  {
1708  for (j = 0; j < ne; j++)
1709  {
1710  if (ind[j] < 0)
1711  {
1712  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1713  }
1714  else
1715  {
1716  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1717  }
1718  }
1719  }
1720  }
1721  if (nf > 0)
1722  // if (dim == 3)
1723  {
1724  ne = nv + ne * E.Size();
1725  ind = fec->DofOrderForOrientation(
1727  for (j = 0; j < nf; j++)
1728  {
1729  if (ind[j] < 0)
1730  {
1731  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[iF]+(-1-ind[j]) );
1732  }
1733  else
1734  {
1735  dofs[ne+j] = nvdofs+nedofs+fdofs[iF]+ind[j];
1736  }
1737  }
1738  }
1739  }
1740 }
1741
1743 {
1744  int j, k, nv, ne, nf, nd, dim = mesh->Dimension();
1745  Array<int> V, E, Eo;
1746  const int *ind;
1747
1748  // for 1D, 2D and 3D faces
1750  ne = (dim > 1) ? fec->DofForGeometry(Geometry::SEGMENT) : 0;
1751  if (nv > 0)
1752  {
1753  mesh->GetFaceVertices(i, V);
1754  }
1755  if (ne > 0)
1756  {
1757  mesh->GetFaceEdges(i, E, Eo);
1758  }
1759  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
1760  nd = V.Size() * nv + E.Size() * ne + nf;
1761  dofs.SetSize(nd);
1762  if (nv > 0)
1763  {
1764  for (k = 0; k < V.Size(); k++)
1765  {
1766  for (j = 0; j < nv; j++)
1767  {
1768  dofs[k*nv+j] = V[k]*nv+j;
1769  }
1770  }
1771  }
1772  nv *= V.Size();
1773  if (ne > 0)
1774  {
1775  for (k = 0; k < E.Size(); k++)
1776  {
1778  for (j = 0; j < ne; j++)
1779  {
1780  if (ind[j] < 0)
1781  {
1782  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1783  }
1784  else
1785  {
1786  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1787  }
1788  }
1789  }
1790  }
1791  ne = nv + ne * E.Size();
1792  if (nf > 0)
1793  {
1794  for (j = nvdofs+nedofs+fdofs[i], k = 0; k < nf; j++, k++)
1795  {
1796  dofs[ne+k] = j;
1797  }
1798  }
1799 }
1800
1802 {
1803  int j, k, nv, ne;
1804  Array<int> V;
1805
1807  if (nv > 0)
1808  {
1809  mesh->GetEdgeVertices(i, V);
1810  }
1812  dofs.SetSize(2*nv+ne);
1813  if (nv > 0)
1814  {
1815  for (k = 0; k < 2; k++)
1816  {
1817  for (j = 0; j < nv; j++)
1818  {
1819  dofs[k*nv+j] = V[k]*nv+j;
1820  }
1821  }
1822  }
1823  nv *= 2;
1824  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
1825  {
1826  dofs[nv+j] = k;
1827  }
1828 }
1829
1831 {
1832  int j, nv;
1833
1835  dofs.SetSize(nv);
1836  for (j = 0; j < nv; j++)
1837  {
1838  dofs[j] = i*nv+j;
1839  }
1840 }
1841
1843 {
1844  int j, k, nb;
1845  if (mesh->Dimension() == 0) { dofs.SetSize(0); return; }
1846  nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (i));
1847  dofs.SetSize (nb);
1848  k = nvdofs + nedofs + nfdofs + bdofs[i];
1849  for (j = 0; j < nb; j++)
1850  {
1851  dofs[j] = k + j;
1852  }
1853 }
1854
1856 {
1857  int j, k, ne;
1858
1859  ne = fec -> DofForGeometry (Geometry::SEGMENT);
1860  dofs.SetSize (ne);
1861  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
1862  {
1863  dofs[j] = k;
1864  }
1865 }
1866
1868 {
1869  int j, k, nf;
1870
1871  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
1872  dofs.SetSize (nf);
1873  if (nf > 0)
1874  {
1875  for (j = 0, k = nvdofs+nedofs+fdofs[i]; j < nf; j++, k++)
1876  {
1877  dofs[j] = k;
1878  }
1879  }
1880 }
1881
1883 {
1884  const FiniteElement *BE;
1885
1886  switch ( mesh->Dimension() )
1887  {
1888  case 1:
1890  break;
1891  case 2:
1893  break;
1894  case 3:
1895  default:
1898  }
1899
1900  if (NURBSext)
1901  {
1903  }
1904
1905  return BE;
1906 }
1907
1909 {
1910  const FiniteElement *fe;
1911
1912  switch (mesh->Dimension())
1913  {
1914  case 1:
1916  break;
1917  case 2:
1919  break;
1920  case 3:
1921  default:
1923  }
1924
1925  // if (NURBSext)
1927
1928  return fe;
1929 }
1930
1932 {
1934 }
1935
1937  int i, Geometry::Type geom_type) const
1938 {
1939  return fec->TraceFiniteElementForGeometry(geom_type);
1940 }
1941
1943 {
1944  Destroy();
1945 }
1946
1948 {
1949  delete cR;
1950  delete cP;
1951  Th.Clear();
1952  L2E_nat.Clear();
1953  L2E_lex.Clear();
1954  for (int i = 0; i < E2Q_array.Size(); i++)
1955  {
1956  delete E2Q_array[i];
1957  }
1958  E2Q_array.SetSize(0);
1959  for (auto &x : L2F)
1960  {
1961  delete x.second;
1962  }
1963  for (int i = 0; i < E2IFQ_array.Size(); i++)
1964  {
1965  delete E2IFQ_array[i];
1966  }
1967  E2IFQ_array.SetSize(0);
1968  for (int i = 0; i < E2BFQ_array.Size(); i++)
1969  {
1970  delete E2BFQ_array[i];
1971  }
1972  E2BFQ_array.SetSize(0);
1973
1976
1977  if (NURBSext)
1978  {
1979  if (own_ext) { delete NURBSext; }
1980  }
1981  else
1982  {
1983  delete elem_dof;
1984  delete bdrElem_dof;
1985
1986  delete [] bdofs;
1987  delete [] fdofs;
1988  }
1989 }
1990
1992  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
1993 {
1994  // Assumptions: see the declaration of the method.
1995
1996  if (T.Type() == Operator::MFEM_SPARSEMAT)
1997  {
1998  Mesh::GeometryList elem_geoms(*mesh);
1999
2001  for (int i = 0; i < elem_geoms.Size(); i++)
2002  {
2003  GetLocalRefinementMatrices(coarse_fes, elem_geoms[i],
2004  localP[elem_geoms[i]]);
2005  }
2006  T.Reset(RefinementMatrix_main(coarse_fes.GetNDofs(),
2007  coarse_fes.GetElementToDofTable(),
2008  localP));
2009  }
2010  else
2011  {
2012  T.Reset(new RefinementOperator(this, &coarse_fes));
2013  }
2014 }
2015
2017  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
2018 {
2019  const SparseMatrix *coarse_P = coarse_fes.GetConformingProlongation();
2020
2021  Operator::Type req_type = T.Type();
2022  GetTransferOperator(coarse_fes, T);
2023
2024  if (req_type == Operator::MFEM_SPARSEMAT)
2025  {
2027  {
2028  T.Reset(mfem::Mult(*cR, *T.As<SparseMatrix>()));
2029  }
2030  if (coarse_P)
2031  {
2032  T.Reset(mfem::Mult(*T.As<SparseMatrix>(), *coarse_P));
2033  }
2034  }
2035  else
2036  {
2037  const int RP_case = bool(GetConformingRestriction()) + 2*bool(coarse_P);
2038  if (RP_case == 0) { return; }
2039  const bool owner = T.OwnsOperator();
2040  T.SetOperatorOwner(false);
2041  switch (RP_case)
2042  {
2043  case 1:
2044  T.Reset(new ProductOperator(cR, T.Ptr(), false, owner));
2045  break;
2046  case 2:
2047  T.Reset(new ProductOperator(T.Ptr(), coarse_P, owner, false));
2048  break;
2049  case 3:
2051  cR, T.Ptr(), coarse_P, false, owner, false));
2052  break;
2053  }
2054  }
2055 }
2056
2057 void FiniteElementSpace::Update(bool want_transform)
2058 {
2059  if (mesh->GetSequence() == sequence)
2060  {
2061  return; // mesh and space are in sync, no-op
2062  }
2063  if (want_transform && mesh->GetSequence() != sequence + 1)
2064  {
2065  MFEM_ABORT("Error in update sequence. Space needs to be updated after "
2066  "each mesh modification.");
2067  }
2068  sequence = mesh->GetSequence();
2069
2070  if (NURBSext)
2071  {
2072  UpdateNURBS();
2073  return;
2074  }
2075
2076  Table* old_elem_dof = NULL;
2077  int old_ndofs;
2078
2079  // save old DOF table
2080  if (want_transform)
2081  {
2082  old_elem_dof = elem_dof;
2083  elem_dof = NULL;
2084  old_ndofs = ndofs;
2085  }
2086
2087  Destroy(); // calls Th.Clear()
2088  Construct();
2090
2091  if (want_transform)
2092  {
2093  // calculate appropriate GridFunction transformation
2094  switch (mesh->GetLastOperation())
2095  {
2096  case Mesh::REFINE:
2097  {
2098  if (Th.Type() != Operator::MFEM_SPARSEMAT)
2099  {
2100  Th.Reset(new RefinementOperator(this, old_elem_dof, old_ndofs));
2101  // The RefinementOperator takes ownership of 'old_elem_dof', so
2102  // we no longer own it:
2103  old_elem_dof = NULL;
2104  }
2105  else
2106  {
2107  // calculate fully assembled matrix
2108  Th.Reset(RefinementMatrix(old_ndofs, old_elem_dof));
2109  }
2110  break;
2111  }
2112
2113  case Mesh::DEREFINE:
2114  {
2116  Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof));
2117  if (cP && cR)
2118  {
2119  Th.SetOperatorOwner(false);
2121  false, false, true));
2122  }
2123  break;
2124  }
2125
2126  default:
2127  break;
2128  }
2129
2130  delete old_elem_dof;
2131  }
2132 }
2133
2134 void FiniteElementSpace::Save(std::ostream &out) const
2135 {
2136  int fes_format = 90; // the original format, v0.9
2137  bool nurbs_unit_weights = false;
2138
2139  // Determine the format that should be used.
2140  if (!NURBSext)
2141  {
2142  // TODO: if this is a variable-order FE space, use fes_format = 100.
2143  }
2144  else
2145  {
2146  const NURBSFECollection *nurbs_fec =
2147  dynamic_cast<const NURBSFECollection *>(fec);
2148  MFEM_VERIFY(nurbs_fec, "invalid FE collection");
2149  nurbs_fec->SetOrder(NURBSext->GetOrder());
2150  const double eps = 5e-14;
2151  nurbs_unit_weights = (NURBSext->GetWeights().Min() >= 1.0-eps &&
2152  NURBSext->GetWeights().Max() <= 1.0+eps);
2154  (NURBSext != mesh->NURBSext && !nurbs_unit_weights) ||
2155  (NURBSext->GetMaster().Size() != 0 ))
2156  {
2157  fes_format = 100; // v1.0 format
2158  }
2159  }
2160
2161  out << (fes_format == 90 ?
2162  "FiniteElementSpace\n" : "MFEM FiniteElementSpace v1.0\n")
2163  << "FiniteElementCollection: " << fec->Name() << '\n'
2164  << "VDim: " << vdim << '\n'
2165  << "Ordering: " << ordering << '\n';
2166
2167  if (fes_format == 100) // v1.0
2168  {
2169  if (!NURBSext)
2170  {
2171  // TODO: this is a variable-order FE space --> write 'element_orders'.
2172  }
2173  else if (NURBSext != mesh->NURBSext)
2174  {
2176  {
2177  out << "NURBS_order\n" << NURBSext->GetOrder() << '\n';
2178  }
2179  else
2180  {
2181  out << "NURBS_orders\n";
2182  // 1 = do not write the size, just the entries:
2183  NURBSext->GetOrders().Save(out, 1);
2184  }
2185  // If periodic BCs are given, write connectivity
2186  if (NURBSext->GetMaster().Size() != 0 )
2187  {
2188  out <<"NURBS_periodic\n";
2189  NURBSext->GetMaster().Save(out);
2190  NURBSext->GetSlave().Save(out);
2191  }
2192  // If the weights are not unit, write them to the output:
2193  if (!nurbs_unit_weights)
2194  {
2195  out << "NURBS_weights\n";
2196  NURBSext->GetWeights().Print(out, 1);
2197  }
2198  }
2199  out << "End: MFEM FiniteElementSpace v1.0\n";
2200  }
2201 }
2202
2204 {
2205  string buff;
2206  int fes_format = 0, ord;
2207  FiniteElementCollection *r_fec;
2208
2209  Destroy();
2210
2211  input >> std::ws;
2212  getline(input, buff); // 'FiniteElementSpace'
2213  filter_dos(buff);
2214  if (buff == "FiniteElementSpace") { fes_format = 90; /* v0.9 */ }
2215  else if (buff == "MFEM FiniteElementSpace v1.0") { fes_format = 100; }
2216  else { MFEM_ABORT("input stream is not a FiniteElementSpace!"); }
2217  getline(input, buff, ' '); // 'FiniteElementCollection:'
2218  input >> std::ws;
2219  getline(input, buff);
2220  filter_dos(buff);
2221  r_fec = FiniteElementCollection::New(buff.c_str());
2222  getline(input, buff, ' '); // 'VDim:'
2223  input >> vdim;
2224  getline(input, buff, ' '); // 'Ordering:'
2225  input >> ord;
2226
2227  NURBSFECollection *nurbs_fec = dynamic_cast<NURBSFECollection*>(r_fec);
2228  NURBSExtension *NURBSext = NULL;
2229  if (fes_format == 90) // original format, v0.9
2230  {
2231  if (nurbs_fec)
2232  {
2233  MFEM_VERIFY(m->NURBSext, "NURBS FE collection requires a NURBS mesh!");
2234  const int order = nurbs_fec->GetOrder();
2235  if (order != m->NURBSext->GetOrder() &&
2237  {
2238  NURBSext = new NURBSExtension(m->NURBSext, order);
2239  }
2240  }
2241  }
2242  else if (fes_format == 100) // v1.0
2243  {
2244  while (1)
2245  {
2246  skip_comment_lines(input, '#');
2247  MFEM_VERIFY(input.good(), "error reading FiniteElementSpace v1.0");
2248  getline(input, buff);
2249  filter_dos(buff);
2250  if (buff == "NURBS_order" || buff == "NURBS_orders")
2251  {
2252  MFEM_VERIFY(nurbs_fec,
2253  buff << ": NURBS FE collection is required!");
2254  MFEM_VERIFY(m->NURBSext, buff << ": NURBS mesh is required!");
2255  MFEM_VERIFY(!NURBSext, buff << ": order redefinition!");
2256  if (buff == "NURBS_order")
2257  {
2258  int order;
2259  input >> order;
2260  NURBSext = new NURBSExtension(m->NURBSext, order);
2261  }
2262  else
2263  {
2264  Array<int> orders;
2266  NURBSext = new NURBSExtension(m->NURBSext, orders);
2267  }
2268  }
2269  else if (buff == "NURBS_periodic")
2270  {
2271  Array<int> master, slave;
2274  NURBSext->ConnectBoundaries(master,slave);
2275  }
2276  else if (buff == "NURBS_weights")
2277  {
2278  MFEM_VERIFY(NURBSext, "NURBS_weights: NURBS_orders have to be "
2279  "specified before NURBS_weights!");
2281  }
2282  else if (buff == "element_orders")
2283  {
2284  MFEM_VERIFY(!nurbs_fec, "section element_orders cannot be used "
2285  "with a NURBS FE collection");
2286  MFEM_ABORT("element_orders: not implemented yet!");
2287  }
2288  else if (buff == "End: MFEM FiniteElementSpace v1.0")
2289  {
2290  break;
2291  }
2292  else
2293  {
2294  MFEM_ABORT("unknown section: " << buff);
2295  }
2296  }
2297  }
2298
2299  Constructor(m, NURBSext, r_fec, vdim, ord);
2300
2301  return r_fec;
2302 }
2303
2304
2306 {
2307  // protected method
2308  int offset = 0;
2309  const int num_elem = mesh->GetNE();
2310  element_offsets = new int[num_elem + 1];
2311  for (int g = 0; g < Geometry::NumGeom; g++)
2312  {
2313  int_rule[g] = NULL;
2314  }
2315  for (int i = 0; i < num_elem; i++)
2316  {
2317  element_offsets[i] = offset;
2318  int geom = mesh->GetElementBaseGeometry(i);
2319  if (int_rule[geom] == NULL)
2320  {
2321  int_rule[geom] = &IntRules.Get(geom, order);
2322  }
2323  offset += int_rule[geom]->GetNPoints();
2324  }
2325  element_offsets[num_elem] = size = offset;
2326 }
2327
2329  : mesh(mesh_)
2330 {
2331  const char *msg = "invalid input stream";
2332  string ident;
2333
2334  in >> ident; MFEM_VERIFY(ident == "QuadratureSpace", msg);
2335  in >> ident; MFEM_VERIFY(ident == "Type:", msg);
2336  in >> ident;
2338  {
2339  in >> ident; MFEM_VERIFY(ident == "Order:", msg);
2340  in >> order;
2341  }
2342  else
2343  {
2344  MFEM_ABORT("unknown QuadratureSpace type: " << ident);
2345  return;
2346  }
2347
2348  Construct();
2349 }
2350
2352 {
2355  << "Order: " << order << '\n';
2356 }
2357
2358
2360  FiniteElementSpace &ran_fes_)
2361  : dom_fes(dom_fes_), ran_fes(ran_fes_),
2362  oper_type(Operator::ANY_TYPE),
2363  fw_t_oper(), bw_t_oper()
2364 {
2365 #ifdef MFEM_USE_MPI
2366  const bool par_dom = dynamic_cast<ParFiniteElementSpace*>(&dom_fes);
2367  const bool par_ran = dynamic_cast<ParFiniteElementSpace*>(&ran_fes);
2368  MFEM_VERIFY(par_dom == par_ran, "the domain and range FE spaces must both"
2369  " be either serial or parallel");
2370  parallel = par_dom;
2371 #endif
2372 }
2373
2375  FiniteElementSpace &fes_in, FiniteElementSpace &fes_out,
2376  const Operator &oper, OperatorHandle &t_oper)
2377 {
2378  if (t_oper.Ptr())
2379  {
2380  return *t_oper.Ptr();
2381  }
2382
2383  if (!Parallel())
2384  {
2385  const SparseMatrix *in_cP = fes_in.GetConformingProlongation();
2386  const SparseMatrix *out_cR = fes_out.GetConformingRestriction();
2388  {
2389  const SparseMatrix *mat = dynamic_cast<const SparseMatrix *>(&oper);
2390  MFEM_VERIFY(mat != NULL, "Operator is not a SparseMatrix");
2391  if (!out_cR)
2392  {
2393  t_oper.Reset(const_cast<SparseMatrix*>(mat), false);
2394  }
2395  else
2396  {
2397  t_oper.Reset(mfem::Mult(*out_cR, *mat));
2398  }
2399  if (in_cP)
2400  {
2401  t_oper.Reset(mfem::Mult(*t_oper.As<SparseMatrix>(), *in_cP));
2402  }
2403  }
2404  else if (oper_type == Operator::ANY_TYPE)
2405  {
2406  const int RP_case = bool(out_cR) + 2*bool(in_cP);
2407  switch (RP_case)
2408  {
2409  case 0:
2410  t_oper.Reset(const_cast<Operator*>(&oper), false);
2411  break;
2412  case 1:
2413  t_oper.Reset(
2414  new ProductOperator(out_cR, &oper, false, false));
2415  break;
2416  case 2:
2417  t_oper.Reset(
2418  new ProductOperator(&oper, in_cP, false, false));
2419  break;
2420  case 3:
2421  t_oper.Reset(
2423  out_cR, &oper, in_cP, false, false, false));
2424  break;
2425  }
2426  }
2427  else
2428  {
2429  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
2430  }
2431  }
2432  else // Parallel() == true
2433  {
2434 #ifdef MFEM_USE_MPI
2435  const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
2437  {
2438  const ParFiniteElementSpace *pfes_in =
2439  dynamic_cast<const ParFiniteElementSpace *>(&fes_in);
2440  const ParFiniteElementSpace *pfes_out =
2441  dynamic_cast<const ParFiniteElementSpace *>(&fes_out);
2442  const SparseMatrix *sp_mat = dynamic_cast<const SparseMatrix *>(&oper);
2443  const HypreParMatrix *hy_mat;
2444  if (sp_mat)
2445  {
2446  SparseMatrix *RA = mfem::Mult(*out_R, *sp_mat);
2447  t_oper.Reset(pfes_in->Dof_TrueDof_Matrix()->
2448  LeftDiagMult(*RA, pfes_out->GetTrueDofOffsets()));
2449  delete RA;
2450  }
2451  else if ((hy_mat = dynamic_cast<const HypreParMatrix *>(&oper)))
2452  {
2453  HypreParMatrix *RA =
2454  hy_mat->LeftDiagMult(*out_R, pfes_out->GetTrueDofOffsets());
2455  t_oper.Reset(mfem::ParMult(RA, pfes_in->Dof_TrueDof_Matrix()));
2456  delete RA;
2457  }
2458  else
2459  {
2460  MFEM_ABORT("unknown Operator type");
2461  }
2462  }
2463  else if (oper_type == Operator::ANY_TYPE)
2464  {
2465  t_oper.Reset(new TripleProductOperator(
2466  out_R, &oper, fes_in.GetProlongationMatrix(),
2467  false, false, false));
2468  }
2469  else
2470  {
2471  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
2472  }
2473 #endif
2474  }
2475
2476  return *t_oper.Ptr();
2477 }
2478
2479
2481 {
2482  if (own_mass_integ) { delete mass_integ; }
2483 }
2484
2486  BilinearFormIntegrator *mass_integ_, bool own_mass_integ_)
2487 {
2488  if (own_mass_integ) { delete mass_integ; }
2489
2490  mass_integ = mass_integ_;
2491  own_mass_integ = own_mass_integ_;
2492 }
2493
2495 {
2496  if (F.Ptr())
2497  {
2498  return *F.Ptr();
2499  }
2500
2501  // Costruct F
2503  {
2505  }
2506  else if (oper_type == Operator::MFEM_SPARSEMAT)
2507  {
2508  Mesh::GeometryList elem_geoms(*ran_fes.GetMesh());
2509
2511  for (int i = 0; i < elem_geoms.Size(); i++)
2512  {
2514  localP[elem_geoms[i]]);
2515  }
2518  }
2519  else
2520  {
2521  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
2522  }
2523
2524  return *F.Ptr();
2525 }
2526
2528 {
2529  if (B.Ptr())
2530  {
2531  return *B.Ptr();
2532  }
2533
2534  // Construct B, if not set, define a suitable mass_integ
2535  if (!mass_integ && ran_fes.GetNE() > 0)
2536  {
2537  const FiniteElement *f_fe_0 = ran_fes.GetFE(0);
2538  const int map_type = f_fe_0->GetMapType();
2539  if (map_type == FiniteElement::VALUE ||
2540  map_type == FiniteElement::INTEGRAL)
2541  {
2542  mass_integ = new MassIntegrator;
2543  }
2544  else if (map_type == FiniteElement::H_DIV ||
2545  map_type == FiniteElement::H_CURL)
2546  {
2548  }
2549  else
2550  {
2551  MFEM_ABORT("unknown type of FE space");
2552  }
2553  own_mass_integ = true;
2554  }
2556  {
2558  &ran_fes, &dom_fes, mass_integ));
2559  }
2560  else
2561  {
2562  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
2563  }
2564
2565  return *B.Ptr();
2566 }
2567
2568
2570  const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
2571  : fes_ho(fes_ho_), fes_lor(fes_lor_)
2572 {
2573  Mesh *mesh_ho = fes_ho.GetMesh();
2574  MFEM_VERIFY(mesh_ho->GetNumGeometries(mesh_ho->Dimension()) <= 1,
2575  "mixed meshes are not supported");
2576
2577  // If the local mesh is empty, skip all computations
2578  if (mesh_ho->GetNE() == 0) { return; }
2579
2580  const FiniteElement *fe_lor = fes_lor.GetFE(0);
2581  const FiniteElement *fe_ho = fes_ho.GetFE(0);
2582  ndof_lor = fe_lor->GetDof();
2583  ndof_ho = fe_ho->GetDof();
2584
2585  const int nel_lor = fes_lor.GetNE();
2586  const int nel_ho = fes_ho.GetNE();
2587
2588  nref = nel_lor/nel_ho;
2589
2590  // Construct the mapping from HO to LOR
2591  // ho2lor.GetRow(iho) will give all the LOR elements contained in iho
2592  ho2lor.SetSize(nel_ho, nref);
2593  const CoarseFineTransformations &cf_tr =
2594  fes_lor.GetMesh()->GetRefinementTransforms();
2595  for (int ilor=0; ilor<nel_lor; ++ilor)
2596  {
2597  int iho = cf_tr.embeddings[ilor].parent;
2599  }
2600  ho2lor.ShiftUpI();
2601
2602  // R will contain the restriction (L^2 projection operator) defined on
2603  // each coarse HO element (and corresponding patch of LOR elements)
2604  R.SetSize(ndof_lor*nref, ndof_ho, nel_ho);
2605  // P will contain the corresponding prolongation operator
2606  P.SetSize(ndof_ho, ndof_lor*nref, nel_ho);
2607
2608  DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
2609  DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
2610
2611  MassIntegrator mi;
2612  DenseMatrix M_lor_el(ndof_lor, ndof_lor);
2613  DenseMatrixInverse Minv_lor_el(&M_lor_el);
2614  DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
2615  DenseMatrix M_mixed_el(ndof_lor, ndof_ho);
2616
2617  Minv_lor = 0.0;
2618  M_lor = 0.0;
2619
2620  DenseMatrix RtMlor(ndof_ho, ndof_lor*nref);
2621  DenseMatrix RtMlorR(ndof_ho, ndof_ho);
2622  DenseMatrixInverse RtMlorR_inv(&RtMlorR);
2623
2625  IsoparametricTransformation &emb_tr = ip_tr.Transf;
2626
2627  Vector shape_ho(ndof_ho);
2628  Vector shape_lor(ndof_lor);
2629
2630  const Geometry::Type geom = fe_ho->GetGeomType();
2631  const DenseTensor &pmats = cf_tr.point_matrices[geom];
2632  emb_tr.SetIdentityTransformation(geom);
2633
2634  for (int iho=0; iho<nel_ho; ++iho)
2635  {
2636  for (int iref=0; iref<nref; ++iref)
2637  {
2638  // Assemble the low-order refined mass matrix and invert locally
2639  int ilor = ho2lor.GetRow(iho)[iref];
2640  ElementTransformation *el_tr = fes_lor.GetElementTransformation(ilor);
2641  mi.AssembleElementMatrix(*fe_lor, *el_tr, M_lor_el);
2642  M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
2643  Minv_lor_el.Factor();
2644  Minv_lor_el.GetInverseMatrix(M_lor_el);
2645  // Insert into the diagonal of the patch LOR mass matrix
2646  Minv_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
2647
2648  // Now assemble the block-row of the mixed mass matrix associated
2649  // with integrating HO functions against LOR functions on the LOR
2650  // sub-element.
2651
2652  // Create the transformation that embeds the fine low-order element
2653  // within the coarse high-order element in reference space
2654  emb_tr.GetPointMat() = pmats(cf_tr.embeddings[ilor].matrix);
2655  emb_tr.FinalizeTransformation();
2656
2657  int order = fe_lor->GetOrder() + fe_ho->GetOrder() + el_tr->OrderW();
2658  const IntegrationRule *ir = &IntRules.Get(geom, order);
2659  M_mixed_el = 0.0;
2660  for (int i = 0; i < ir->GetNPoints(); i++)
2661  {
2662  const IntegrationPoint &ip_lor = ir->IntPoint(i);
2663  IntegrationPoint ip_ho;
2664  ip_tr.Transform(ip_lor, ip_ho);
2665  fe_lor->CalcShape(ip_lor, shape_lor);
2666  fe_ho->CalcShape(ip_ho, shape_ho);
2667  el_tr->SetIntPoint(&ip_lor);
2668  // For now we use the geometry information from the LOR space
2669  // which means we won't be mass conservative if the mesh is curved
2670  double w = el_tr->Weight()*ip_lor.weight;
2671  shape_lor *= w;
2673  }
2674  M_mixed.CopyMN(M_mixed_el, iref*ndof_lor, 0);
2675  }
2676  mfem::Mult(Minv_lor, M_mixed, R(iho));
2677
2678  mfem::MultAtB(R(iho), M_lor, RtMlor);
2679  mfem::Mult(RtMlor, R(iho), RtMlorR);
2680  RtMlorR_inv.Factor();
2681  RtMlorR_inv.Mult(RtMlor, P(iho));
2682  }
2683 }
2684
2686  const Vector &x, Vector &y) const
2687 {
2688  int vdim = fes_ho.GetVDim();
2689  Array<int> vdofs;
2690  DenseMatrix xel_mat(ndof_ho, vdim);
2691  DenseMatrix yel_mat(ndof_lor*nref, vdim);
2692  for (int iho=0; iho<fes_ho.GetNE(); ++iho)
2693  {
2694  fes_ho.GetElementVDofs(iho, vdofs);
2695  x.GetSubVector(vdofs, xel_mat.GetData());
2696  mfem::Mult(R(iho), xel_mat, yel_mat);
2697  // Place result correctly into the low-order vector
2698  for (int iref=0; iref<nref; ++iref)
2699  {
2700  int ilor = ho2lor.GetRow(iho)[iref];
2701  for (int vd=0; vd<vdim; ++vd)
2702  {
2703  fes_lor.GetElementDofs(ilor, vdofs);
2704  fes_lor.DofsToVDofs(vd, vdofs);
2705  y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
2706  }
2707  }
2708  }
2709 }
2710
2712  const Vector &x, Vector &y) const
2713 {
2714  int vdim = fes_ho.GetVDim();
2715  Array<int> vdofs;
2716  DenseMatrix xel_mat(ndof_lor*nref, vdim);
2717  DenseMatrix yel_mat(ndof_ho, vdim);
2718  for (int iho=0; iho<fes_ho.GetNE(); ++iho)
2719  {
2720  // Extract the LOR DOFs
2721  for (int iref=0; iref<nref; ++iref)
2722  {
2723  int ilor = ho2lor.GetRow(iho)[iref];
2724  for (int vd=0; vd<vdim; ++vd)
2725  {
2726  fes_lor.GetElementDofs(ilor, vdofs);
2727  fes_lor.DofsToVDofs(vd, vdofs);
2728  x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
2729  }
2730  }
2731  // Locally prolongate
2732  mfem::Mult(P(iho), xel_mat, yel_mat);
2733  // Place the result in the HO vector
2734  fes_ho.GetElementVDofs(iho, vdofs);
2735  y.SetSubVector(vdofs, yel_mat.GetData());
2736  }
2737 }
2738
2740 {
2741  if (!F) { F = new L2Projection(dom_fes, ran_fes); }
2742  return *F;
2743 }
2744
2746 {
2747  if (!B)
2748  {
2749  if (!F) { F = new L2Projection(dom_fes, ran_fes); }
2750  B = new L2Prolongation(*F);
2751  }
2752  return *B;
2753 }
2754
2755 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Abstract class for Finite Elements.
Definition: fe.hpp:232
int GetOrder() const
Get the order of the NURBS collection: either a positive number, when using fixed order...
Definition: fe_coll.hpp:320
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:289
const Operator & MakeTrueOperator(FiniteElementSpace &fes_in, FiniteElementSpace &fes_out, const Operator &oper, OperatorHandle &t_oper)
Definition: fespace.cpp:2374
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:393
int Size() const
Logical size of the array.
Definition: array.hpp:124
void GetFaceEdges(int i, Array< int > &, Array< int > &) const
Definition: mesh.cpp:4401
For scalar fields; preserves point values.
Definition: fe.hpp:276
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:255
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
void GetVertexVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:196
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:2281
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:96
const Vector & GetWeights() const
Definition: nurbs.hpp:381
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2776
Read an Array from the stream in using format fmt. The format fmt can be:
Definition: array.cpp:53
void BuildTranspose() const
Build and store internally the transpose of this matrix which will be used in the methods AddMultTran...
Definition: sparsemat.cpp:671
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:378
int * GetJ()
Definition: table.hpp:114
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
Definition: fespace.hpp:107
L2FaceValues
Definition: restriction.hpp:26
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
DenseMatrix & GetPointMat()
Definition: eltrans.hpp:321
const Array< int > & GetMaster() const
Definition: nurbs.hpp:324
void GetElementInteriorVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:202
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:143
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:4369
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:813
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:3954
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:2057
void SetSize(int dim, int connections_per_row)
Set the size and the number of connections for the table.
Definition: table.cpp:128
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:2934
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:890
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:4431
Operator that extracts Face degrees of freedom.
Definition: restriction.hpp:75
BilinearFormIntegrator * mass_integ
Ownership depends on own_mass_integ.
Definition: fespace.hpp:810
void BuildElementToDofTable() const
Definition: fespace.cpp:214
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Definition: fe.cpp:130
DerefinementOperator(const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes, BilinearFormIntegrator *mass_integ)
Definition: fespace.cpp:1144
void Prolongate(const Vector &x, Vector &y) const
Definition: fespace.cpp:2711
static const int NumGeom
Definition: geom.hpp:41
FiniteElementSpace & dom_fes
Domain FE space.
Definition: fespace.hpp:718
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:1936
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:422
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:4323
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_Int *row_starts=NULL) const
Multiply the HypreParMatrix on the left by a block-diagonal parallel matrix D and return the result a...
Definition: hypre.cpp:1103
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof)
Definition: fespace.cpp:1036
SparseMatrix * DerefinementMatrix(int old_ndofs, const Table *old_elem_dof)
Calculate GridFunction restriction matrix after mesh derefinement.
Definition: fespace.cpp:1303
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:172
Ordering::Type ordering
Definition: fespace.hpp:104
const Geometry::Type geom
Definition: ex1.cpp:40
static int DecodeDof(int dof, double &sign)
Helper to remove encoded sign from a DOF.
Definition: fespace.hpp:162
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:302
Definition: array.hpp:277
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const DenseTensor localP[]) const
Definition: fespace.cpp:950
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:729
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:193
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:4563
int SizeK() const
Definition: densemat.hpp:764
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:91
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:812
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:472
int GetFaceOrder(int i) const
Returns the order of the i&#39;th face finite element.
Definition: fespace.cpp:101
int GetNKV() const
Definition: nurbs.hpp:350
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:318
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Definition: fespace.cpp:297
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition: mesh.hpp:852
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int vdim
Vector dimension (number of unknowns per degree of freedom).
Definition: fespace.hpp:99
Geometry::Type Geom() const
Definition: ncmesh.hpp:163
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:191
OperatorHandle L2E_lex
Definition: fespace.hpp:132
Operator::Type oper_type
Desired Operator::Type for the construction of all operators defined by the underlying transfer algor...
Definition: fespace.hpp:724
void Save(std::ostream &out, int fmt=0) const
Save the Array to the stream out using the format fmt. The format fmt can be:
Definition: array.cpp:40
int GetBdrAttribute(int i) const
Definition: fespace.hpp:455
void SetSize(int i, int j, int k)
Definition: densemat.hpp:768
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: fespace.cpp:366
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
Definition: ncmesh.cpp:2792
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.cpp:111
void GetDegenerateFaceDofs(int index, Array< int > &dofs, Geometry::Type master_geom) const
Definition: fespace.cpp:570
T Sum()
Sum all entries.
Definition: array.cpp:115
Geometry::Type GetFaceBaseGeometry(int i) const
Definition: mesh.hpp:785
Abstract parallel finite element space.
Definition: pfespace.hpp:28
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1908
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:2016
int GetMapType() const
Definition: fe.hpp:334
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3250
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
Definition: fespace.cpp:2203
virtual const Operator & ForwardOperator()
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
Definition: fespace.cpp:2739
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3238
int HasFaceDofs(Geometry::Type GeomType) const
Definition: fe_coll.cpp:25
static bool IsEnabled()
Return true if any backend other than Backend::CPU is enabled.
Definition: device.hpp:225
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:100
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2504
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:316
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:166
const FiniteElementCollection * fec
Associated FE collection (not owned).
Definition: fespace.hpp:96
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
Definition: fespace.hpp:129
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:1610
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:790
int Size_of_connections() const
Definition: table.hpp:98
void GetLocalDerefinementMatrices(Geometry::Type geom, DenseTensor &localR) const
Definition: fespace.cpp:1277
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
Definition: mesh.hpp:802
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: text.hpp:27
void DeleteAll()
Delete whole array.
Definition: array.hpp:802
int GetNDof() const
Definition: nurbs.hpp:360
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:759
Array< int > dof_elem_array
Definition: fespace.hpp:115
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:385
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1867
virtual int DofForGeometry(Geometry::Type GeomType) const =0
virtual void GetFaceDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1742
FiniteElementSpace & ran_fes
Range FE space.
Definition: fespace.hpp:719
void ConvertFromConformingVDofs(const Array< int > &cdofs, Array< int > &dofs)
Definition: fespace.cpp:422
const IntegrationRule * IntRule
Not owned.
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:833
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3278
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:820
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:311
ID for class SparseMatrix.
Definition: operator.hpp:232
const Array< int > & GetSlave() const
Definition: nurbs.hpp:326
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: fespace.cpp:1246
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:51
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
long GetSequence() const
Definition: mesh.hpp:1181
void GetVertexDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1830
virtual const Operator & BackwardOperator()
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Definition: fespace.cpp:2527
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:787
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:588
FaceType
Definition: mesh.hpp:42
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:231
const FiniteElement * GetEdgeElement(int i) const
Definition: fespace.cpp:1931
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:2417
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: fespace.cpp:1095
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2057
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:417
DenseTensor point_matrices[Geometry::NumGeom]
Matrices for IsoparametricTransformation organized by Geometry::Type.
Definition: ncmesh.hpp:63
Data type sparse matrix.
Definition: sparsemat.hpp:40
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:707
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
Operator that extracts Face degrees of freedom.
Definition: restriction.hpp:99
For scalar fields; preserves volume integrals.
Definition: fe.hpp:277
OperatorHandle F
Forward, coarse-to-fine, operator.
Definition: fespace.hpp:813
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:153
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
GridTransfer(FiniteElementSpace &dom_fes_, FiniteElementSpace &ran_fes_)
Definition: fespace.cpp:2359
std::vector< Master > masters
Definition: ncmesh.hpp:196
Type
Ordering methods:
Definition: fespace.hpp:32
SparseMatrix * cR
Conforming restriction matrix such that cR.cP=I.
Definition: fespace.hpp:125
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1271
const IntegrationRule * IntRule
Not owned.
SparseMatrix * D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Definition: fespace.cpp:463
static const int NumVerts[NumGeom]
Definition: geom.hpp:48
SparseMatrix * D2C_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Definition: fespace.cpp:431
Definition: table.hpp:80
void GetCoarseToFineMap(const Mesh &fine_mesh, Table &coarse_to_fine, Array< int > &coarse_to_ref_type, Table &ref_type_to_matrix, Array< Geometry::Type > &ref_type_to_geom) const
Definition: ncmesh.cpp:4043
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:281
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:311
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:1540
void BuildConformingInterpolation() const
Calculate the cP and cR matrices for a nonconforming mesh.
Definition: fespace.cpp:621
NURBSExtension * StealNURBSext()
Definition: fespace.cpp:1448
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:142
Array< Embedding > embeddings
Fine element positions in their parents.
Definition: ncmesh.hpp:65
virtual const FiniteElement * TraceFiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.hpp:57
L2Projection(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
Definition: fespace.cpp:2569
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I)
Definition: fespace.cpp:532
Array< int > dof_ldof_array
Definition: fespace.hpp:115
const Array< int > & GetOrders() const
Definition: nurbs.hpp:345
bool Conforming() const
Definition: fespace.hpp:301
SparseMatrix * cP
Definition: fespace.hpp:123
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
int Dimension() const
Definition: mesh.hpp:744
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:4603
Definition: fespace.hpp:148
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
void GetEdgeDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1801
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
OperatorHandle L2E_nat
The element restriction operators, see GetElementRestriction().
Definition: fespace.hpp:132
int slaves_end
slave faces
Definition: ncmesh.hpp:170
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:1991
int SizeI() const
Definition: densemat.hpp:762
OperatorHandle B
Backward, fine-to-coarse, operator.
Definition: fespace.hpp:814
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1842
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
Definition: fespace.cpp:1012
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:348
List of mesh geometries stored as Array&lt;Geometry::Type&gt;.
Definition: mesh.hpp:817
virtual const Operator & BackwardOperator()
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Definition: fespace.cpp:2745
virtual const Operator * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Definition: fespace.cpp:862
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:4601
void SetMassIntegrator(BilinearFormIntegrator *mass_integ_, bool own_mass_integ_=true)
Assign a mass integrator to be used in the construction of the backward, fine-to-coarse, transfer operator.
Definition: fespace.cpp:2485
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1072
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:441
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:229
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
RefinementOperator(const FiniteElementSpace *fespace, Table *old_elem_dof, int old_ndofs)
Definition: fespace.cpp:1054
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:2913
void SetOrder(int Order) const
Set the order and the name, based on the given Order: either a positive number for fixed order...
Definition: fe_coll.cpp:2568
void SetIdentityTransformation(Geometry::Type GeomType)
Definition: eltrans.cpp:370
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:31
std::vector< Slave > slaves
Definition: ncmesh.hpp:197
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:649
Nonconforming edge/face within a bigger edge/face.
Definition: ncmesh.hpp:178
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:314
virtual ~FiniteElementSpace()
Definition: fespace.cpp:1942
bool OwnsOperator() const
Return true if the OperatorHandle owns the held Operator.
Definition: handle.hpp:109
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Given a particular Finite Element computes the element matrix elmat.
Definition: bilininteg.cpp:66
SparseMatrix * H2L_GlobalRestrictionMatrix(FiniteElementSpace *lfes)
Definition: fespace.cpp:494
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
L2Prolongation * B
Backward, fine-to-coarse, operator.
Definition: fespace.hpp:905
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:537
int GetOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:95
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Definition: fespace.cpp:2351
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:116
Table * GetElementDofTable()
Definition: nurbs.hpp:371
void MakeVDimMatrix(SparseMatrix &mat) const
Replicate &#39;mat&#39; in the vector dimension, according to vdim ordering mode.
Definition: fespace.cpp:785
Table * GetBdrElementDofTable()
Definition: nurbs.hpp:372
bool own_mass_integ
Ownership flag for mass_integ.
Definition: fespace.hpp:811
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:28
void GetEdgeInteriorVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:208
void GetEntityDofs(int entity, int index, Array< int > &dofs, Geometry::Type master_geom=Geometry::INVALID) const
Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
Definition: fespace.cpp:608
void ShiftUpI()
Definition: table.cpp:119
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:12621
const CoarseFineTransformations & GetRefinementTransforms()
Definition: mesh.cpp:8094
L2Projection * F
Forward, coarse-to-fine, operator.
Definition: fespace.hpp:904
Operation GetLastOperation() const
Return type of last modification of the mesh.
Definition: mesh.hpp:1175
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:47
bool Nonconforming() const
Definition: fespace.hpp:302
int SizeJ() const
Definition: densemat.hpp:763
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:1659
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:690
virtual const Operator & ForwardOperator()
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
Definition: fespace.cpp:2494
void ConvertToConformingVDofs(const Array< int > &dofs, Array< int > &cdofs)
Definition: fespace.cpp:414
virtual const char * Name() const
Definition: fe_coll.hpp:53
Class for integration point with weight.
Definition: intrules.hpp:25
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:112
bool Parallel() const
Definition: fespace.hpp:732
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:4347
Mesh * mesh
The mesh that FE space lives on (not owned).
Definition: fespace.hpp:93
void ReorderElementToDofTable()
Reorder the scalar DOFs based on the element ordering.
Definition: fespace.cpp:243
GridFunction interpolation operator applicable after mesh refinement.
Definition: fespace.hpp:185
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:799
static const int DimStart[MaxDim+2]
Definition: geom.hpp:47
const NCList & GetNCList(int entity)
Return vertex/edge/face list (entity = 0/1/2, respectively).
Definition: ncmesh.hpp:234
Definition: fespace.hpp:150
int GetNConformingDofs() const
Definition: fespace.cpp:827
void ConnectBoundaries()
Definition: nurbs.cpp:1793
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1855
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition: fespace.cpp:890
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:65
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:635
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.cpp:117
int dim
Definition: ex24.cpp:43
Definition: fespace.hpp:149
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:307
signed char geom
Definition: ncmesh.hpp:156
IsoparametricTransformation Transf
Definition: eltrans.hpp:347
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
Lexicographic ordering for tensor-product FiniteElements.
int parent
Element index in the coarse mesh.
Definition: ncmesh.hpp:49
FiniteElementSpace()
Default constructor: the object is invalid until initialized using the method Load().
Definition: fespace.cpp:59
Create a QuadratureSpace based on the global rules from IntRules.
Definition: fespace.hpp:691
void Constructor(Mesh *mesh, NURBSExtension *ext, const FiniteElementCollection *fec, int vdim=1, int ordering=Ordering::byNODES)
Help function for constructors + Load().
Definition: fespace.cpp:1402
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:192
Geometry::Type GetBdrElementBaseGeometry(int i) const
Definition: mesh.hpp:795
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th face element (2D and 3D).
Definition: fespace.cpp:184
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2457
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th edge.
Definition: fespace.cpp:190
HYPRE_Int * GetTrueDofOffsets() const
Definition: pfespace.hpp:248
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Definition: fespace.cpp:403
virtual void Mult(const Vector &x, Vector &y) const
Perform the L2 projection onto the LOR space.
Definition: fespace.cpp:2685
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:699
void filter_dos(std::string &line)
Definition: text.hpp:41
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:820
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:702
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2648
Vector data type.
Definition: vector.hpp:48
std::tuple< bool, ElementDofOrdering, FaceType, L2FaceValues > key_face
The face restriction operators, see GetFaceRestriction().
Definition: fespace.hpp:134
ID for class HypreParMatrix.
Definition: operator.hpp:233
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:1534
const Table & GetElementToDofTable() const
Definition: fespace.hpp:524
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:60
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.cpp:123
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Definition: fespace.cpp:556
Class representing the storage layout of a QuadratureFunction.
Definition: fespace.hpp:671
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:65
int RowSize(int i) const
Definition: table.hpp:108
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i&#39;th boundary element.
Definition: fespace.cpp:1882
NURBSExtension * NURBSext
Definition: fespace.hpp:117
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:178
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:187
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:315
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:3656
int index
Mesh number.
Definition: ncmesh.hpp:153
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:725
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:376
void Save(std::ostream &out) const
Definition: fespace.cpp:2134
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:838
Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:59
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:834
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:159
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:124
Defines the position of a fine element within a coarse element.
Definition: ncmesh.hpp:46
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:137
int matrix
Index into the DenseTensor corresponding to the parent Geometry::Type stored in CoarseFineTransformat...
Definition: ncmesh.hpp:52
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143