MFEM  v3.4
Finite element discretization library
fespace.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of FiniteElementSpace
13 
14 #include "../general/text.hpp"
15 #include "../mesh/mesh_headers.hpp"
16 #include "fem.hpp"
17 
18 #include <cmath>
19 #include <cstdarg>
20 #include <limits>
21 
22 using namespace std;
23 
24 namespace mfem
25 {
26 
27 template <> void Ordering::
28 DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
29 {
30  // static method
31  int size = dofs.Size();
32  dofs.SetSize(size*vdim);
33  for (int vd = 1; vd < vdim; vd++)
34  {
35  for (int i = 0; i < size; i++)
36  {
37  dofs[i+size*vd] = Map<byNODES>(ndofs, vdim, dofs[i], vd);
38  }
39  }
40 }
41 
42 template <> void Ordering::
43 DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
44 {
45  // static method
46  int size = dofs.Size();
47  dofs.SetSize(size*vdim);
48  for (int vd = vdim-1; vd >= 0; vd--)
49  {
50  for (int i = 0; i < size; i++)
51  {
52  dofs[i+size*vd] = Map<byVDIM>(ndofs, vdim, dofs[i], vd);
53  }
54  }
55 }
56 
57 
58 FiniteElementSpace::FiniteElementSpace()
59  : mesh(NULL), fec(NULL), vdim(0), ordering(Ordering::byNODES),
60  ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0),
61  fdofs(NULL), bdofs(NULL),
62  elem_dof(NULL), bdrElem_dof(NULL),
63  NURBSext(NULL), own_ext(false),
64  cP(NULL), cR(NULL), cP_is_set(false),
65  Th(Operator::ANY_TYPE),
66  sequence(0)
67 { }
68 
70  Mesh *mesh,
71  const FiniteElementCollection *fec)
72 {
73  mesh = mesh ? mesh : orig.mesh;
74  fec = fec ? fec : orig.fec;
75  NURBSExtension *NURBSext = NULL;
76  if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
77  {
78 #ifdef MFEM_USE_MPI
79  ParNURBSExtension *pNURBSext =
80  dynamic_cast<ParNURBSExtension *>(orig.NURBSext);
81  if (pNURBSext)
82  {
83  NURBSext = new ParNURBSExtension(*pNURBSext);
84  }
85  else
86 #endif
87  {
88  NURBSext = new NURBSExtension(*orig.NURBSext);
89  }
90  }
91  Constructor(mesh, NURBSext, fec, orig.vdim, orig.ordering);
92 }
93 
95 {
96  int GeomType = mesh->GetElementBaseGeometry(i);
97  return fec->FiniteElementForGeometry(GeomType)->GetOrder();
98 }
99 
101 {
102  int GeomType = mesh->GetFaceBaseGeometry(i);
103  return fec->FiniteElementForGeometry(GeomType)->GetOrder();
104 }
105 
106 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs, int ndofs) const
107 {
108  if (vdim == 1) { return; }
109  if (ndofs < 0) { ndofs = this->ndofs; }
110 
112  {
113  Ordering::DofsToVDofs<Ordering::byNODES>(ndofs, vdim, dofs);
114  }
115  else
116  {
117  Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs, vdim, dofs);
118  }
119 }
120 
121 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs, int ndofs) const
122 {
123  if (vdim == 1) { return; }
124  if (ndofs < 0) { ndofs = this->ndofs; }
125 
127  {
128  for (int i = 0; i < dofs.Size(); i++)
129  {
130  dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs, vdim, dofs[i], vd);
131  }
132  }
133  else
134  {
135  for (int i = 0; i < dofs.Size(); i++)
136  {
137  dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dofs[i], vd);
138  }
139  }
140 }
141 
142 int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs) const
143 {
144  if (vdim == 1) { return dof; }
145  if (ndofs < 0) { ndofs = this->ndofs; }
146 
148  {
149  return Ordering::Map<Ordering::byNODES>(ndofs, vdim, dof, vd);
150  }
151  else
152  {
153  return Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dof, vd);
154  }
155 }
156 
157 // static function
159 {
160  int n = vdofs.Size(), *vdof = vdofs;
161  for (int i = 0; i < n; i++)
162  {
163  int j;
164  if ((j = vdof[i]) < 0)
165  {
166  vdof[i] = -1-j;
167  }
168  }
169 }
170 
172 {
173  GetElementDofs(i, vdofs);
174  DofsToVDofs(vdofs);
175 }
176 
178 {
179  GetBdrElementDofs(i, vdofs);
180  DofsToVDofs(vdofs);
181 }
182 
184 {
185  GetFaceDofs(i, vdofs);
186  DofsToVDofs(vdofs);
187 }
188 
190 {
191  GetEdgeDofs(i, vdofs);
192  DofsToVDofs(vdofs);
193 }
194 
196 {
197  GetVertexDofs(i, vdofs);
198  DofsToVDofs(vdofs);
199 }
200 
202 {
203  GetElementInteriorDofs(i, vdofs);
204  DofsToVDofs(vdofs);
205 }
206 
208 {
209  GetEdgeInteriorDofs(i, vdofs);
210  DofsToVDofs(vdofs);
211 }
212 
214 {
215  if (elem_dof) { return; }
216 
217  Table *el_dof = new Table;
218  Array<int> dofs;
219  el_dof -> MakeI (mesh -> GetNE());
220  for (int i = 0; i < mesh -> GetNE(); i++)
221  {
222  GetElementDofs (i, dofs);
223  el_dof -> AddColumnsInRow (i, dofs.Size());
224  }
225  el_dof -> MakeJ();
226  for (int i = 0; i < mesh -> GetNE(); i++)
227  {
228  GetElementDofs (i, dofs);
229  el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
230  }
231  el_dof -> ShiftUpI();
232  elem_dof = el_dof;
233 }
234 
236 {
237  delete elem_dof;
238  elem_dof = NULL;
240 }
241 
243 {
244  Array<int> dof_marker(ndofs);
245 
246  dof_marker = -1;
247 
248  int *J = elem_dof->GetJ(), nnz = elem_dof->Size_of_connections();
249  for (int k = 0, dof_counter = 0; k < nnz; k++)
250  {
251  const int sdof = J[k]; // signed dof
252  const int dof = (sdof < 0) ? -1-sdof : sdof;
253  int new_dof = dof_marker[dof];
254  if (new_dof < 0)
255  {
256  dof_marker[dof] = new_dof = dof_counter++;
257  }
258  J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
259  }
260 }
261 
263 {
264  if (dof_elem_array.Size()) { return; }
265 
267 
270  dof_elem_array = -1;
271  for (int i = 0; i < mesh -> GetNE(); i++)
272  {
273  const int *dofs = elem_dof -> GetRow(i);
274  const int n = elem_dof -> RowSize(i);
275  for (int j = 0; j < n; j++)
276  {
277  if (dof_elem_array[dofs[j]] < 0)
278  {
279  dof_elem_array[dofs[j]] = i;
280  dof_ldof_array[dofs[j]] = j;
281  }
282  }
283  }
284 }
285 
286 static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
287 {
288  for (int i = 0; i < dofs.Size(); i++)
289  {
290  int k = dofs[i];
291  if (k < 0) { k = -1 - k; }
292  mark_array[k] = -1;
293  }
294 }
295 
297  Array<int> &ess_vdofs,
298  int component) const
299 {
300  Array<int> vdofs, dofs;
301 
302  ess_vdofs.SetSize(GetVSize());
303  ess_vdofs = 0;
304 
305  for (int i = 0; i < GetNBE(); i++)
306  {
307  if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
308  {
309  if (component < 0)
310  {
311  // Mark all components.
312  GetBdrElementVDofs(i, vdofs);
313  mark_dofs(vdofs, ess_vdofs);
314  }
315  else
316  {
317  GetBdrElementDofs(i, dofs);
318  for (int d = 0; d < dofs.Size(); d++)
319  { dofs[d] = DofToVDof(dofs[d], component); }
320  mark_dofs(dofs, ess_vdofs);
321  }
322  }
323  }
324 
325  // mark possible hidden boundary edges in a non-conforming mesh, also
326  // local DOFs affected by boundary elements on other processors
327  if (mesh->ncmesh)
328  {
329  Array<int> bdr_verts, bdr_edges;
330  mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges);
331 
332  for (int i = 0; i < bdr_verts.Size(); i++)
333  {
334  if (component < 0)
335  {
336  GetVertexVDofs(bdr_verts[i], vdofs);
337  mark_dofs(vdofs, ess_vdofs);
338  }
339  else
340  {
341  GetVertexDofs(bdr_verts[i], dofs);
342  for (int d = 0; d < dofs.Size(); d++)
343  { dofs[d] = DofToVDof(dofs[d], component); }
344  mark_dofs(dofs, ess_vdofs);
345  }
346  }
347  for (int i = 0; i < bdr_edges.Size(); i++)
348  {
349  if (component < 0)
350  {
351  GetEdgeVDofs(bdr_edges[i], vdofs);
352  mark_dofs(vdofs, ess_vdofs);
353  }
354  else
355  {
356  GetEdgeDofs(bdr_edges[i], dofs);
357  for (int d = 0; d < dofs.Size(); d++)
358  { dofs[d] = DofToVDof(dofs[d], component); }
359  mark_dofs(dofs, ess_vdofs);
360  }
361  }
362  }
363 }
364 
366  Array<int> &ess_tdof_list,
367  int component)
368 {
369  Array<int> ess_vdofs, ess_tdofs;
370  GetEssentialVDofs(bdr_attr_is_ess, ess_vdofs, component);
372  if (!R)
373  {
374  ess_tdofs.MakeRef(ess_vdofs);
375  }
376  else
377  {
378  R->BooleanMult(ess_vdofs, ess_tdofs);
379  }
380  MarkerToList(ess_tdofs, ess_tdof_list);
381 }
382 
383 // static method
385  Array<int> &list)
386 {
387  int num_marked = 0;
388  for (int i = 0; i < marker.Size(); i++)
389  {
390  if (marker[i]) { num_marked++; }
391  }
392  list.SetSize(0);
393  list.Reserve(num_marked);
394  for (int i = 0; i < marker.Size(); i++)
395  {
396  if (marker[i]) { list.Append(i); }
397  }
398 }
399 
400 // static method
401 void FiniteElementSpace::ListToMarker(const Array<int> &list, int marker_size,
402  Array<int> &marker, int mark_val)
403 {
404  marker.SetSize(marker_size);
405  marker = 0;
406  for (int i = 0; i < list.Size(); i++)
407  {
408  marker[list[i]] = mark_val;
409  }
410 }
411 
413  Array<int> &cdofs)
414 {
416  if (cP) { cP->BooleanMultTranspose(dofs, cdofs); }
417  else { dofs.Copy(cdofs); }
418 }
419 
421  Array<int> &dofs)
422 {
424  if (cR) { cR->BooleanMultTranspose(cdofs, dofs); }
425  else { cdofs.Copy(dofs); }
426 }
427 
428 SparseMatrix *
430 {
431  int i, j;
432  Array<int> d_vdofs, c_vdofs;
433  SparseMatrix *R;
434 
435  R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
436 
437  for (i = 0; i < mesh -> GetNE(); i++)
438  {
439  this -> GetElementVDofs (i, d_vdofs);
440  cfes -> GetElementVDofs (i, c_vdofs);
441 
442 #ifdef MFEM_DEBUG
443  if (d_vdofs.Size() != c_vdofs.Size())
444  {
445  mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
446  }
447 #endif
448 
449  for (j = 0; j < d_vdofs.Size(); j++)
450  {
451  R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
452  }
453  }
454 
455  R -> Finalize();
456 
457  return R;
458 }
459 
460 SparseMatrix *
462 {
463  int i, j;
464  Array<int> d_dofs, c_dofs;
465  SparseMatrix *R;
466 
467  R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
468 
469  for (i = 0; i < mesh -> GetNE(); i++)
470  {
471  this -> GetElementDofs (i, d_dofs);
472  cfes -> GetElementDofs (i, c_dofs);
473 
474 #ifdef MFEM_DEBUG
475  if (c_dofs.Size() != 1)
476  mfem_error ("FiniteElementSpace::"
477  "D2Const_GlobalRestrictionMatrix (...)");
478 #endif
479 
480  for (j = 0; j < d_dofs.Size(); j++)
481  {
482  R -> Set (c_dofs[0], d_dofs[j], 1.0);
483  }
484  }
485 
486  R -> Finalize();
487 
488  return R;
489 }
490 
491 SparseMatrix *
493 {
494  SparseMatrix *R;
495  DenseMatrix loc_restr;
496  Array<int> l_dofs, h_dofs;
497 
498  R = new SparseMatrix (lfes -> GetNDofs(), ndofs);
499 
500  if (!lfes->GetNE())
501  {
502  R->Finalize();
503  return R;
504  }
505 
506  const FiniteElement *h_fe = this -> GetFE (0);
507  const FiniteElement *l_fe = lfes -> GetFE (0);
510  h_fe->Project(*l_fe, T, loc_restr);
511 
512  for (int i = 0; i < mesh -> GetNE(); i++)
513  {
514  this -> GetElementDofs (i, h_dofs);
515  lfes -> GetElementDofs (i, l_dofs);
516 
517  R -> SetSubMatrix (l_dofs, h_dofs, loc_restr, 1);
518  }
519 
520  R -> Finalize();
521 
522  return R;
523 }
524 
525 void
527  Array<int>& slave_dofs, DenseMatrix& I)
528 {
529  for (int i = 0; i < slave_dofs.Size(); i++)
530  {
531  int sdof = slave_dofs[i];
532  if (!deps.RowSize(sdof)) // not processed yet?
533  {
534  for (int j = 0; j < master_dofs.Size(); j++)
535  {
536  double coef = I(i, j);
537  if (std::abs(coef) > 1e-12)
538  {
539  int mdof = master_dofs[j];
540  if (mdof != sdof && mdof != (-1-sdof))
541  {
542  deps.Add(sdof, mdof, coef);
543  }
544  }
545  }
546  }
547  }
548 }
549 
550 bool FiniteElementSpace::DofFinalizable(int dof, const Array<bool>& finalized,
551  const SparseMatrix& deps)
552 {
553  const int* dep = deps.GetRowColumns(dof);
554  int ndep = deps.RowSize(dof);
555 
556  // are all constraining DOFs finalized?
557  for (int i = 0; i < ndep; i++)
558  {
559  if (!finalized[dep[i]]) { return false; }
560  }
561  return true;
562 }
563 
564 void
565 FiniteElementSpace::GetEntityDofs(int entity, int index, Array<int> &dofs) const
566 {
567  switch (entity)
568  {
569  case 0: GetVertexDofs(index, dofs); break;
570  case 1: GetEdgeDofs(index, dofs); break;
571  case 2: GetFaceDofs(index, dofs); break;
572  }
573 }
574 
575 
577 {
578 #ifdef MFEM_USE_MPI
579  MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
580  "This method should not be used with a ParFiniteElementSpace!");
581 #endif
582 
583  if (cP_is_set) { return; }
584  cP_is_set = true;
585 
586  // For each slave DOF, the dependency matrix will contain a row that
587  // expresses the slave DOF as a linear combination of its immediate master
588  // DOFs. Rows of independent DOFs will remain empty.
589  SparseMatrix deps(ndofs);
590 
591  // collect local edge/face dependencies
592  for (int entity = 1; entity <= 2; entity++)
593  {
594  const NCMesh::NCList &list = (entity > 1) ? mesh->ncmesh->GetFaceList()
595  /* */ : mesh->ncmesh->GetEdgeList();
596  if (!list.masters.size()) { continue; }
597 
599  if (entity > 1) { T.SetFE(&QuadrilateralFE); }
600  else { T.SetFE(&SegmentFE); }
601 
602  int geom = (entity > 1) ? Geometry::SQUARE : Geometry::SEGMENT;
604  if (!fe) { continue; }
605 
606  Array<int> master_dofs, slave_dofs;
607  DenseMatrix I(fe->GetDof());
608 
609  // loop through all master edges/faces, constrain their slave edges/faces
610  for (unsigned mi = 0; mi < list.masters.size(); mi++)
611  {
612  const NCMesh::Master &master = list.masters[mi];
613  GetEntityDofs(entity, master.index, master_dofs);
614  if (!master_dofs.Size()) { continue; }
615 
616  for (int si = master.slaves_begin; si < master.slaves_end; si++)
617  {
618  const NCMesh::Slave &slave = list.slaves[si];
619  GetEntityDofs(entity, slave.index, slave_dofs);
620  if (!slave_dofs.Size()) { continue; }
621 
622  slave.OrientedPointMatrix(T.GetPointMat());
624  fe->GetLocalInterpolation(T, I);
625 
626  // make each slave DOF dependent on all master DOFs
627  AddDependencies(deps, master_dofs, slave_dofs, I);
628  }
629  }
630  }
631 
632  deps.Finalize();
633 
634  // DOFs that stayed independent are true DOFs
635  int n_true_dofs = 0;
636  for (int i = 0; i < ndofs; i++)
637  {
638  if (!deps.RowSize(i)) { n_true_dofs++; }
639  }
640 
641  // if all dofs are true dofs leave cP and cR NULL
642  if (n_true_dofs == ndofs)
643  {
644  cP = cR = NULL; // will be treated as identities
645  return;
646  }
647 
648  // create the conforming restriction matrix cR
649  int *cR_J;
650  {
651  int *cR_I = new int[n_true_dofs+1];
652  double *cR_A = new double[n_true_dofs];
653  cR_J = new int[n_true_dofs];
654  for (int i = 0; i < n_true_dofs; i++)
655  {
656  cR_I[i] = i;
657  cR_A[i] = 1.0;
658  }
659  cR_I[n_true_dofs] = n_true_dofs;
660  cR = new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, ndofs);
661  }
662 
663  // create the conforming prolongation matrix cP
664  cP = new SparseMatrix(ndofs, n_true_dofs);
665 
666  Array<bool> finalized(ndofs);
667  finalized = false;
668 
669  // put identity in the restriction and prolongation matrices for true DOFs
670  for (int i = 0, true_dof = 0; i < ndofs; i++)
671  {
672  if (!deps.RowSize(i))
673  {
674  cR_J[true_dof] = i;
675  cP->Add(i, true_dof++, 1.0);
676  finalized[i] = true;
677  }
678  }
679 
680  // Now calculate cP rows of slave DOFs as combinations of cP rows of their
681  // master DOFs. It is possible that some slave DOFs depend on DOFs that are
682  // themselves slaves. Here we resolve such indirect constraints by first
683  // calculating rows of the cP matrix for DOFs whose master DOF cP rows are
684  // already known (in the first iteration these are the true DOFs). In the
685  // second iteration, slaves of slaves can be 'finalized' (given a row in the
686  // cP matrix), in the third iteration slaves of slaves of slaves, etc.
687  bool finished;
688  int n_finalized = n_true_dofs;
689  Array<int> cols;
690  Vector srow;
691  do
692  {
693  finished = true;
694  for (int dof = 0; dof < ndofs; dof++)
695  {
696  if (!finalized[dof] && DofFinalizable(dof, finalized, deps))
697  {
698  const int* dep_col = deps.GetRowColumns(dof);
699  const double* dep_coef = deps.GetRowEntries(dof);
700  int n_dep = deps.RowSize(dof);
701 
702  for (int j = 0; j < n_dep; j++)
703  {
704  cP->GetRow(dep_col[j], cols, srow);
705  srow *= dep_coef[j];
706  cP->AddRow(dof, cols, srow);
707  }
708 
709  finalized[dof] = true;
710  n_finalized++;
711  finished = false;
712  }
713  }
714  }
715  while (!finished);
716 
717  // if everything is consistent (mesh, face orientations, etc.), we should
718  // be able to finalize all slave DOFs, otherwise it's a serious error
719  if (n_finalized != ndofs)
720  {
721  MFEM_ABORT("Error creating cP matrix.");
722  }
723 
724  cP->Finalize();
725 
726  if (vdim > 1)
727  {
728  MakeVDimMatrix(*cP);
729  MakeVDimMatrix(*cR);
730  }
731 }
732 
734 {
735  if (vdim == 1) { return; }
736 
737  int height = mat.Height();
738  int width = mat.Width();
739 
740  SparseMatrix *vmat = new SparseMatrix(vdim*height, vdim*width);
741 
742  Array<int> dofs, vdofs;
743  Vector srow;
744  for (int i = 0; i < height; i++)
745  {
746  mat.GetRow(i, dofs, srow);
747  for (int vd = 0; vd < vdim; vd++)
748  {
749  dofs.Copy(vdofs);
750  DofsToVDofs(vd, vdofs, width);
751  vmat->SetRow(DofToVDof(i, vd, height), vdofs, srow);
752  }
753  }
754  vmat->Finalize();
755 
756  mat.Swap(*vmat);
757  delete vmat;
758 }
759 
760 
762 {
763  if (Conforming()) { return NULL; }
765  return cP;
766 }
767 
769 {
770  if (Conforming()) { return NULL; }
772  return cR;
773 }
774 
776 {
778  return P ? (P->Width() / vdim) : ndofs;
779 }
780 
782  const int coarse_ndofs, const Table &coarse_elem_dof,
783  const DenseTensor &localP) const
784 {
785  MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
786 
787  Array<int> dofs, coarse_dofs, coarse_vdofs;
788  Vector row;
789 
790  const int coarse_ldof = localP.SizeJ();
791  const int fine_ldof = localP.SizeI();
792  SparseMatrix *P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim,
793  coarse_ldof);
794 
795  Array<int> mark(P->Height());
796  mark = 0;
797 
799 
800  for (int k = 0; k < mesh->GetNE(); k++)
801  {
802  const Embedding &emb = rtrans.embeddings[k];
803  const DenseMatrix &lP = localP(emb.matrix);
804 
805  elem_dof->GetRow(k, dofs);
806  coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
807 
808  for (int vd = 0; vd < vdim; vd++)
809  {
810  coarse_dofs.Copy(coarse_vdofs);
811  DofsToVDofs(vd, coarse_vdofs, coarse_ndofs);
812 
813  for (int i = 0; i < fine_ldof; i++)
814  {
815  int r = DofToVDof(dofs[i], vd);
816  int m = (r >= 0) ? r : (-1 - r);
817 
818  if (!mark[m])
819  {
820  lP.GetRow(i, row);
821  P->SetRow(r, coarse_vdofs, row);
822  mark[m] = 1;
823  }
824  }
825  }
826  }
827 
828  MFEM_ASSERT(mark.Sum() == P->Height(), "Not all rows of P set.");
829  return P;
830 }
831 
833 {
834  int geom = mesh->GetElementBaseGeometry(); // assuming the same geom
836 
838 
839  int nmat = rtrans.point_matrices.SizeK();
840  int ldof = fe->GetDof(); // assuming the same FE everywhere
841 
844 
845  // calculate local interpolation matrices for all refinement types
846  localP.SetSize(ldof, ldof, nmat);
847  for (int i = 0; i < nmat; i++)
848  {
849  isotr.GetPointMat() = rtrans.point_matrices(i);
850  isotr.FinalizeTransformation();
851  fe->GetLocalInterpolation(isotr, localP(i));
852  }
853 }
854 
856  const Table* old_elem_dof)
857 {
858  MFEM_VERIFY(ndofs >= old_ndofs, "Previous space is not coarser.");
859 
860  DenseTensor localP;
862 
863  return RefinementMatrix_main(old_ndofs, *old_elem_dof, localP);
864 }
865 
867 (const FiniteElementSpace* fespace, Table* old_elem_dof, int old_ndofs)
868  : fespace(fespace)
869  , old_elem_dof(old_elem_dof)
870 {
871  MFEM_VERIFY(fespace->GetNDofs() >= old_ndofs,
872  "Previous space is not coarser.");
873 
874  width = old_ndofs * fespace->GetVDim();
875  height = fespace->GetVSize();
876 
877  fespace->GetLocalRefinementMatrices(localP);
878 }
879 
881  const FiniteElementSpace *fespace, const FiniteElementSpace *coarse_fes)
882  : Operator(fespace->GetVSize(), coarse_fes->GetVSize()),
883  fespace(fespace), old_elem_dof(NULL)
884 {
885  fespace->GetLocalRefinementMatrices(*coarse_fes, localP);
886  // Make a copy of the coarse elem_dof Table.
887  old_elem_dof = new Table(coarse_fes->GetElementToDofTable());
888 }
889 
891 {
892  delete old_elem_dof;
893 }
894 
896 ::Mult(const Vector &x, Vector &y) const
897 {
898  Mesh* mesh = fespace->GetMesh();
900 
901  Array<int> dofs, old_dofs, old_vdofs;
902 
903  Array<char> processed(fespace->GetVSize());
904  processed = 0;
905 
906  int vdim = fespace->GetVDim();
907  int old_ndofs = width / vdim;
908 
909  for (int k = 0; k < mesh->GetNE(); k++)
910  {
911  const Embedding &emb = rtrans.embeddings[k];
912  const DenseMatrix &lP = localP(emb.matrix);
913 
914  fespace->GetElementDofs(k, dofs);
915  old_elem_dof->GetRow(emb.parent, old_dofs);
916 
917  for (int vd = 0; vd < vdim; vd++)
918  {
919  old_dofs.Copy(old_vdofs);
920  fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
921 
922  for (int i = 0; i < dofs.Size(); i++)
923  {
924  double rsign, osign;
925  int r = fespace->DofToVDof(dofs[i], vd);
926  r = DecodeDof(r, rsign);
927 
928  if (!processed[r])
929  {
930  double value = 0.0;
931  for (int j = 0; j < old_vdofs.Size(); j++)
932  {
933  int o = DecodeDof(old_vdofs[j], osign);
934  value += x[o] * lP(i, j) * osign;
935  }
936  y[r] = value * rsign;
937  processed[r] = 1;
938  }
939  }
940  }
941  }
942 }
943 
945  const DenseMatrix &invdfdx,
946  const IntegrationPoint &pt, Vector &x)
947 {
948  // invert a linear transform with one Newton step
949  IntegrationPoint p0;
950  p0.Set3(0, 0, 0);
951  trans.Transform(p0, x);
952 
953  double store[3];
954  Vector v(store, x.Size());
955  pt.Get(v, x.Size());
956  v -= x;
957 
958  invdfdx.Mult(v, x);
959 }
960 
962 {
963  int geom = mesh->GetElementBaseGeometry(); // assuming the same geom
965  const IntegrationRule &nodes = fe->GetNodes();
966 
967  const CoarseFineTransformations &dtrans =
969 
970  int nmat = dtrans.point_matrices.SizeK();
971  int ldof = fe->GetDof();
972  int dim = mesh->Dimension();
973 
974  LinearFECollection linfec;
976  isotr.SetFE(linfec.FiniteElementForGeometry(geom));
977 
978  DenseMatrix invdfdx(dim);
979  IntegrationPoint ipt;
980  Vector pt(&ipt.x, dim), shape(ldof);
981 
982  // calculate local restriction matrices for all refinement types
983  localR.SetSize(ldof, ldof, nmat);
984  for (int i = 0; i < nmat; i++)
985  {
986  DenseMatrix &lR = localR(i);
987  lR = infinity(); // marks invalid rows
988 
989  isotr.GetPointMat() = dtrans.point_matrices(i);
990  isotr.FinalizeTransformation();
991  isotr.SetIntPoint(&nodes[0]);
992  CalcInverse(isotr.Jacobian(), invdfdx);
993 
994  for (int j = 0; j < nodes.Size(); j++)
995  {
996  InvertLinearTrans(isotr, invdfdx, nodes[j], pt);
997  if (Geometries.CheckPoint(geom, ipt)) // do we need an epsilon here?
998  {
999  IntegrationPoint ip;
1000  ip.Set(pt, dim);
1001  MFEM_ASSERT(dynamic_cast<const NodalFiniteElement*>(fe),
1002  "only nodal FEs are implemented");
1003  fe->CalcShape(ip, shape); // TODO: H(curl), etc.?
1004  lR.SetRow(j, shape);
1005  }
1006  }
1007  lR.Threshold(1e-12);
1008  }
1009 }
1010 
1012  const Table* old_elem_dof)
1013 {
1014  MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
1015  MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
1016  MFEM_VERIFY(ndofs <= old_ndofs, "Previous space is not finer.");
1017 
1018  Array<int> dofs, old_dofs, old_vdofs;
1019  Vector row;
1020 
1021  DenseTensor localR;
1023 
1024  SparseMatrix *R = new SparseMatrix(ndofs*vdim, old_ndofs*vdim,
1025  localR.SizeI());
1026  Array<int> mark(R->Height());
1027  mark = 0;
1028 
1029  const CoarseFineTransformations &dtrans =
1031 
1032  for (int k = 0; k < dtrans.embeddings.Size(); k++)
1033  {
1034  const Embedding &emb = dtrans.embeddings[k];
1035  DenseMatrix &lR = localR(emb.matrix);
1036 
1037  elem_dof->GetRow(emb.parent, dofs);
1038  old_elem_dof->GetRow(k, old_dofs);
1039 
1040  for (int vd = 0; vd < vdim; vd++)
1041  {
1042  old_dofs.Copy(old_vdofs);
1043  DofsToVDofs(vd, old_vdofs, old_ndofs);
1044 
1045  for (int i = 0; i < lR.Height(); i++)
1046  {
1047  if (lR(i, 0) == infinity()) { continue; }
1048 
1049  int r = DofToVDof(dofs[i], vd);
1050  int m = (r >= 0) ? r : (-1 - r);
1051 
1052  if (!mark[m])
1053  {
1054  lR.GetRow(i, row);
1055  R->SetRow(r, old_vdofs, row);
1056  mark[m] = 1;
1057  }
1058  }
1059  }
1060  }
1061 
1062  MFEM_ASSERT(mark.Sum() == R->Height(), "Not all rows of R set.");
1063  return R;
1064 }
1065 
1067  const FiniteElementSpace &coarse_fes, DenseTensor &localP) const
1068 {
1069  // Assumptions: see the declaration of the method.
1070 
1071  int fine_geom = mesh->GetElementBaseGeometry(0);
1072  const FiniteElement *fine_fe = fec->FiniteElementForGeometry(fine_geom);
1073  const FiniteElement *coarse_fe = coarse_fes.GetFE(0);
1074 
1076 
1077  int nmat = rtrans.point_matrices.SizeK();
1078 
1080  isotr.SetIdentityTransformation(fine_geom);
1081 
1082  // Calculate the local interpolation matrices for all refinement types
1083  localP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), nmat);
1084  for (int i = 0; i < nmat; i++)
1085  {
1086  isotr.GetPointMat() = rtrans.point_matrices(i);
1087  isotr.FinalizeTransformation();
1088  fine_fe->GetTransferMatrix(*coarse_fe, isotr, localP(i));
1089  }
1090 }
1091 
1094  int vdim, int ordering)
1095 {
1096  this->mesh = mesh;
1097  this->fec = fec;
1098  this->vdim = vdim;
1099  this->ordering = (Ordering::Type) ordering;
1100 
1101  elem_dof = NULL;
1102  sequence = mesh->GetSequence();
1104 
1105  const NURBSFECollection *nurbs_fec =
1106  dynamic_cast<const NURBSFECollection *>(fec);
1107  if (nurbs_fec)
1108  {
1109  if (!mesh->NURBSext)
1110  {
1111  mfem_error("FiniteElementSpace::FiniteElementSpace :\n"
1112  " NURBS FE space requires NURBS mesh.");
1113  }
1114 
1115  if (NURBSext == NULL)
1116  {
1117  this->NURBSext = mesh->NURBSext;
1118  own_ext = 0;
1119  }
1120  else
1121  {
1122  this->NURBSext = NURBSext;
1123  own_ext = 1;
1124  }
1125  UpdateNURBS();
1126  cP = cR = NULL;
1127  cP_is_set = false;
1128  }
1129  else
1130  {
1131  this->NURBSext = NULL;
1132  own_ext = 0;
1133  Construct();
1134  }
1136 }
1137 
1139 {
1140  if (NURBSext && !own_ext)
1141  {
1142  mfem_error("FiniteElementSpace::StealNURBSext");
1143  }
1144  own_ext = 0;
1145 
1146  return NURBSext;
1147 }
1148 
1150 {
1151  nvdofs = 0;
1152  nedofs = 0;
1153  nfdofs = 0;
1154  nbdofs = 0;
1155  fdofs = NULL;
1156  bdofs = NULL;
1157 
1158  dynamic_cast<const NURBSFECollection *>(fec)->Reset();
1159 
1160  ndofs = NURBSext->GetNDof();
1163 }
1164 
1166 {
1167  // This method should be used only for non-NURBS spaces.
1168  MFEM_ASSERT(!NURBSext, "internal error");
1169 
1170  elem_dof = NULL;
1171  bdrElem_dof = NULL;
1172 
1174 
1175  if ( mesh->Dimension() > 1 )
1176  {
1178  }
1179  else
1180  {
1181  nedofs = 0;
1182  }
1183 
1184  ndofs = 0;
1185  nfdofs = 0;
1186  nbdofs = 0;
1187  bdofs = NULL;
1188  fdofs = NULL;
1189  cP = NULL;
1190  cR = NULL;
1191  cP_is_set = false;
1192  // Th is initialized/destroyed before this method is called.
1193 
1194  if (mesh->Dimension() == 3 && mesh->GetNE())
1195  {
1196  // Here we assume that all faces in the mesh have the same base
1197  // geometry -- the base geometry of the 0-th face element.
1198  // The class Mesh assumes the same inside GetFaceBaseGeometry(...).
1199  // Thus we do not need to generate all the faces in the mesh
1200  // if we do not need them.
1201  int fdof = fec->DofForGeometry(mesh->GetFaceBaseGeometry(0));
1202  if (fdof > 0)
1203  {
1204  fdofs = new int[mesh->GetNFaces()+1];
1205  fdofs[0] = 0;
1206  for (int i = 0; i < mesh->GetNFaces(); i++)
1207  {
1208  nfdofs += fdof;
1209  // nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(i));
1210  fdofs[i+1] = nfdofs;
1211  }
1212  }
1213  }
1214 
1215  if (mesh->Dimension() > 0)
1216  {
1217  bdofs = new int[mesh->GetNE()+1];
1218  bdofs[0] = 0;
1219  for (int i = 0; i < mesh->GetNE(); i++)
1220  {
1221  int geom = mesh->GetElementBaseGeometry(i);
1223  bdofs[i+1] = nbdofs;
1224  }
1225  }
1226 
1227  ndofs = nvdofs + nedofs + nfdofs + nbdofs;
1228 
1229  // Do not build elem_dof Table here: in parallel it has to be constructed
1230  // later.
1231 }
1232 
1234 {
1235  if (elem_dof)
1236  {
1237  elem_dof -> GetRow (i, dofs);
1238  }
1239  else
1240  {
1241  Array<int> V, E, Eo, F, Fo;
1242  int k, j, nv, ne, nf, nb, nfd, nd;
1243  int *ind, dim;
1244 
1245  dim = mesh->Dimension();
1247  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
1248  nb = (dim > 0) ? fec->DofForGeometry(mesh->GetElementBaseGeometry(i)) : 0;
1249  if (nv > 0)
1250  {
1251  mesh->GetElementVertices(i, V);
1252  }
1253  if (ne > 0)
1254  {
1255  mesh->GetElementEdges(i, E, Eo);
1256  }
1257  nfd = 0;
1258  if (dim == 3)
1259  {
1261  {
1262  mesh->GetElementFaces(i, F, Fo);
1263  for (k = 0; k < F.Size(); k++)
1264  {
1265  nfd += fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
1266  }
1267  }
1268  }
1269  nd = V.Size() * nv + E.Size() * ne + nfd + nb;
1270  dofs.SetSize(nd);
1271  if (nv > 0)
1272  {
1273  for (k = 0; k < V.Size(); k++)
1274  {
1275  for (j = 0; j < nv; j++)
1276  {
1277  dofs[k*nv+j] = V[k]*nv+j;
1278  }
1279  }
1280  nv *= V.Size();
1281  }
1282  if (ne > 0)
1283  {
1284  // if (dim > 1)
1285  for (k = 0; k < E.Size(); k++)
1286  {
1288  for (j = 0; j < ne; j++)
1289  {
1290  if (ind[j] < 0)
1291  {
1292  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1293  }
1294  else
1295  {
1296  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1297  }
1298  }
1299  }
1300  }
1301  ne = nv + ne * E.Size();
1302  if (nfd > 0)
1303  // if (dim == 3)
1304  {
1305  for (k = 0; k < F.Size(); k++)
1306  {
1308  Fo[k]);
1310  for (j = 0; j < nf; j++)
1311  {
1312  if (ind[j] < 0)
1313  {
1314  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[F[k]]+(-1-ind[j]) );
1315  }
1316  else
1317  {
1318  dofs[ne+j] = nvdofs+nedofs+fdofs[F[k]]+ind[j];
1319  }
1320  }
1321  ne += nf;
1322  }
1323  }
1324  if (nb > 0)
1325  {
1326  k = nvdofs + nedofs + nfdofs + bdofs[i];
1327  for (j = 0; j < nb; j++)
1328  {
1329  dofs[ne+j] = k + j;
1330  }
1331  }
1332  }
1333 }
1334 
1336 {
1337  const FiniteElement *FE =
1339 
1340  if (NURBSext)
1341  {
1342  NURBSext->LoadFE(i, FE);
1343  }
1344 
1345  return FE;
1346 }
1347 
1349 {
1350  if (bdrElem_dof)
1351  {
1352  bdrElem_dof->GetRow(i, dofs);
1353  }
1354  else
1355  {
1356  Array<int> V, E, Eo;
1357  int k, j, nv, ne, nf, nd, iF, oF;
1358  int *ind, dim;
1359 
1360  dim = mesh->Dimension();
1362  if (nv > 0)
1363  {
1364  mesh->GetBdrElementVertices(i, V);
1365  }
1366  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
1367  if (ne > 0)
1368  {
1369  mesh->GetBdrElementEdges(i, E, Eo);
1370  }
1371  nd = V.Size() * nv + E.Size() * ne;
1372  nf = (dim == 3) ? (fec->DofForGeometry(
1373  mesh->GetBdrElementBaseGeometry(i))) : (0);
1374  if (nf > 0)
1375  {
1376  nd += nf;
1377  mesh->GetBdrElementFace(i, &iF, &oF);
1378  }
1379  dofs.SetSize(nd);
1380  if (nv > 0)
1381  {
1382  for (k = 0; k < V.Size(); k++)
1383  {
1384  for (j = 0; j < nv; j++)
1385  {
1386  dofs[k*nv+j] = V[k]*nv+j;
1387  }
1388  }
1389  nv *= V.Size();
1390  }
1391  if (ne > 0)
1392  {
1393  // if (dim > 1)
1394  for (k = 0; k < E.Size(); k++)
1395  {
1397  for (j = 0; j < ne; j++)
1398  {
1399  if (ind[j] < 0)
1400  {
1401  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1402  }
1403  else
1404  {
1405  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1406  }
1407  }
1408  }
1409  }
1410  if (nf > 0)
1411  // if (dim == 3)
1412  {
1413  ne = nv + ne * E.Size();
1414  ind = (fec->DofOrderForOrientation(
1415  mesh->GetBdrElementBaseGeometry(i), oF));
1416  for (j = 0; j < nf; j++)
1417  {
1418  if (ind[j] < 0)
1419  {
1420  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[iF]+(-1-ind[j]) );
1421  }
1422  else
1423  {
1424  dofs[ne+j] = nvdofs+nedofs+fdofs[iF]+ind[j];
1425  }
1426  }
1427  }
1428  }
1429 }
1430 
1432 {
1433  int j, k, nv, ne, nf, nd, dim = mesh->Dimension();
1434  Array<int> V, E, Eo;
1435  const int *ind;
1436 
1437  // for 1D, 2D and 3D faces
1439  ne = (dim > 1) ? fec->DofForGeometry(Geometry::SEGMENT) : 0;
1440  if (nv > 0)
1441  {
1442  mesh->GetFaceVertices(i, V);
1443  }
1444  if (ne > 0)
1445  {
1446  mesh->GetFaceEdges(i, E, Eo);
1447  }
1448  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
1449  nd = V.Size() * nv + E.Size() * ne + nf;
1450  dofs.SetSize(nd);
1451  if (nv > 0)
1452  {
1453  for (k = 0; k < V.Size(); k++)
1454  {
1455  for (j = 0; j < nv; j++)
1456  {
1457  dofs[k*nv+j] = V[k]*nv+j;
1458  }
1459  }
1460  }
1461  nv *= V.Size();
1462  if (ne > 0)
1463  {
1464  for (k = 0; k < E.Size(); k++)
1465  {
1467  for (j = 0; j < ne; j++)
1468  {
1469  if (ind[j] < 0)
1470  {
1471  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1472  }
1473  else
1474  {
1475  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1476  }
1477  }
1478  }
1479  }
1480  ne = nv + ne * E.Size();
1481  if (nf > 0)
1482  {
1483  for (j = nvdofs+nedofs+fdofs[i], k = 0; k < nf; j++, k++)
1484  {
1485  dofs[ne+k] = j;
1486  }
1487  }
1488 }
1489 
1491 {
1492  int j, k, nv, ne;
1493  Array<int> V;
1494 
1496  if (nv > 0)
1497  {
1498  mesh->GetEdgeVertices(i, V);
1499  }
1501  dofs.SetSize(2*nv+ne);
1502  if (nv > 0)
1503  {
1504  for (k = 0; k < 2; k++)
1505  {
1506  for (j = 0; j < nv; j++)
1507  {
1508  dofs[k*nv+j] = V[k]*nv+j;
1509  }
1510  }
1511  }
1512  nv *= 2;
1513  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
1514  {
1515  dofs[nv+j] = k;
1516  }
1517 }
1518 
1520 {
1521  int j, nv;
1522 
1524  dofs.SetSize(nv);
1525  for (j = 0; j < nv; j++)
1526  {
1527  dofs[j] = i*nv+j;
1528  }
1529 }
1530 
1532 {
1533  int j, k, nb;
1534  if (mesh->Dimension() == 0) { dofs.SetSize(0); return; }
1535  nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (i));
1536  dofs.SetSize (nb);
1537  k = nvdofs + nedofs + nfdofs + bdofs[i];
1538  for (j = 0; j < nb; j++)
1539  {
1540  dofs[j] = k + j;
1541  }
1542 }
1543 
1545 {
1546  int j, k, ne;
1547 
1548  ne = fec -> DofForGeometry (Geometry::SEGMENT);
1549  dofs.SetSize (ne);
1550  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
1551  {
1552  dofs[j] = k;
1553  }
1554 }
1555 
1557 {
1558  int j, k, nf;
1559 
1560  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
1561  dofs.SetSize (nf);
1562  if (nf > 0)
1563  {
1564  for (j = 0, k = nvdofs+nedofs+fdofs[i]; j < nf; j++, k++)
1565  {
1566  dofs[j] = k;
1567  }
1568  }
1569 }
1570 
1572 {
1573  const FiniteElement *BE;
1574 
1575  switch ( mesh->Dimension() )
1576  {
1577  case 1:
1579  break;
1580  case 2:
1582  break;
1583  case 3:
1584  default:
1587  }
1588 
1589  if (NURBSext)
1590  {
1591  NURBSext->LoadBE(i, BE);
1592  }
1593 
1594  return BE;
1595 }
1596 
1598 {
1599  const FiniteElement *fe;
1600 
1601  switch (mesh->Dimension())
1602  {
1603  case 1:
1605  break;
1606  case 2:
1608  break;
1609  case 3:
1610  default:
1612  }
1613 
1614  // if (NURBSext)
1615  // NURBSext->LoadFaceElement(i, fe);
1616 
1617  return fe;
1618 }
1619 
1621 {
1623 }
1624 
1626  int i, int geom_type) const
1627 {
1628  return fec->TraceFiniteElementForGeometry(geom_type);
1629 }
1630 
1632 {
1633  Destroy();
1634 }
1635 
1637 {
1638  delete cR;
1639  delete cP;
1640  Th.Clear();
1641 
1644 
1645  if (NURBSext)
1646  {
1647  if (own_ext) { delete NURBSext; }
1648  }
1649  else
1650  {
1651  delete elem_dof;
1652  delete bdrElem_dof;
1653 
1654  delete [] bdofs;
1655  delete [] fdofs;
1656  }
1657 }
1658 
1660  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
1661 {
1662  // Assumptions: see the declaration of the method.
1663 
1664  if (T.Type() == Operator::MFEM_SPARSEMAT)
1665  {
1666  DenseTensor localP;
1667  GetLocalRefinementMatrices(coarse_fes, localP);
1668  T.Reset(RefinementMatrix_main(coarse_fes.GetNDofs(),
1669  coarse_fes.GetElementToDofTable(),
1670  localP));
1671  }
1672  else
1673  {
1674  T.Reset(new RefinementOperator(this, &coarse_fes));
1675  }
1676 }
1677 
1679  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
1680 {
1681  const SparseMatrix *coarse_P = coarse_fes.GetConformingProlongation();
1682 
1683  Operator::Type req_type = T.Type();
1684  GetTransferOperator(coarse_fes, T);
1685 
1686  if (req_type == Operator::MFEM_SPARSEMAT)
1687  {
1689  {
1690  T.Reset(mfem::Mult(*cR, *T.As<SparseMatrix>()));
1691  }
1692  if (coarse_P)
1693  {
1694  T.Reset(mfem::Mult(*T.As<SparseMatrix>(), *coarse_P));
1695  }
1696  }
1697  else
1698  {
1699  const int RP_case = bool(GetConformingRestriction()) + 2*bool(coarse_P);
1700  if (RP_case == 0) { return; }
1701  const bool owner = T.OwnsOperator();
1702  T.SetOperatorOwner(false);
1703  switch (RP_case)
1704  {
1705  case 1:
1706  T.Reset(new ProductOperator(cR, T.Ptr(), false, owner));
1707  break;
1708  case 2:
1709  T.Reset(new ProductOperator(T.Ptr(), coarse_P, owner, false));
1710  break;
1711  case 3:
1713  cR, T.Ptr(), coarse_P, false, owner, false));
1714  break;
1715  }
1716  }
1717 }
1718 
1719 void FiniteElementSpace::Update(bool want_transform)
1720 {
1721  if (mesh->GetSequence() == sequence)
1722  {
1723  return; // mesh and space are in sync, no-op
1724  }
1725  if (want_transform && mesh->GetSequence() != sequence + 1)
1726  {
1727  MFEM_ABORT("Error in update sequence. Space needs to be updated after "
1728  "each mesh modification.");
1729  }
1730  sequence = mesh->GetSequence();
1731 
1732  if (NURBSext)
1733  {
1734  UpdateNURBS();
1735  return;
1736  }
1737 
1738  Table* old_elem_dof = NULL;
1739  int old_ndofs;
1740 
1741  // save old DOF table
1742  if (want_transform)
1743  {
1744  old_elem_dof = elem_dof;
1745  elem_dof = NULL;
1746  old_ndofs = ndofs;
1747  }
1748 
1749  Destroy(); // calls Th.Clear()
1750  Construct();
1752 
1753  if (want_transform)
1754  {
1755  // calculate appropriate GridFunction transformation
1756  switch (mesh->GetLastOperation())
1757  {
1758  case Mesh::REFINE:
1759  {
1760  if (Th.Type() != Operator::MFEM_SPARSEMAT)
1761  {
1762  Th.Reset(new RefinementOperator(this, old_elem_dof, old_ndofs));
1763  // The RefinementOperator takes ownership of 'old_elem_dof', so
1764  // we no longer own it:
1765  old_elem_dof = NULL;
1766  }
1767  else
1768  {
1769  // calculate fully assembled matrix
1770  Th.Reset(RefinementMatrix(old_ndofs, old_elem_dof));
1771  }
1772  break;
1773  }
1774 
1775  case Mesh::DEREFINE:
1776  {
1778  Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof));
1779  if (cP && cR)
1780  {
1781  Th.SetOperatorOwner(false);
1783  false, false, true));
1784  }
1785  break;
1786  }
1787 
1788  default:
1789  break;
1790  }
1791 
1792  delete old_elem_dof;
1793  }
1794 }
1795 
1796 void FiniteElementSpace::Save(std::ostream &out) const
1797 {
1798  int fes_format = 90; // the original format, v0.9
1799  bool nurbs_unit_weights = false;
1800 
1801  // Determine the format that should be used.
1802  if (!NURBSext)
1803  {
1804  // TODO: if this is a variable-order FE space, use fes_format = 100.
1805  }
1806  else
1807  {
1808  const NURBSFECollection *nurbs_fec =
1809  dynamic_cast<const NURBSFECollection *>(fec);
1810  MFEM_VERIFY(nurbs_fec, "invalid FE collection");
1811  nurbs_fec->SetOrder(NURBSext->GetOrder());
1812  const double eps = 5e-14;
1813  nurbs_unit_weights = (NURBSext->GetWeights().Min() >= 1.0-eps &&
1814  NURBSext->GetWeights().Max() <= 1.0+eps);
1816  (NURBSext != mesh->NURBSext && !nurbs_unit_weights))
1817  {
1818  fes_format = 100; // v1.0 format
1819  }
1820  }
1821 
1822  out << (fes_format == 90 ?
1823  "FiniteElementSpace\n" : "MFEM FiniteElementSpace v1.0\n")
1824  << "FiniteElementCollection: " << fec->Name() << '\n'
1825  << "VDim: " << vdim << '\n'
1826  << "Ordering: " << ordering << '\n';
1827 
1828  if (fes_format == 100) // v1.0
1829  {
1830  if (!NURBSext)
1831  {
1832  // TODO: this is a variable-order FE space --> write 'element_orders'.
1833  }
1834  else if (NURBSext != mesh->NURBSext)
1835  {
1837  {
1838  out << "NURBS_order\n" << NURBSext->GetOrder() << '\n';
1839  }
1840  else
1841  {
1842  out << "NURBS_orders\n";
1843  // 1 = do not write the size, just the entries:
1844  NURBSext->GetOrders().Save(out, 1);
1845  }
1846  // If the weights are not unit, write them to the output:
1847  if (!nurbs_unit_weights)
1848  {
1849  out << "NURBS_weights\n";
1850  NURBSext->GetWeights().Print(out, 1);
1851  }
1852  }
1853  out << "End: MFEM FiniteElementSpace v1.0\n";
1854  }
1855 }
1856 
1858 {
1859  string buff;
1860  int fes_format = 0, ord;
1861  FiniteElementCollection *r_fec;
1862 
1863  Destroy();
1864 
1865  input >> std::ws;
1866  getline(input, buff); // 'FiniteElementSpace'
1867  filter_dos(buff);
1868  if (buff == "FiniteElementSpace") { fes_format = 90; /* v0.9 */ }
1869  else if (buff == "MFEM FiniteElementSpace v1.0") { fes_format = 100; }
1870  else { MFEM_ABORT("input stream is not a FiniteElementSpace!"); }
1871  getline(input, buff, ' '); // 'FiniteElementCollection:'
1872  input >> std::ws;
1873  getline(input, buff);
1874  filter_dos(buff);
1875  r_fec = FiniteElementCollection::New(buff.c_str());
1876  getline(input, buff, ' '); // 'VDim:'
1877  input >> vdim;
1878  getline(input, buff, ' '); // 'Ordering:'
1879  input >> ord;
1880 
1881  NURBSFECollection *nurbs_fec = dynamic_cast<NURBSFECollection*>(r_fec);
1882  NURBSExtension *NURBSext = NULL;
1883  if (fes_format == 90) // original format, v0.9
1884  {
1885  if (nurbs_fec)
1886  {
1887  MFEM_VERIFY(m->NURBSext, "NURBS FE collection requires a NURBS mesh!");
1888  const int order = nurbs_fec->GetOrder();
1889  if (order != m->NURBSext->GetOrder() &&
1891  {
1892  NURBSext = new NURBSExtension(m->NURBSext, order);
1893  }
1894  }
1895  }
1896  else if (fes_format == 100) // v1.0
1897  {
1898  while (1)
1899  {
1900  skip_comment_lines(input, '#');
1901  MFEM_VERIFY(input.good(), "error reading FiniteElementSpace v1.0");
1902  getline(input, buff);
1903  filter_dos(buff);
1904  if (buff == "NURBS_order" || buff == "NURBS_orders")
1905  {
1906  MFEM_VERIFY(nurbs_fec,
1907  buff << ": NURBS FE collection is required!");
1908  MFEM_VERIFY(m->NURBSext, buff << ": NURBS mesh is required!");
1909  MFEM_VERIFY(!NURBSext, buff << ": order redefinition!");
1910  if (buff == "NURBS_order")
1911  {
1912  int order;
1913  input >> order;
1914  NURBSext = new NURBSExtension(m->NURBSext, order);
1915  }
1916  else
1917  {
1918  Array<int> orders;
1919  orders.Load(m->NURBSext->GetNKV(), input);
1920  NURBSext = new NURBSExtension(m->NURBSext, orders);
1921  }
1922  }
1923  else if (buff == "NURBS_weights")
1924  {
1925  MFEM_VERIFY(NURBSext, "NURBS_weights: NURBS_orders have to be "
1926  "specified before NURBS_weights!");
1927  NURBSext->GetWeights().Load(input, NURBSext->GetNDof());
1928  }
1929  else if (buff == "element_orders")
1930  {
1931  MFEM_VERIFY(!nurbs_fec, "section element_orders cannot be used "
1932  "with a NURBS FE collection");
1933  MFEM_ABORT("element_orders: not implemented yet!");
1934  }
1935  else if (buff == "End: MFEM FiniteElementSpace v1.0")
1936  {
1937  break;
1938  }
1939  else
1940  {
1941  MFEM_ABORT("unknown section: " << buff);
1942  }
1943  }
1944  }
1945 
1946  Constructor(m, NURBSext, r_fec, vdim, ord);
1947 
1948  return r_fec;
1949 }
1950 
1951 
1953 {
1954  // protected method
1955  int offset = 0;
1956  const int num_elem = mesh->GetNE();
1957  element_offsets = new int[num_elem + 1];
1958  for (int g = 0; g < Geometry::NumGeom; g++)
1959  {
1960  int_rule[g] = NULL;
1961  }
1962  for (int i = 0; i < num_elem; i++)
1963  {
1964  element_offsets[i] = offset;
1965  int geom = mesh->GetElementBaseGeometry(i);
1966  if (int_rule[geom] == NULL)
1967  {
1968  int_rule[geom] = &IntRules.Get(geom, order);
1969  }
1970  offset += int_rule[geom]->GetNPoints();
1971  }
1972  element_offsets[num_elem] = size = offset;
1973 }
1974 
1975 QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
1976  : mesh(mesh_)
1977 {
1978  const char *msg = "invalid input stream";
1979  string ident;
1980 
1981  in >> ident; MFEM_VERIFY(ident == "QuadratureSpace", msg);
1982  in >> ident; MFEM_VERIFY(ident == "Type:", msg);
1983  in >> ident;
1984  if (ident == "default_quadrature")
1985  {
1986  in >> ident; MFEM_VERIFY(ident == "Order:", msg);
1987  in >> order;
1988  }
1989  else
1990  {
1991  MFEM_ABORT("unknown QuadratureSpace type: " << ident);
1992  return;
1993  }
1994 
1995  Construct();
1996 }
1997 
1998 void QuadratureSpace::Save(std::ostream &out) const
1999 {
2000  out << "QuadratureSpace\n"
2001  << "Type: default_quadrature\n"
2002  << "Order: " << order << '\n';
2003 }
2004 
2005 
2006 } // namespace mfem
Abstract class for Finite Elements.
Definition: fe.hpp:140
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:266
void Set(const double *p, const int dim)
Definition: intrules.hpp:32
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1544
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:2035
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:232
void Load(std::istream &in, int fmt=0)
Read an Array from the stream in using format fmt. The format fmt can be:
Definition: array.cpp:94
int * GetJ()
Definition: table.hpp:114
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
Definition: fespace.hpp:84
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
DenseMatrix & GetPointMat()
Read and write access to the underlying point matrix describing the transformation.
Definition: eltrans.hpp:312
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:177
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:3698
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:2915
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:1719
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
virtual const FiniteElement * TraceFiniteElementForGeometry(int GeomType) const
Definition: fe_coll.hpp:54
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:2509
static const int NumGeom
Definition: geom.hpp:34
bool Nonconforming() const
Definition: fespace.hpp:231
int Dimension() const
Definition: mesh.hpp:645
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:305
void InvertLinearTrans(IsoparametricTransformation &trans, const DenseMatrix &invdfdx, const IntegrationPoint &pt, Vector &x)
Definition: fespace.cpp:944
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof)
Definition: fespace.cpp:855
SparseMatrix * DerefinementMatrix(int old_ndofs, const Table *old_elem_dof)
Calculate GridFunction restriction matrix after mesh derefinement.
Definition: fespace.cpp:1011
Ordering::Type ordering
Definition: fespace.hpp:81
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:118
void GetEdgeInteriorVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:207
void GetVertexVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:195
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
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:271
const SparseMatrix * GetConformingRestriction() const
Definition: fespace.cpp:768
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:224
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:169
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:183
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:605
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
void SetIdentityTransformation(int GeomType)
Definition: eltrans.cpp:370
long GetSequence() const
Definition: mesh.hpp:1027
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
void GetFaceEdges(int i, Array< int > &, Array< int > &) const
Definition: mesh.cpp:3752
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:319
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:76
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:699
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:250
void SetSize(int i, int j, int k)
Definition: densemat.hpp:682
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: fespace.cpp:365
int SizeJ() const
Definition: densemat.hpp:679
const Vector & GetWeights() const
Definition: nurbs.hpp:350
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:624
T Sum()
Sum all entries.
Definition: array.cpp:152
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:2171
void BuildConformingInterpolation() const
Calculate the cP and cR matrices for a nonconforming mesh.
Definition: fespace.cpp:576
int GetFaceBaseGeometry(int i) const
Definition: mesh.cpp:3975
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
Definition: fespace.cpp:1857
STL namespace.
virtual void Transform(const IntegrationPoint &, Vector &)
Definition: eltrans.cpp:461
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: fespace.cpp:896
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2258
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:285
const FiniteElementCollection * fec
Associated FE collection (not owned).
Definition: fespace.hpp:73
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
Definition: fespace.hpp:106
void GetVertexDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1519
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: text.hpp:27
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:214
void DeleteAll()
Delete whole array.
Definition: array.hpp:709
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1335
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Definition: ncmesh.hpp:188
Array< int > dof_elem_array
Definition: fespace.hpp:92
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:384
virtual void GetFaceDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1431
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1597
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Definition: fe.cpp:124
void ConvertFromConformingVDofs(const Array< int > &cdofs, Array< int > &dofs)
Definition: fespace.cpp:420
void GetElementInteriorVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:201
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:615
ID for class SparseMatrix.
Definition: operator.hpp:127
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:49
Geometry Geometries
Definition: geom.cpp:760
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
Definition: ncmesh.cpp:2048
General product operator: x -> (A*B)(x) = A(B(x)).
Definition: operator.hpp:329
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:126
int dim
Definition: ex3.cpp:47
void BuildElementToDofTable() const
Definition: fespace.cpp:213
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:81
DenseTensor point_matrices
matrices for IsoparametricTransformation
Definition: ncmesh.hpp:54
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:614
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:317
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:146
const Array< int > & GetOrders() const
Read-only access to the orders of all knot vectors.
Definition: nurbs.hpp:314
std::vector< Master > masters
Definition: ncmesh.hpp:172
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:67
Type
Ordering methods:
Definition: fespace.hpp:30
SparseMatrix * cR
Conforming restriction matrix such that cR.cP=I.
Definition: fespace.hpp:102
SparseMatrix * D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Definition: fespace.cpp:461
SparseMatrix * D2C_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Definition: fespace.cpp:429
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:3782
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...
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition: mesh.hpp:713
void Get(double *p, const int dim) const
Definition: intrules.hpp:45
NURBSExtension * StealNURBSext()
Definition: fespace.cpp:1138
virtual const char * Name() const
Definition: fe_coll.hpp:50
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:146
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:695
const SparseMatrix * GetConformingProlongation() const
Definition: fespace.cpp:761
Array< Embedding > embeddings
fine element positions in their parents
Definition: ncmesh.hpp:55
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:3952
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I)
Definition: fespace.cpp:526
Array< int > dof_ldof_array
Definition: fespace.hpp:92
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:3720
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:117
SparseMatrix * cP
Definition: fespace.hpp:100
void SetRow(int r, const Vector &row)
Definition: densemat.cpp:2817
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Definition: fespace.cpp:1998
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:189
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:168
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:277
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1556
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:3100
const FiniteElement * GetEdgeElement(int i) const
Definition: fespace.cpp:1620
int slaves_end
slave faces
Definition: ncmesh.hpp:148
void Threshold(double eps)
Replace small entries, abs(a_ij) <= eps, with zero.
Definition: densemat.cpp:2833
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3150
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:242
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:63
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const DenseTensor &localP) const
Definition: fespace.cpp:781
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:124
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:224
RefinementOperator(const FiniteElementSpace *fespace, Table *old_elem_dof, int old_ndofs)
Definition: fespace.cpp:867
int GetBdrElementBaseGeometry(int i=0) const
Definition: mesh.hpp:691
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:2530
std::vector< Slave > slaves
Definition: ncmesh.hpp:173
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:171
Nonconforming edge/face within a bigger edge/face.
Definition: ncmesh.hpp:155
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:66
virtual ~FiniteElementSpace()
Definition: fespace.cpp:1631
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:2296
SparseMatrix * H2L_GlobalRestrictionMatrix(FiniteElementSpace *lfes)
Definition: fespace.cpp:492
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:785
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Definition: fespace.cpp:296
int GetOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:94
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:110
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:217
Table * GetElementDofTable()
Definition: nurbs.hpp:340
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
int GetElementBaseGeometry(int i=0) const
Definition: mesh.hpp:688
Table * GetBdrElementDofTable()
Definition: nurbs.hpp:341
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:772
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:26
void GetLocalRefinementMatrices(DenseTensor &localP) const
Definition: fespace.cpp:832
const CoarseFineTransformations & GetRefinementTransforms()
Definition: mesh.cpp:6853
bool OwnsOperator() const
Return true if the OperatorHandle owns the held Operator.
Definition: handle.hpp:103
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:644
virtual int * DofOrderForOrientation(int GeomType, int Or) const =0
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:1233
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:44
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:1678
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:618
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:176
void ConvertToConformingVDofs(const Array< int > &dofs, Array< int > &cdofs)
Definition: fespace.cpp:412
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:85
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:106
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:1348
int GetOrder() const
Get the order of the NURBS collection: either a positive number, when using fixed order...
Definition: fe_coll.hpp:297
Mesh * mesh
The mesh that FE space lives on (not owned).
Definition: fespace.hpp:70
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:90
void ReorderElementToDofTable()
Reorder the scalar DOFs based on the element ordering.
Definition: fespace.cpp:242
GridFunction interpolation operator applicable after mesh refinement.
Definition: fespace.hpp:136
int GetBdrAttribute(int i) const
Definition: fespace.hpp:315
void Save(std::ostream &out) const
Definition: fespace.cpp:1796
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition: operator.hpp:373
void GetEdgeDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1490
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:298
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th edge.
Definition: fespace.cpp:189
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:286
void MakeVDimMatrix(SparseMatrix &mat) const
Definition: fespace.cpp:733
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:190
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:253
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:3914
int SizeI() const
Definition: densemat.hpp:678
void GetLocalDerefinementMatrices(DenseTensor &localR) const
Definition: fespace.cpp:961
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:42
int GetNKV() const
Definition: nurbs.hpp:319
int parent
element index in the coarse mesh
Definition: ncmesh.hpp:45
FiniteElementSpace()
Default constructor: the object is invalid until initialized using the method Load().
Definition: fespace.cpp:58
QuadratureSpace(Mesh *mesh_, int order_)
Create a QuadratureSpace based on the global rules from IntRules.
Definition: fespace.hpp:551
void Constructor(Mesh *mesh, NURBSExtension *ext, const FiniteElementCollection *fec, int vdim=1, int ordering=Ordering::byNODES)
Help function for constructors + Load().
Definition: fespace.cpp:1092
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
int Size_of_connections() const
Definition: table.hpp:98
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:177
int Size() const
Logical size of the array.
Definition: array.hpp:133
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2210
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Definition: fespace.cpp:401
void filter_dos(std::string &line)
Definition: text.hpp:41
void GetEntityDofs(int entity, int index, Array< int > &dofs) const
Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
Definition: fespace.cpp:565
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:720
int GetFaceOrder(int i) const
Returns the order of the i&#39;th face finite element.
Definition: fespace.cpp:100
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:142
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:627
Vector data type.
Definition: vector.hpp:48
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:531
int SizeK() const
Definition: densemat.hpp:680
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:297
const FiniteElement * GetTraceElement(int i, int geom_type) const
Return the trace element from element &#39;i&#39; to the given &#39;geom_type&#39;.
Definition: fespace.cpp:1625
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:52
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition: ncmesh.hpp:195
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1531
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Definition: fespace.cpp:550
NURBSExtension * NURBSext
Definition: fespace.hpp:94
virtual int DofForGeometry(int GeomType) const =0
Operation GetLastOperation() const
Return type of last modification of the mesh.
Definition: mesh.hpp:1021
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:64
Abstract operator.
Definition: operator.hpp:21
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:3260
const IntegrationRule & GetNodes() const
Definition: fe.hpp:270
int GetNConformingDofs() const
Definition: fespace.cpp:775
int index
Mesh number.
Definition: ncmesh.hpp:136
BiLinear2DFiniteElement QuadrilateralFE
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:638
int HasFaceDofs(int GeomType) const
Definition: fe_coll.cpp:25
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:2439
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:353
const Table & GetElementToDofTable() const
Definition: fespace.hpp:384
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i&#39;th boundary element.
Definition: fespace.cpp:1571
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
bool Conforming() const
Definition: fespace.hpp:230
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:106
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:221
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:158
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:118
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:1659
Defines the position of a fine element within a coarse element.
Definition: ncmesh.hpp:43
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:131
int GetNDof() const
Definition: nurbs.hpp:329
int matrix
index into CoarseFineTransformations::point_matrices
Definition: ncmesh.hpp:46
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const =0
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, but treat all elements as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:667