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