MFEM  v3.3.2
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 "../mesh/mesh_headers.hpp"
15 #include "fem.hpp"
16 
17 #include <cmath>
18 #include <cstdarg>
19 #include <limits>
20 
21 using namespace std;
22 
23 namespace mfem
24 {
25 
26 template <> void Ordering::
27 DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
28 {
29  // static method
30  int size = dofs.Size();
31  dofs.SetSize(size*vdim);
32  for (int vd = 1; vd < vdim; vd++)
33  {
34  for (int i = 0; i < size; i++)
35  {
36  dofs[i+size*vd] = Map<byNODES>(ndofs, vdim, dofs[i], vd);
37  }
38  }
39 }
40 
41 template <> void Ordering::
42 DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
43 {
44  // static method
45  int size = dofs.Size();
46  dofs.SetSize(size*vdim);
47  for (int vd = vdim-1; vd >= 0; vd--)
48  {
49  for (int i = 0; i < size; i++)
50  {
51  dofs[i+size*vd] = Map<byVDIM>(ndofs, vdim, dofs[i], vd);
52  }
53  }
54 }
55 
56 int FiniteElementSpace::GetOrder(int i) const
57 {
58  int GeomType = mesh->GetElementBaseGeometry(i);
59  return fec->FiniteElementForGeometry(GeomType)->GetOrder();
60 }
61 
62 int FiniteElementSpace::GetFaceOrder(int i) const
63 {
64  int GeomType = mesh->GetFaceBaseGeometry(i);
65  return fec->FiniteElementForGeometry(GeomType)->GetOrder();
66 }
67 
68 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs, int ndofs) const
69 {
70  if (vdim == 1) { return; }
71  if (ndofs < 0) { ndofs = this->ndofs; }
72 
73  if (ordering == Ordering::byNODES)
74  {
75  Ordering::DofsToVDofs<Ordering::byNODES>(ndofs, vdim, dofs);
76  }
77  else
78  {
79  Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs, vdim, dofs);
80  }
81 }
82 
83 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs, int ndofs) const
84 {
85  if (vdim == 1) { return; }
86  if (ndofs < 0) { ndofs = this->ndofs; }
87 
88  if (ordering == Ordering::byNODES)
89  {
90  for (int i = 0; i < dofs.Size(); i++)
91  {
92  dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs, vdim, dofs[i], vd);
93  }
94  }
95  else
96  {
97  for (int i = 0; i < dofs.Size(); i++)
98  {
99  dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dofs[i], vd);
100  }
101  }
102 }
103 
104 int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs) const
105 {
106  if (vdim == 1) { return dof; }
107  if (ndofs < 0) { ndofs = this->ndofs; }
108 
109  if (ordering == Ordering::byNODES)
110  {
111  return Ordering::Map<Ordering::byNODES>(ndofs, vdim, dof, vd);
112  }
113  else
114  {
115  return Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dof, vd);
116  }
117 }
118 
119 // static function
120 void FiniteElementSpace::AdjustVDofs (Array<int> &vdofs)
121 {
122  int n = vdofs.Size(), *vdof = vdofs;
123  for (int i = 0; i < n; i++)
124  {
125  int j;
126  if ((j = vdof[i]) < 0)
127  {
128  vdof[i] = -1-j;
129  }
130  }
131 }
132 
133 void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs) const
134 {
135  GetElementDofs(i, vdofs);
136  DofsToVDofs(vdofs);
137 }
138 
139 void FiniteElementSpace::GetBdrElementVDofs(int i, Array<int> &vdofs) const
140 {
141  GetBdrElementDofs(i, vdofs);
142  DofsToVDofs(vdofs);
143 }
144 
145 void FiniteElementSpace::GetFaceVDofs(int i, Array<int> &vdofs) const
146 {
147  GetFaceDofs(i, vdofs);
148  DofsToVDofs(vdofs);
149 }
150 
151 void FiniteElementSpace::GetEdgeVDofs(int i, Array<int> &vdofs) const
152 {
153  GetEdgeDofs(i, vdofs);
154  DofsToVDofs(vdofs);
155 }
156 
157 void FiniteElementSpace::GetVertexVDofs(int i, Array<int> &vdofs) const
158 {
159  GetVertexDofs(i, vdofs);
160  DofsToVDofs(vdofs);
161 }
162 
163 void FiniteElementSpace::GetElementInteriorVDofs(int i, Array<int> &vdofs) const
164 {
165  GetElementInteriorDofs(i, vdofs);
166  DofsToVDofs(vdofs);
167 }
168 
169 void FiniteElementSpace::GetEdgeInteriorVDofs(int i, Array<int> &vdofs) const
170 {
171  GetEdgeInteriorDofs(i, vdofs);
172  DofsToVDofs(vdofs);
173 }
174 
175 void FiniteElementSpace::BuildElementToDofTable() const
176 {
177  if (elem_dof) { return; }
178 
179  Table *el_dof = new Table;
180  Array<int> dofs;
181  el_dof -> MakeI (mesh -> GetNE());
182  for (int i = 0; i < mesh -> GetNE(); i++)
183  {
184  GetElementDofs (i, dofs);
185  el_dof -> AddColumnsInRow (i, dofs.Size());
186  }
187  el_dof -> MakeJ();
188  for (int i = 0; i < mesh -> GetNE(); i++)
189  {
190  GetElementDofs (i, dofs);
191  el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
192  }
193  el_dof -> ShiftUpI();
194  elem_dof = el_dof;
195 }
196 
197 void FiniteElementSpace::RebuildElementToDofTable()
198 {
199  delete elem_dof;
200  elem_dof = NULL;
201  BuildElementToDofTable();
202 }
203 
204 void FiniteElementSpace::ReorderElementToDofTable()
205 {
206  Array<int> dof_marker(ndofs);
207 
208  dof_marker = -1;
209 
210  int *J = elem_dof->GetJ(), nnz = elem_dof->Size_of_connections();
211  for (int k = 0, dof_counter = 0; k < nnz; k++)
212  {
213  const int sdof = J[k]; // signed dof
214  const int dof = (sdof < 0) ? -1-sdof : sdof;
215  int new_dof = dof_marker[dof];
216  if (new_dof < 0)
217  {
218  dof_marker[dof] = new_dof = dof_counter++;
219  }
220  J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
221  }
222 }
223 
224 void FiniteElementSpace::BuildDofToArrays()
225 {
226  if (dof_elem_array.Size()) { return; }
227 
228  BuildElementToDofTable();
229 
230  dof_elem_array.SetSize (ndofs);
231  dof_ldof_array.SetSize (ndofs);
232  dof_elem_array = -1;
233  for (int i = 0; i < mesh -> GetNE(); i++)
234  {
235  const int *dofs = elem_dof -> GetRow(i);
236  const int n = elem_dof -> RowSize(i);
237  for (int j = 0; j < n; j++)
238  {
239  if (dof_elem_array[dofs[j]] < 0)
240  {
241  dof_elem_array[dofs[j]] = i;
242  dof_ldof_array[dofs[j]] = j;
243  }
244  }
245  }
246 }
247 
248 static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
249 {
250  for (int i = 0; i < dofs.Size(); i++)
251  {
252  int k = dofs[i];
253  if (k < 0) { k = -1 - k; }
254  mark_array[k] = -1;
255  }
256 }
257 
258 void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
259  Array<int> &ess_vdofs,
260  int component) const
261 {
262  Array<int> vdofs, dofs;
263 
264  ess_vdofs.SetSize(GetVSize());
265  ess_vdofs = 0;
266 
267  for (int i = 0; i < GetNBE(); i++)
268  {
269  if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
270  {
271  if (component < 0)
272  {
273  // Mark all components.
274  GetBdrElementVDofs(i, vdofs);
275  mark_dofs(vdofs, ess_vdofs);
276  }
277  else
278  {
279  GetBdrElementDofs(i, dofs);
280  for (int d = 0; d < dofs.Size(); d++)
281  { dofs[d] = DofToVDof(dofs[d], component); }
282  mark_dofs(dofs, ess_vdofs);
283  }
284  }
285  }
286 
287  // mark possible hidden boundary edges in a non-conforming mesh, also
288  // local DOFs affected by boundary elements on other processors
289  if (mesh->ncmesh)
290  {
291  Array<int> bdr_verts, bdr_edges;
292  mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges);
293 
294  for (int i = 0; i < bdr_verts.Size(); i++)
295  {
296  if (component < 0)
297  {
298  GetVertexVDofs(bdr_verts[i], vdofs);
299  mark_dofs(vdofs, ess_vdofs);
300  }
301  else
302  {
303  GetVertexDofs(bdr_verts[i], dofs);
304  for (int d = 0; d < dofs.Size(); d++)
305  { dofs[d] = DofToVDof(dofs[d], component); }
306  mark_dofs(dofs, ess_vdofs);
307  }
308  }
309  for (int i = 0; i < bdr_edges.Size(); i++)
310  {
311  if (component < 0)
312  {
313  GetEdgeVDofs(bdr_edges[i], vdofs);
314  mark_dofs(vdofs, ess_vdofs);
315  }
316  else
317  {
318  GetEdgeDofs(bdr_edges[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 
327 void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
328  Array<int> &ess_tdof_list,
329  int component)
330 {
331  Array<int> ess_vdofs, ess_tdofs;
332  GetEssentialVDofs(bdr_attr_is_ess, ess_vdofs, component);
333  const SparseMatrix *R = GetConformingRestriction();
334  if (!R)
335  {
336  ess_tdofs.MakeRef(ess_vdofs);
337  }
338  else
339  {
340  R->BooleanMult(ess_vdofs, ess_tdofs);
341  }
342  MarkerToList(ess_tdofs, ess_tdof_list);
343 }
344 
345 // static method
346 void FiniteElementSpace::MarkerToList(const Array<int> &marker,
347  Array<int> &list)
348 {
349  int num_marked = 0;
350  for (int i = 0; i < marker.Size(); i++)
351  {
352  if (marker[i]) { num_marked++; }
353  }
354  list.SetSize(0);
355  list.Reserve(num_marked);
356  for (int i = 0; i < marker.Size(); i++)
357  {
358  if (marker[i]) { list.Append(i); }
359  }
360 }
361 
362 // static method
363 void FiniteElementSpace::ListToMarker(const Array<int> &list, int marker_size,
364  Array<int> &marker, int mark_val)
365 {
366  marker.SetSize(marker_size);
367  marker = 0;
368  for (int i = 0; i < list.Size(); i++)
369  {
370  marker[list[i]] = mark_val;
371  }
372 }
373 
374 void FiniteElementSpace::ConvertToConformingVDofs(const Array<int> &dofs,
375  Array<int> &cdofs)
376 {
377  GetConformingProlongation();
378  if (cP) { cP->BooleanMultTranspose(dofs, cdofs); }
379  else { dofs.Copy(cdofs); }
380 }
381 
382 void FiniteElementSpace::ConvertFromConformingVDofs(const Array<int> &cdofs,
383  Array<int> &dofs)
384 {
385  GetConformingRestriction();
386  if (cR) { cR->BooleanMultTranspose(cdofs, dofs); }
387  else { cdofs.Copy(dofs); }
388 }
389 
390 SparseMatrix *
391 FiniteElementSpace::D2C_GlobalRestrictionMatrix (FiniteElementSpace *cfes)
392 {
393  int i, j;
394  Array<int> d_vdofs, c_vdofs;
395  SparseMatrix *R;
396 
397  R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
398 
399  for (i = 0; i < mesh -> GetNE(); i++)
400  {
401  this -> GetElementVDofs (i, d_vdofs);
402  cfes -> GetElementVDofs (i, c_vdofs);
403 
404 #ifdef MFEM_DEBUG
405  if (d_vdofs.Size() != c_vdofs.Size())
406  {
407  mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
408  }
409 #endif
410 
411  for (j = 0; j < d_vdofs.Size(); j++)
412  {
413  R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
414  }
415  }
416 
417  R -> Finalize();
418 
419  return R;
420 }
421 
422 SparseMatrix *
423 FiniteElementSpace::D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
424 {
425  int i, j;
426  Array<int> d_dofs, c_dofs;
427  SparseMatrix *R;
428 
429  R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
430 
431  for (i = 0; i < mesh -> GetNE(); i++)
432  {
433  this -> GetElementDofs (i, d_dofs);
434  cfes -> GetElementDofs (i, c_dofs);
435 
436 #ifdef MFEM_DEBUG
437  if (c_dofs.Size() != 1)
438  mfem_error ("FiniteElementSpace::"
439  "D2Const_GlobalRestrictionMatrix (...)");
440 #endif
441 
442  for (j = 0; j < d_dofs.Size(); j++)
443  {
444  R -> Set (c_dofs[0], d_dofs[j], 1.0);
445  }
446  }
447 
448  R -> Finalize();
449 
450  return R;
451 }
452 
453 SparseMatrix *
454 FiniteElementSpace::H2L_GlobalRestrictionMatrix (FiniteElementSpace *lfes)
455 {
456  SparseMatrix *R;
457  DenseMatrix loc_restr;
458  Array<int> l_dofs, h_dofs;
459 
460  R = new SparseMatrix (lfes -> GetNDofs(), ndofs);
461 
462  if (!lfes->GetNE())
463  {
464  R->Finalize();
465  return R;
466  }
467 
468  const FiniteElement *h_fe = this -> GetFE (0);
469  const FiniteElement *l_fe = lfes -> GetFE (0);
472  h_fe->Project(*l_fe, T, loc_restr);
473 
474  for (int i = 0; i < mesh -> GetNE(); i++)
475  {
476  this -> GetElementDofs (i, h_dofs);
477  lfes -> GetElementDofs (i, l_dofs);
478 
479  R -> SetSubMatrix (l_dofs, h_dofs, loc_restr, 1);
480  }
481 
482  R -> Finalize();
483 
484  return R;
485 }
486 
487 static void AddDependencies(SparseMatrix& deps, Array<int>& master_dofs,
488  Array<int>& slave_dofs, DenseMatrix& I)
489 {
490  for (int i = 0; i < slave_dofs.Size(); i++)
491  {
492  int sdof = slave_dofs[i];
493  if (!deps.RowSize(sdof)) // not processed yet?
494  {
495  for (int j = 0; j < master_dofs.Size(); j++)
496  {
497  double coef = I(i, j);
498  if (std::abs(coef) > 1e-12)
499  {
500  int mdof = master_dofs[j];
501  if (mdof != sdof && mdof != (-1-sdof))
502  {
503  deps.Add(sdof, mdof, coef);
504  }
505  }
506  }
507  }
508  }
509 }
510 
511 static bool DofFinalizable(int dof, const Array<bool>& finalized,
512  const SparseMatrix& deps)
513 {
514  const int* dep = deps.GetRowColumns(dof);
515  int ndep = deps.RowSize(dof);
516 
517  // are all constraining DOFs finalized?
518  for (int i = 0; i < ndep; i++)
519  {
520  if (!finalized[dep[i]]) { return false; }
521  }
522  return true;
523 }
524 
525 /** This is a helper function to get edge (type == 0) or face (type == 1) DOFs.
526  The function is aware of ghost edges/faces in parallel, for which an empty
527  DOF list is returned. */
528 void FiniteElementSpace::GetEdgeFaceDofs(int type, int index, Array<int> &dofs)
529 const
530 {
531  dofs.SetSize(0);
532  if (type)
533  {
534  if (index < mesh->GetNFaces()) { GetFaceDofs(index, dofs); }
535  }
536  else
537  {
538  if (index < mesh->GetNEdges()) { GetEdgeDofs(index, dofs); }
539  }
540 }
541 
542 void FiniteElementSpace::GetConformingInterpolation() const
543 {
544 #ifdef MFEM_USE_MPI
545  MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
546  "This method should not be used with a ParFiniteElementSpace!");
547 #endif
548  if (cP_is_set) { return; }
549  cP_is_set = true;
550 
551  // For each slave DOF, the dependency matrix will contain a row that
552  // expresses the slave DOF as a linear combination of its immediate master
553  // DOFs. Rows of independent DOFs will remain empty.
554  SparseMatrix deps(ndofs);
555 
556  // collect local edge/face dependencies
557  for (int type = 0; type <= 1; type++)
558  {
559  const NCMesh::NCList &list = type ? mesh->ncmesh->GetFaceList()
560  /**/ : mesh->ncmesh->GetEdgeList();
561  if (!list.masters.size()) { continue; }
562 
564  if (type) { T.SetFE(&QuadrilateralFE); }
565  else { T.SetFE(&SegmentFE); }
566 
567  int geom = type ? Geometry::SQUARE : Geometry::SEGMENT;
568  const FiniteElement* fe = fec->FiniteElementForGeometry(geom);
569  if (!fe) { continue; }
570 
571  Array<int> master_dofs, slave_dofs;
572  DenseMatrix I(fe->GetDof());
573 
574  // loop through all master edges/faces, constrain their slave edges/faces
575  for (unsigned mi = 0; mi < list.masters.size(); mi++)
576  {
577  const NCMesh::Master &master = list.masters[mi];
578  GetEdgeFaceDofs(type, master.index, master_dofs);
579  if (!master_dofs.Size()) { continue; }
580 
581  for (int si = master.slaves_begin; si < master.slaves_end; si++)
582  {
583  const NCMesh::Slave &slave = list.slaves[si];
584  GetEdgeFaceDofs(type, slave.index, slave_dofs);
585  if (!slave_dofs.Size()) { continue; }
586 
587  slave.OrientedPointMatrix(T.GetPointMat());
589  fe->GetLocalInterpolation(T, I);
590 
591  // make each slave DOF dependent on all master DOFs
592  AddDependencies(deps, master_dofs, slave_dofs, I);
593  }
594  }
595  }
596 
597  deps.Finalize();
598 
599  // DOFs that stayed independent are true DOFs
600  int n_true_dofs = 0;
601  for (int i = 0; i < ndofs; i++)
602  {
603  if (!deps.RowSize(i)) { n_true_dofs++; }
604  }
605 
606  // if all dofs are true dofs leave cP and cR NULL
607  if (n_true_dofs == ndofs)
608  {
609  cP = cR = NULL; // will be treated as identities
610  return;
611  }
612 
613  // create the conforming restriction matrix cR
614  int *cR_J;
615  {
616  int *cR_I = new int[n_true_dofs+1];
617  double *cR_A = new double[n_true_dofs];
618  cR_J = new int[n_true_dofs];
619  for (int i = 0; i < n_true_dofs; i++)
620  {
621  cR_I[i] = i;
622  cR_A[i] = 1.0;
623  }
624  cR_I[n_true_dofs] = n_true_dofs;
625  cR = new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, ndofs);
626  }
627 
628  // create the conforming prolongation matrix cP
629  cP = new SparseMatrix(ndofs, n_true_dofs);
630 
631  Array<bool> finalized(ndofs);
632  finalized = false;
633 
634  // put identity in the restriction and prolongation matrices for true DOFs
635  for (int i = 0, true_dof = 0; i < ndofs; i++)
636  {
637  if (!deps.RowSize(i))
638  {
639  cR_J[true_dof] = i;
640  cP->Add(i, true_dof++, 1.0);
641  finalized[i] = true;
642  }
643  }
644 
645  // Now calculate cP rows of slave DOFs as combinations of cP rows of their
646  // master DOFs. It is possible that some slave DOFs depend on DOFs that are
647  // themselves slaves. Here we resolve such indirect constraints by first
648  // calculating rows of the cP matrix for DOFs whose master DOF cP rows are
649  // already known (in the first iteration these are the true DOFs). In the
650  // second iteration, slaves of slaves can be 'finalized' (given a row in the
651  // cP matrix), in the third iteration slaves of slaves of slaves, etc.
652  bool finished;
653  int n_finalized = n_true_dofs;
654  Array<int> cols;
655  Vector srow;
656  do
657  {
658  finished = true;
659  for (int dof = 0; dof < ndofs; dof++)
660  {
661  if (!finalized[dof] && DofFinalizable(dof, finalized, deps))
662  {
663  const int* dep_col = deps.GetRowColumns(dof);
664  const double* dep_coef = deps.GetRowEntries(dof);
665  int n_dep = deps.RowSize(dof);
666 
667  for (int j = 0; j < n_dep; j++)
668  {
669  cP->GetRow(dep_col[j], cols, srow);
670  srow *= dep_coef[j];
671  cP->AddRow(dof, cols, srow);
672  }
673 
674  finalized[dof] = true;
675  n_finalized++;
676  finished = false;
677  }
678  }
679  }
680  while (!finished);
681 
682  // if everything is consistent (mesh, face orientations, etc.), we should
683  // be able to finalize all slave DOFs, otherwise it's a serious error
684  if (n_finalized != ndofs)
685  {
686  MFEM_ABORT("Error creating cP matrix.");
687  }
688 
689  cP->Finalize();
690 
691  if (vdim > 1)
692  {
693  MakeVDimMatrix(*cP);
694  MakeVDimMatrix(*cR);
695  }
696 }
697 
698 void FiniteElementSpace::MakeVDimMatrix(SparseMatrix &mat) const
699 {
700  if (vdim == 1) { return; }
701 
702  int height = mat.Height();
703  int width = mat.Width();
704 
705  SparseMatrix *vmat = new SparseMatrix(vdim*height, vdim*width);
706 
707  Array<int> dofs, vdofs;
708  Vector srow;
709  for (int i = 0; i < height; i++)
710  {
711  mat.GetRow(i, dofs, srow);
712  for (int vd = 0; vd < vdim; vd++)
713  {
714  dofs.Copy(vdofs);
715  DofsToVDofs(vd, vdofs, width);
716  vmat->SetRow(DofToVDof(i, vd, height), vdofs, srow);
717  }
718  }
719  vmat->Finalize();
720 
721  mat.Swap(*vmat);
722  delete vmat;
723 }
724 
725 const SparseMatrix* FiniteElementSpace::GetConformingProlongation() const
726 {
727  if (Conforming()) { return NULL; }
728  if (!cP_is_set) { GetConformingInterpolation(); }
729  return cP;
730 }
731 
732 const SparseMatrix* FiniteElementSpace::GetConformingRestriction() const
733 {
734  if (Conforming()) { return NULL; }
735  if (!cP_is_set) { GetConformingInterpolation(); }
736  return cR;
737 }
738 
739 int FiniteElementSpace::GetNConformingDofs() const
740 {
741  const SparseMatrix* P = GetConformingProlongation();
742  return P ? (P->Width() / vdim) : ndofs;
743 }
744 
745 SparseMatrix* FiniteElementSpace::RefinementMatrix(int old_ndofs,
746  const Table* old_elem_dof)
747 {
748  MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
749  MFEM_VERIFY(ndofs >= old_ndofs, "Previous space is not coarser.");
750 
751  Array<int> dofs, old_dofs, old_vdofs;
752  Vector row;
753 
754  const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
755 
756  int geom = mesh->GetElementBaseGeometry(); // assuming the same geom
757  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
758 
760  isotr.SetIdentityTransformation(geom);
761 
762  int nmat = rtrans.point_matrices.SizeK();
763  int ldof = fe->GetDof(); // assuming the same FE everywhere
764 
765  // calculate local interpolation matrices for all refinement types
766  DenseTensor localP(ldof, ldof, nmat);
767  for (int i = 0; i < nmat; i++)
768  {
769  isotr.GetPointMat() = rtrans.point_matrices(i);
770  isotr.FinalizeTransformation();
771  fe->GetLocalInterpolation(isotr, localP(i));
772  }
773 
774  SparseMatrix *P = new SparseMatrix(ndofs*vdim, old_ndofs*vdim, ldof);
775 
776  Array<int> mark(P->Height());
777  mark = 0;
778  for (int k = 0; k < mesh->GetNE(); k++)
779  {
780  const Embedding &emb = rtrans.embeddings[k];
781  DenseMatrix &lP = localP(emb.matrix);
782 
783  elem_dof->GetRow(k, dofs);
784  old_elem_dof->GetRow(emb.parent, old_dofs);
785 
786  for (int vd = 0; vd < vdim; vd++)
787  {
788  old_dofs.Copy(old_vdofs);
789  DofsToVDofs(vd, old_vdofs, old_ndofs);
790 
791  for (int i = 0; i < ldof; i++)
792  {
793  int r = DofToVDof(dofs[i], vd);
794  int m = (r >= 0) ? r : (-1 - r);
795 
796  if (!mark[m])
797  {
798  lP.GetRow(i, row);
799  P->SetRow(r, old_vdofs, row); // TODO: profiler says this is slow!
800  mark[m] = 1;
801  }
802  }
803  }
804  }
805 
806  MFEM_ASSERT(mark.Sum() == P->Height(), "Not all rows of P set.");
807  return P;
808 }
809 
811  const DenseMatrix &invdfdx,
812  const IntegrationPoint &pt, Vector &x)
813 {
814  // invert a linear transform with one Newton step
815  IntegrationPoint p0;
816  p0.Set3(0, 0, 0);
817  trans.Transform(p0, x);
818 
819  double store[3];
820  Vector v(store, x.Size());
821  pt.Get(v, x.Size());
822  v -= x;
823 
824  invdfdx.Mult(v, x);
825 }
826 
827 void FiniteElementSpace::GetLocalDerefinementMatrices(
828  int geom, const CoarseFineTransformations &dt, DenseTensor &localR)
829 {
830  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
831  const IntegrationRule &nodes = fe->GetNodes();
832 
833  LinearFECollection linfec;
835  isotr.SetFE(linfec.FiniteElementForGeometry(geom));
836 
837  int nmat = dt.point_matrices.SizeK();
838  int ldof = fe->GetDof();
839  int dim = mesh->Dimension();
840 
841  DenseMatrix invdfdx(dim);
842  IntegrationPoint ipt;
843  Vector pt(&ipt.x, dim), shape(ldof);
844 
845  // calculate local restriction matrices for all refinement types
846  localR.SetSize(ldof, ldof, nmat);
847  for (int i = 0; i < nmat; i++)
848  {
849  DenseMatrix &lR = localR(i);
850  lR = numeric_limits<double>::infinity(); // marks invalid rows
851 
852  isotr.GetPointMat() = dt.point_matrices(i);
853  isotr.FinalizeTransformation();
854  isotr.SetIntPoint(&nodes[0]);
855  CalcInverse(isotr.Jacobian(), invdfdx);
856 
857  for (int j = 0; j < nodes.Size(); j++)
858  {
859  InvertLinearTrans(isotr, invdfdx, nodes[j], pt);
860  if (Geometries.CheckPoint(geom, ipt)) // do we need an epsilon here?
861  {
862  IntegrationPoint ip;
863  ip.Set(pt, dim);
864  fe->CalcShape(ip, shape); // TODO: H(curl), etc.?
865  MFEM_ASSERT(dynamic_cast<const NodalFiniteElement*>(fe),
866  "only nodal FEs are implemented");
867  lR.SetRow(j, shape);
868  }
869  }
870  }
871 }
872 
873 SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
874  const Table* old_elem_dof)
875 {
876  MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
877  MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
878  MFEM_VERIFY(ndofs <= old_ndofs, "Previous space is not finer.");
879 
880  Array<int> dofs, old_dofs, old_vdofs;
881  Vector row;
882 
883  const CoarseFineTransformations &dtrans =
884  mesh->ncmesh->GetDerefinementTransforms();
885 
886  int geom = mesh->ncmesh->GetElementGeometry();
887  int ldof = fec->FiniteElementForGeometry(geom)->GetDof();
888 
889  DenseTensor localR;
890  GetLocalDerefinementMatrices(geom, dtrans, localR);
891 
892  SparseMatrix *R = new SparseMatrix(ndofs*vdim, old_ndofs*vdim, ldof);
893 
894  Array<int> mark(R->Height());
895  mark = 0;
896  for (int k = 0; k < dtrans.embeddings.Size(); k++)
897  {
898  const Embedding &emb = dtrans.embeddings[k];
899  DenseMatrix &lR = localR(emb.matrix);
900 
901  elem_dof->GetRow(emb.parent, dofs);
902  old_elem_dof->GetRow(k, old_dofs);
903 
904  for (int vd = 0; vd < vdim; vd++)
905  {
906  old_dofs.Copy(old_vdofs);
907  DofsToVDofs(vd, old_vdofs, old_ndofs);
908 
909  for (int i = 0; i < lR.Height(); i++)
910  {
911  if (lR(i, 0) == numeric_limits<double>::infinity()) { continue; }
912 
913  int r = DofToVDof(dofs[i], vd);
914  int m = (r >= 0) ? r : (-1 - r);
915 
916  if (!mark[m])
917  {
918  lR.GetRow(i, row);
919  R->SetRow(r, old_vdofs, row);
920  mark[m] = 1;
921  }
922  }
923  }
924  }
925 
926  MFEM_ASSERT(mark.Sum() == R->Height(), "Not all rows of R set.");
927  return R;
928 }
929 
930 FiniteElementSpace::FiniteElementSpace(Mesh *mesh,
931  const FiniteElementCollection *fec,
932  int vdim, int ordering)
933 {
934  this->mesh = mesh;
935  this->fec = fec;
936  this->vdim = vdim;
937  this->ordering = (Ordering::Type) ordering;
938 
939  elem_dof = NULL;
940  sequence = mesh->GetSequence();
941 
942  const NURBSFECollection *nurbs_fec =
943  dynamic_cast<const NURBSFECollection *>(fec);
944  if (nurbs_fec)
945  {
946  if (!mesh->NURBSext)
947  {
948  mfem_error("FiniteElementSpace::FiniteElementSpace :\n"
949  " NURBS FE space requires NURBS mesh.");
950  }
951  else
952  {
953  int Order = nurbs_fec->GetOrder();
954  if (mesh->NURBSext->GetOrder() == Order)
955  {
956  NURBSext = mesh->NURBSext;
957  own_ext = 0;
958  }
959  else
960  {
961  NURBSext = new NURBSExtension(mesh->NURBSext, Order);
962  own_ext = 1;
963  }
964  UpdateNURBS();
965  cP = cR = NULL;
966  cP_is_set = false;
967  T = NULL;
968  own_T = true;
969  }
970  }
971  else
972  {
973  NURBSext = NULL;
974  own_ext = 0;
975  Construct();
976  }
977 
978  BuildElementToDofTable();
979 }
980 
981 NURBSExtension *FiniteElementSpace::StealNURBSext()
982 {
983  if (NURBSext && !own_ext)
984  {
985  mfem_error("FiniteElementSpace::StealNURBSext");
986  }
987  own_ext = 0;
988 
989  return NURBSext;
990 }
991 
992 void FiniteElementSpace::UpdateNURBS()
993 {
994  nvdofs = 0;
995  nedofs = 0;
996  nfdofs = 0;
997  nbdofs = 0;
998  fdofs = NULL;
999  bdofs = NULL;
1000 
1001  dynamic_cast<const NURBSFECollection *>(fec)->Reset();
1002 
1003  ndofs = NURBSext->GetNDof();
1004  elem_dof = NURBSext->GetElementDofTable();
1005  bdrElem_dof = NURBSext->GetBdrElementDofTable();
1006 }
1007 
1008 void FiniteElementSpace::Construct()
1009 {
1010  int i;
1011 
1012  elem_dof = NULL;
1013  bdrElem_dof = NULL;
1014 
1015  nvdofs = mesh->GetNV() * fec->DofForGeometry(Geometry::POINT);
1016 
1017  if ( mesh->Dimension() > 1 )
1018  {
1019  nedofs = mesh->GetNEdges() * fec->DofForGeometry(Geometry::SEGMENT);
1020  }
1021  else
1022  {
1023  nedofs = 0;
1024  }
1025 
1026  ndofs = 0;
1027  nfdofs = 0;
1028  nbdofs = 0;
1029  bdofs = NULL;
1030  fdofs = NULL;
1031  cP = NULL;
1032  cR = NULL;
1033  cP_is_set = false;
1034  T = NULL;
1035  own_T = true;
1036 
1037  if (mesh->Dimension() == 3 && mesh->GetNE())
1038  {
1039  // Here we assume that all faces in the mesh have the same base
1040  // geometry -- the base geometry of the 0-th face element.
1041  // The class Mesh assumes the same inside GetFaceBaseGeometry(...).
1042  // Thus we do not need to generate all the faces in the mesh
1043  // if we do not need them.
1044  int fdof = fec->DofForGeometry(mesh->GetFaceBaseGeometry(0));
1045  if (fdof > 0)
1046  {
1047  fdofs = new int[mesh->GetNFaces()+1];
1048  fdofs[0] = 0;
1049  for (i = 0; i < mesh->GetNFaces(); i++)
1050  {
1051  nfdofs += fdof;
1052  // nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(i));
1053  fdofs[i+1] = nfdofs;
1054  }
1055  }
1056  }
1057 
1058  bdofs = new int[mesh->GetNE()+1];
1059  bdofs[0] = 0;
1060  for (i = 0; i < mesh->GetNE(); i++)
1061  {
1062  nbdofs += fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
1063  bdofs[i+1] = nbdofs;
1064  }
1065 
1066  ndofs = nvdofs + nedofs + nfdofs + nbdofs;
1067 
1068  // Do not build elem_dof Table here: in parallel it has to be constructed
1069  // later.
1070 }
1071 
1072 void FiniteElementSpace::GetElementDofs (int i, Array<int> &dofs) const
1073 {
1074  if (elem_dof)
1075  {
1076  elem_dof -> GetRow (i, dofs);
1077  }
1078  else
1079  {
1080  Array<int> V, E, Eo, F, Fo;
1081  int k, j, nv, ne, nf, nb, nfd, nd;
1082  int *ind, dim;
1083 
1084  dim = mesh->Dimension();
1085  nv = fec->DofForGeometry(Geometry::POINT);
1086  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
1087  nb = fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
1088  if (nv > 0)
1089  {
1090  mesh->GetElementVertices(i, V);
1091  }
1092  if (ne > 0)
1093  {
1094  mesh->GetElementEdges(i, E, Eo);
1095  }
1096  nfd = 0;
1097  if (dim == 3)
1098  {
1099  if (fec->HasFaceDofs(mesh->GetElementBaseGeometry(i)))
1100  {
1101  mesh->GetElementFaces(i, F, Fo);
1102  for (k = 0; k < F.Size(); k++)
1103  {
1104  nfd += fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
1105  }
1106  }
1107  }
1108  nd = V.Size() * nv + E.Size() * ne + nfd + nb;
1109  dofs.SetSize(nd);
1110  if (nv > 0)
1111  {
1112  for (k = 0; k < V.Size(); k++)
1113  {
1114  for (j = 0; j < nv; j++)
1115  {
1116  dofs[k*nv+j] = V[k]*nv+j;
1117  }
1118  }
1119  nv *= V.Size();
1120  }
1121  if (ne > 0)
1122  {
1123  // if (dim > 1)
1124  for (k = 0; k < E.Size(); k++)
1125  {
1126  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
1127  for (j = 0; j < ne; j++)
1128  {
1129  if (ind[j] < 0)
1130  {
1131  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1132  }
1133  else
1134  {
1135  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1136  }
1137  }
1138  }
1139  }
1140  ne = nv + ne * E.Size();
1141  if (nfd > 0)
1142  // if (dim == 3)
1143  {
1144  for (k = 0; k < F.Size(); k++)
1145  {
1146  ind = fec->DofOrderForOrientation(mesh->GetFaceBaseGeometry(F[k]),
1147  Fo[k]);
1148  nf = fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
1149  for (j = 0; j < nf; j++)
1150  {
1151  if (ind[j] < 0)
1152  {
1153  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[F[k]]+(-1-ind[j]) );
1154  }
1155  else
1156  {
1157  dofs[ne+j] = nvdofs+nedofs+fdofs[F[k]]+ind[j];
1158  }
1159  }
1160  ne += nf;
1161  }
1162  }
1163  k = nvdofs + nedofs + nfdofs + bdofs[i];
1164  for (j = 0; j < nb; j++)
1165  {
1166  dofs[ne+j] = k + j;
1167  }
1168  }
1169 }
1170 
1171 const FiniteElement *FiniteElementSpace::GetFE(int i) const
1172 {
1173  const FiniteElement *FE =
1174  fec->FiniteElementForGeometry(mesh->GetElementBaseGeometry(i));
1175 
1176  if (NURBSext)
1177  {
1178  NURBSext->LoadFE(i, FE);
1179  }
1180 
1181  return FE;
1182 }
1183 
1184 void FiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
1185 {
1186  if (bdrElem_dof)
1187  {
1188  bdrElem_dof->GetRow(i, dofs);
1189  }
1190  else
1191  {
1192  Array<int> V, E, Eo;
1193  int k, j, nv, ne, nf, nd, iF, oF;
1194  int *ind, dim;
1195 
1196  dim = mesh->Dimension();
1197  nv = fec->DofForGeometry(Geometry::POINT);
1198  if (nv > 0)
1199  {
1200  mesh->GetBdrElementVertices(i, V);
1201  }
1202  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
1203  if (ne > 0)
1204  {
1205  mesh->GetBdrElementEdges(i, E, Eo);
1206  }
1207  nd = V.Size() * nv + E.Size() * ne;
1208  nf = (dim == 3) ? (fec->DofForGeometry(
1209  mesh->GetBdrElementBaseGeometry(i))) : (0);
1210  if (nf > 0)
1211  {
1212  nd += nf;
1213  mesh->GetBdrElementFace(i, &iF, &oF);
1214  }
1215  dofs.SetSize(nd);
1216  if (nv > 0)
1217  {
1218  for (k = 0; k < V.Size(); k++)
1219  {
1220  for (j = 0; j < nv; j++)
1221  {
1222  dofs[k*nv+j] = V[k]*nv+j;
1223  }
1224  }
1225  nv *= V.Size();
1226  }
1227  if (ne > 0)
1228  {
1229  // if (dim > 1)
1230  for (k = 0; k < E.Size(); k++)
1231  {
1232  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
1233  for (j = 0; j < ne; j++)
1234  {
1235  if (ind[j] < 0)
1236  {
1237  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1238  }
1239  else
1240  {
1241  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1242  }
1243  }
1244  }
1245  }
1246  if (nf > 0)
1247  // if (dim == 3)
1248  {
1249  ne = nv + ne * E.Size();
1250  ind = (fec->DofOrderForOrientation(
1251  mesh->GetBdrElementBaseGeometry(i), oF));
1252  for (j = 0; j < nf; j++)
1253  {
1254  if (ind[j] < 0)
1255  {
1256  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[iF]+(-1-ind[j]) );
1257  }
1258  else
1259  {
1260  dofs[ne+j] = nvdofs+nedofs+fdofs[iF]+ind[j];
1261  }
1262  }
1263  }
1264  }
1265 }
1266 
1267 void FiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs) const
1268 {
1269  int j, k, nv, ne, nf, nd, dim = mesh->Dimension();
1270  Array<int> V, E, Eo;
1271  const int *ind;
1272 
1273  // for 1D, 2D and 3D faces
1274  nv = fec->DofForGeometry(Geometry::POINT);
1275  ne = (dim > 1) ? fec->DofForGeometry(Geometry::SEGMENT) : 0;
1276  if (nv > 0)
1277  {
1278  mesh->GetFaceVertices(i, V);
1279  }
1280  if (ne > 0)
1281  {
1282  mesh->GetFaceEdges(i, E, Eo);
1283  }
1284  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
1285  nd = V.Size() * nv + E.Size() * ne + nf;
1286  dofs.SetSize(nd);
1287  if (nv > 0)
1288  {
1289  for (k = 0; k < V.Size(); k++)
1290  {
1291  for (j = 0; j < nv; j++)
1292  {
1293  dofs[k*nv+j] = V[k]*nv+j;
1294  }
1295  }
1296  }
1297  nv *= V.Size();
1298  if (ne > 0)
1299  {
1300  for (k = 0; k < E.Size(); k++)
1301  {
1302  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
1303  for (j = 0; j < ne; j++)
1304  {
1305  if (ind[j] < 0)
1306  {
1307  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
1308  }
1309  else
1310  {
1311  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
1312  }
1313  }
1314  }
1315  }
1316  ne = nv + ne * E.Size();
1317  if (nf > 0)
1318  {
1319  for (j = nvdofs+nedofs+fdofs[i], k = 0; k < nf; j++, k++)
1320  {
1321  dofs[ne+k] = j;
1322  }
1323  }
1324 }
1325 
1326 void FiniteElementSpace::GetEdgeDofs(int i, Array<int> &dofs) const
1327 {
1328  int j, k, nv, ne;
1329  Array<int> V;
1330 
1331  nv = fec->DofForGeometry(Geometry::POINT);
1332  if (nv > 0)
1333  {
1334  mesh->GetEdgeVertices(i, V);
1335  }
1336  ne = fec->DofForGeometry(Geometry::SEGMENT);
1337  dofs.SetSize(2*nv+ne);
1338  if (nv > 0)
1339  {
1340  for (k = 0; k < 2; k++)
1341  {
1342  for (j = 0; j < nv; j++)
1343  {
1344  dofs[k*nv+j] = V[k]*nv+j;
1345  }
1346  }
1347  }
1348  nv *= 2;
1349  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
1350  {
1351  dofs[nv+j] = k;
1352  }
1353 }
1354 
1355 void FiniteElementSpace::GetVertexDofs(int i, Array<int> &dofs) const
1356 {
1357  int j, nv;
1358 
1359  nv = fec->DofForGeometry(Geometry::POINT);
1360  dofs.SetSize(nv);
1361  for (j = 0; j < nv; j++)
1362  {
1363  dofs[j] = i*nv+j;
1364  }
1365 }
1366 
1367 void FiniteElementSpace::GetElementInteriorDofs (int i, Array<int> &dofs) const
1368 {
1369  int j, k, nb;
1370  nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (i));
1371  dofs.SetSize (nb);
1372  k = nvdofs + nedofs + nfdofs + bdofs[i];
1373  for (j = 0; j < nb; j++)
1374  {
1375  dofs[j] = k + j;
1376  }
1377 }
1378 
1379 void FiniteElementSpace::GetEdgeInteriorDofs (int i, Array<int> &dofs) const
1380 {
1381  int j, k, ne;
1382 
1383  ne = fec -> DofForGeometry (Geometry::SEGMENT);
1384  dofs.SetSize (ne);
1385  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
1386  {
1387  dofs[j] = k;
1388  }
1389 }
1390 
1391 void FiniteElementSpace::GetFaceInteriorDofs (int i, Array<int> &dofs) const
1392 {
1393  int j, k, nf;
1394 
1395  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
1396  dofs.SetSize (nf);
1397  if (nf > 0)
1398  {
1399  for (j = 0, k = nvdofs+nedofs+fdofs[i]; j < nf; j++, k++)
1400  {
1401  dofs[j] = k;
1402  }
1403  }
1404 }
1405 
1406 const FiniteElement *FiniteElementSpace::GetBE (int i) const
1407 {
1408  const FiniteElement *BE;
1409 
1410  switch ( mesh->Dimension() )
1411  {
1412  case 1:
1413  BE = fec->FiniteElementForGeometry(Geometry::POINT);
1414  break;
1415  case 2:
1416  BE = fec->FiniteElementForGeometry(Geometry::SEGMENT);
1417  break;
1418  case 3:
1419  default:
1420  BE = fec->FiniteElementForGeometry(
1421  mesh->GetBdrElementBaseGeometry(i));
1422  }
1423 
1424  if (NURBSext)
1425  {
1426  NURBSext->LoadBE(i, BE);
1427  }
1428 
1429  return BE;
1430 }
1431 
1432 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
1433 {
1434  const FiniteElement *fe;
1435 
1436  switch (mesh->Dimension())
1437  {
1438  case 1:
1439  fe = fec->FiniteElementForGeometry(Geometry::POINT);
1440  break;
1441  case 2:
1442  fe = fec->FiniteElementForGeometry(Geometry::SEGMENT);
1443  break;
1444  case 3:
1445  default:
1446  fe = fec->FiniteElementForGeometry(mesh->GetFaceBaseGeometry(i));
1447  }
1448 
1449  // if (NURBSext)
1450  // NURBSext->LoadFaceElement(i, fe);
1451 
1452  return fe;
1453 }
1454 
1455 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i) const
1456 {
1457  return fec->FiniteElementForGeometry(Geometry::SEGMENT);
1458 }
1459 
1460 const FiniteElement *FiniteElementSpace::GetTraceElement(
1461  int i, int geom_type) const
1462 {
1463  return fec->TraceFiniteElementForGeometry(geom_type);
1464 }
1465 
1466 FiniteElementSpace::~FiniteElementSpace()
1467 {
1468  Destroy();
1469 }
1470 
1471 void FiniteElementSpace::Destroy()
1472 {
1473  delete cR;
1474  delete cP;
1475  if (own_T) { delete T; }
1476 
1477  dof_elem_array.DeleteAll();
1478  dof_ldof_array.DeleteAll();
1479 
1480  if (NURBSext)
1481  {
1482  if (own_ext) { delete NURBSext; }
1483  }
1484  else
1485  {
1486  delete elem_dof;
1487  delete bdrElem_dof;
1488 
1489  delete [] bdofs;
1490  delete [] fdofs;
1491  }
1492 }
1493 
1494 void FiniteElementSpace::Update(bool want_transform)
1495 {
1496  if (mesh->GetSequence() == sequence)
1497  {
1498  return; // mesh and space are in sync, no-op
1499  }
1500  if (want_transform && mesh->GetSequence() != sequence + 1)
1501  {
1502  MFEM_ABORT("Error in update sequence. Space needs to be updated after "
1503  "each mesh modification.");
1504  }
1505  sequence = mesh->GetSequence();
1506 
1507  if (NURBSext)
1508  {
1509  UpdateNURBS();
1510  return;
1511  }
1512 
1513  Table* old_elem_dof= NULL;
1514  int old_ndofs;
1515 
1516  // save old DOF table
1517  if (want_transform)
1518  {
1519  old_elem_dof = elem_dof;
1520  elem_dof = NULL;
1521  old_ndofs = ndofs;
1522  }
1523 
1524  Destroy();
1525  Construct(); // sets T to NULL, own_T to true
1526  BuildElementToDofTable();
1527 
1528  if (want_transform)
1529  {
1530  // calculate appropriate GridFunction transformation
1531  switch (mesh->GetLastOperation())
1532  {
1533  case Mesh::REFINE:
1534  {
1535  T = RefinementMatrix(old_ndofs, old_elem_dof);
1536  break;
1537  }
1538 
1539  case Mesh::DEREFINE:
1540  {
1541  GetConformingInterpolation();
1542  T = DerefinementMatrix(old_ndofs, old_elem_dof);
1543  if (cP && cR)
1544  {
1545  T = new TripleProductOperator(cP, cR, T, false, false, true);
1546  }
1547  break;
1548  }
1549 
1550  default:
1551  break; // T stays NULL
1552  }
1553  }
1554 
1555  delete old_elem_dof;
1556 }
1557 
1558 void FiniteElementSpace::Save(std::ostream &out) const
1559 {
1560  out << "FiniteElementSpace\n"
1561  << "FiniteElementCollection: " << fec->Name() << '\n'
1562  << "VDim: " << vdim << '\n'
1563  << "Ordering: " << ordering << '\n';
1564 }
1565 
1566 
1567 void QuadratureSpace::Construct()
1568 {
1569  // protected method
1570  int offset = 0;
1571  const int num_elem = mesh->GetNE();
1572  element_offsets = new int[num_elem + 1];
1573  for (int g = 0; g < Geometry::NumGeom; g++)
1574  {
1575  int_rule[g] = NULL;
1576  }
1577  for (int i = 0; i < num_elem; i++)
1578  {
1579  element_offsets[i] = offset;
1580  int geom = mesh->GetElementBaseGeometry(i);
1581  if (int_rule[geom] == NULL)
1582  {
1583  int_rule[geom] = &IntRules.Get(geom, order);
1584  }
1585  offset += int_rule[geom]->GetNPoints();
1586  }
1587  element_offsets[num_elem] = size = offset;
1588 }
1589 
1590 QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
1591  : mesh(mesh_)
1592 {
1593  const char *msg = "invalid input stream";
1594  string ident;
1595 
1596  in >> ident; MFEM_VERIFY(ident == "QuadratureSpace", msg);
1597  in >> ident; MFEM_VERIFY(ident == "Type:", msg);
1598  in >> ident;
1599  if (ident == "default_quadrature")
1600  {
1601  in >> ident; MFEM_VERIFY(ident == "Order:", msg);
1602  in >> order;
1603  }
1604  else
1605  {
1606  MFEM_ABORT("unknown QuadratureSpace type: " << ident);
1607  return;
1608  }
1609 
1610  Construct();
1611 }
1612 
1613 void QuadratureSpace::Save(std::ostream &out) const
1614 {
1615  out << "QuadratureSpace\n"
1616  << "Type: default_quadrature\n"
1617  << "Order: " << order << '\n';
1618 }
1619 
1620 
1621 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
Abstract class for Finite Elements.
Definition: fe.hpp:47
int GetOrder() const
Definition: fe_coll.hpp:348
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:333
void Set(const double *p, const int dim)
Definition: intrules.hpp:32
int Size() const
Logical size of the array.
Definition: array.hpp:110
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:213
void Get(double *p, const int dim) const
Definition: intrules.hpp:45
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:1907
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
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Definition: fe.cpp:115
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:278
void InvertLinearTrans(IsoparametricTransformation &trans, const DenseMatrix &invdfdx, const IntegrationPoint &pt, Vector &x)
Definition: fespace.cpp:810
const Geometry::Type geom
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:260
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:633
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:166
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
int SizeK() const
Definition: densemat.hpp:660
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:161
void SetIdentityTransformation(int GeomType)
Definition: eltrans.cpp:370
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:376
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
void SetSize(int i, int j, int k)
Definition: densemat.hpp:662
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
Definition: ncmesh.cpp:1986
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe.cpp:109
T Sum()
Sum all entries.
Definition: array.cpp:151
STL namespace.
virtual void Transform(const IntegrationPoint &, Vector &)
Definition: eltrans.cpp:461
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:274
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:184
Geometry Geometries
Definition: geom.cpp:759
long GetSequence() const
Definition: mesh.hpp:1012
int dim
Definition: ex3.cpp:47
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:2043
DenseTensor point_matrices
matrices for IsoparametricTransformation
Definition: ncmesh.hpp:54
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:548
std::vector< Master > masters
Definition: ncmesh.hpp:169
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:67
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...
const IntegrationRule & GetNodes() const
Definition: fe.hpp:167
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:123
Array< Embedding > embeddings
fine element positions in their parents
Definition: ncmesh.hpp:55
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:523
void SetRow(int r, const Vector &row)
Definition: densemat.cpp:2817
int slaves_end
slave faces
Definition: ncmesh.hpp:144
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3156
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:119
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:63
std::vector< Slave > slaves
Definition: ncmesh.hpp:170
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:122
void mfem_error(const char *msg)
Definition: error.cpp:107
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:503
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Definition: fespace.cpp:1613
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
Definition: solvers.cpp:501
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:174
Class for integration point with weight.
Definition: intrules.hpp:25
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition: operator.hpp:346
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:298
void GetRow(int r, Vector &row)
Definition: densemat.cpp:2296
int parent
element index in the coarse mesh
Definition: ncmesh.hpp:45
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2082
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:644
Vector data type.
Definition: vector.hpp:41
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:143
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:296
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:52
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
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:3132
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:2456
int index
Mesh number.
Definition: ncmesh.hpp:132
BiLinear2DFiniteElement QuadrilateralFE
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:634
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:343
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
Defines the position of a fine element within a coarse element.
Definition: ncmesh.hpp:43
int matrix
index into CoarseFineTransformations::point_matrices
Definition: ncmesh.hpp:46