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