MFEM  v3.0
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
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
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 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs) const
30 {
31  int i, j, size;
32
33  if (vdim == 1) return;
34
35  size = dofs.Size();
36  dofs.SetSize (size * vdim);
37
38  switch(ordering)
39  {
40  case Ordering::byNODES:
41  for (i = 1; i < vdim; i++)
42  for (j = 0; j < size; j++)
43  if (dofs[j] < 0)
44  dofs[size * i + j] = -1 - ( ndofs * i + (-1-dofs[j]) );
45  else
46  dofs[size * i + j] = ndofs * i + dofs[j];
47  break;
48
49  case Ordering::byVDIM:
50  for (i = vdim-1; i >= 0; i--)
51  for (j = 0; j < size; j++)
52  if (dofs[j] < 0)
53  dofs[size * i + j] = -1 - ( (-1-dofs[j]) * vdim + i );
54  else
55  dofs[size * i + j] = dofs[j] * vdim + i;
56  break;
57  }
58 }
59
60 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs) const
61 {
62  if (vdim == 1)
63  return;
64  if (ordering == Ordering::byNODES)
65  {
66  for (int i = 0; i < dofs.Size(); i++)
67  {
68  int dof = dofs[i];
69  if (dof < 0)
70  dofs[i] = -1 - ((-1-dof) + vd * ndofs);
71  else
72  dofs[i] = dof + vd * ndofs;
73  }
74  }
75  else
76  {
77  for (int i = 0; i < dofs.Size(); i++)
78  {
79  int dof = dofs[i];
80  if (dof < 0)
81  dofs[i] = -1 - ((-1-dof) * vdim + vd);
82  else
83  dofs[i] = dof * vdim + vd;
84  }
85  }
86 }
87
88 int FiniteElementSpace::DofToVDof (int dof, int vd) const
89 {
90  if (vdim == 1)
91  return dof;
92  if (ordering == Ordering::byNODES)
93  {
94  if (dof < 0)
95  return -1 - ((-1-dof) + vd * ndofs);
96  else
97  return dof + vd * ndofs;
98  }
99  if (dof < 0)
100  return -1 - ((-1-dof) * vdim + vd);
101  else
102  return dof * vdim + vd;
103 }
104
105 // static function
107 {
108  int n = vdofs.Size(), *vdof = vdofs;
109  for (int i = 0; i < n; i++)
110  {
111  int j;
112  if ((j=vdof[i]) < 0)
113  vdof[i] = -1-j;
114  }
115 }
116
117 void FiniteElementSpace::GetElementVDofs(int iE, Array<int> &dofs) const
118 {
119  GetElementDofs(iE, dofs);
120  DofsToVDofs (dofs);
121 }
122
123 void FiniteElementSpace::GetBdrElementVDofs (int iE, Array<int> &dofs) const
124 {
125  GetBdrElementDofs(iE, dofs);
126  DofsToVDofs (dofs);
127 }
128
129 void FiniteElementSpace::GetFaceVDofs (int iF, Array<int> &dofs) const
130 {
131  GetFaceDofs (iF, dofs);
132  DofsToVDofs (dofs);
133 }
134
135 void FiniteElementSpace::GetEdgeVDofs (int iE, Array<int> &dofs) const
136 {
137  GetEdgeDofs (iE, dofs);
138  DofsToVDofs (dofs);
139 }
140
141 void FiniteElementSpace::GetElementInteriorVDofs (int i, Array<int> &vdofs)
142  const
143 {
144  GetElementInteriorDofs (i, vdofs);
145  DofsToVDofs (vdofs);
146 }
147
148 void FiniteElementSpace::GetEdgeInteriorVDofs (int i, Array<int> &vdofs)
149  const
150 {
151  GetEdgeInteriorDofs (i, vdofs);
152  DofsToVDofs (vdofs);
153 }
154
155 void FiniteElementSpace::BuildElementToDofTable()
156 {
157  if (elem_dof)
158  return;
159
160  Table *el_dof = new Table;
161  Array<int> dofs;
162  el_dof -> MakeI (mesh -> GetNE());
163  for (int i = 0; i < mesh -> GetNE(); i++)
164  {
165  GetElementDofs (i, dofs);
166  el_dof -> AddColumnsInRow (i, dofs.Size());
167  }
168  el_dof -> MakeJ();
169  for (int i = 0; i < mesh -> GetNE(); i++)
170  {
171  GetElementDofs (i, dofs);
172  el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
173  }
174  el_dof -> ShiftUpI();
175  elem_dof = el_dof;
176 }
177
178 void FiniteElementSpace::BuildDofToArrays()
179 {
180  if (dof_elem_array.Size())
181  return;
182  BuildElementToDofTable();
183
184  dof_elem_array.SetSize (ndofs);
185  dof_ldof_array.SetSize (ndofs);
186  dof_elem_array = -1;
187  for (int i = 0; i < mesh -> GetNE(); i++)
188  {
189  const int *dofs = elem_dof -> GetRow(i);
190  const int n = elem_dof -> RowSize(i);
191  for (int j = 0; j < n; j++)
192  if (dof_elem_array[dofs[j]] < 0)
193  {
194  dof_elem_array[dofs[j]] = i;
195  dof_ldof_array[dofs[j]] = j;
196  }
197  }
198 }
199
200 DenseMatrix * FiniteElementSpace::LocalInterpolation
201 (int k, int num_c_dofs, RefinementType type, Array<int> &rows)
202 {
203  int i,j,l=-1;
204
205  // generate refinement data if necessary
206  for (i = 0; i < RefData.Size(); i++)
207  if (RefData[i] -> type == type)
208  { l = i; break; }
209
210  if (l == -1)
211  { ConstructRefinementData (k,num_c_dofs,type); l = i; }
212
213  // determine the global indices of the fine dofs on the coarse element
214  Array<int> l_dofs, g_dofs;
215  int num_fine_elems = RefData[l] -> num_fine_elems;
216  rows.SetSize(RefData[l] -> num_fine_dofs);
217
218  for (i = 0; i < num_fine_elems; i++) {
219  GetElementDofs(mesh->GetFineElem(k,i),g_dofs); // TWO_LEVEL_FINE
220  RefData[l] -> fl_to_fc -> GetRow (i, l_dofs);
221  for (j = 0; j < l_dofs.Size(); j++)
222  rows[l_dofs[j]] = g_dofs[j];
223  }
224
225  return RefData[l] -> I;
226 }
227
228 // helper to set submatrix of A repeated vdim times
229 static void SetVDofSubMatrixTranspose(SparseMatrix& A,
230  Array<int>& rows, Array<int>& cols,
231  const DenseMatrix& subm, int vdim)
232 {
233  if (vdim == 1)
234  {
235  A.SetSubMatrixTranspose(rows, cols, subm, 1);
236  }
237  else
238  {
239  int nr = subm.Width(), nc = subm.Height();
240  for (int d = 0; d < vdim; d++)
241  {
242  Array<int> rows_sub(rows.GetData() + d*nr, nr); // (not owner)
243  Array<int> cols_sub(cols.GetData() + d*nc, nc); // (not owner)
244  A.SetSubMatrixTranspose(rows_sub, cols_sub, subm, 1);
245  }
246  }
247 }
248
249 SparseMatrix * FiniteElementSpace::GlobalRestrictionMatrix
250 (FiniteElementSpace * cfes, int one_vdim)
251 {
252  int k;
253  SparseMatrix * R;
254  DenseMatrix * r;
255  Array<int> rows, cols;
256
257  if (one_vdim == -1)
258  one_vdim = (ordering == Ordering::byNODES) ? 1 : 0;
259
260  if (mesh->ncmesh)
261  {
262  MFEM_VERIFY(vdim == 1 || one_vdim == 0,
263  "parameter 'one_vdim' must be 0 for nonconforming mesh.");
264  return NC_GlobalRestrictionMatrix(cfes, mesh->ncmesh);
265  }
266
267  mesh->SetState(Mesh::TWO_LEVEL_COARSE);
268  int vdim_or_1 = (one_vdim ? 1 : vdim);
269  R = new SparseMatrix(vdim_or_1 * cfes->GetNDofs(), vdim_or_1 * ndofs);
270
271  for (k = 0; k < mesh -> GetNE(); k++)
272  {
273  cfes->GetElementDofs(k, rows);
274
275  mesh->SetState(Mesh::TWO_LEVEL_FINE);
276  r = LocalInterpolation(k, rows.Size(), mesh->GetRefinementType(k), cols);
277
278  if (vdim_or_1 != 1)
279  {
280  cfes->DofsToVDofs(rows);
281  DofsToVDofs(cols);
282  }
283  SetVDofSubMatrixTranspose(*R, rows, cols, *r, vdim_or_1);
284
285  mesh->SetState(Mesh::TWO_LEVEL_COARSE);
286  }
287
288  return R;
289 }
290
291 SparseMatrix* FiniteElementSpace::NC_GlobalRestrictionMatrix
292 (FiniteElementSpace* cfes, NCMesh* ncmesh)
293 {
294  Array<int> rows, cols, rs, cs;
295  LinearFECollection linfec;
296
297  NCMesh::FineTransform* transforms = ncmesh->GetFineTransforms();
298
299  SparseMatrix* R = new SparseMatrix(cfes->GetVSize(), this->GetVSize());
300
301  // we mark each fine DOF the first time its column is set so that slave node
302  // values don't get represented twice in R
303  Array<int> mark(this->GetNDofs());
304  mark = 0;
305
306  // loop over the fine elements, get interpolations of the coarse elements
307  for (int k = 0; k < mesh->GetNE(); k++)
308  {
309  mesh->SetState(Mesh::TWO_LEVEL_COARSE);
310  cfes->GetElementDofs(transforms[k].coarse_index, rows);
311
312  mesh->SetState(Mesh::TWO_LEVEL_FINE);
313  this->GetElementDofs(k, cols);
314
315  if (!transforms[k].IsIdentity())
316  {
317  int geom = mesh->GetElementBaseGeometry(k);
318  const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
319
321  trans.SetFE(linfec.FiniteElementForGeometry(geom));
322  trans.GetPointMat() = transforms[k].point_matrix;
323
324  DenseMatrix I(fe->GetDof());
325  fe->GetLocalInterpolation(trans, I);
326  // TODO: use std::unordered_map to cache I matrices (point_matrix as key)
327
328  // make sure we don't set any column of R more than once
329  for (int i = 0; i < I.Height(); i++)
330  {
331  int col = cols[i];
332  if (col < 0)
333  col = -1 - col;
334  if (mark[col]++)
335  I.SetRow(i, 0); // zero the i-th row of I
336  }
337
338  cfes->DofsToVDofs(rows);
339  this->DofsToVDofs(cols);
340  SetVDofSubMatrixTranspose(*R, rows, cols, I, vdim);
341  }
342  else // optimization: insert identity for elements that were not refined
343  {
344  MFEM_ASSERT(rows.Size() == cols.Size(), "");
345  for (int i = 0; i < rows.Size(); i++)
346  {
347  int col = cols[i];
348  if (col < 0)
349  col = -1 - col;
350  if (!mark[col]++)
351  for (int vd = 0; vd < vdim; vd++)
352  R->Set(cfes->DofToVDof(rows[i], vd),
353  this->DofToVDof(cols[i], vd), 1.0);
354  }
355  }
356  }
357
358  delete [] transforms;
359  return R;
360 }
361
362 void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
363  Array<int> &ess_dofs) const
364 {
365  int i, j, k;
366  Array<int> vdofs;
367
368  ess_dofs.SetSize(GetVSize());
369  ess_dofs = 0;
370
371  for (i = 0; i < GetNBE(); i++)
372  if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
373  {
374  GetBdrElementVDofs(i, vdofs);
375  for (j = 0; j < vdofs.Size(); j++)
376  if ( (k = vdofs[j]) >= 0 )
377  ess_dofs[k] = -1;
378  else
379  ess_dofs[-1-k] = -1;
380  }
381 }
382
383 void FiniteElementSpace::MarkDependency(const SparseMatrix *D,
384  const Array<int> &row_marker,
385  Array<int> &col_marker)
386 {
387  if (D)
388  {
389  col_marker.SetSize(D->Width());
390  col_marker = 0;
391
392  for (int i = 0; i < D->Height(); i++)
393  if (row_marker[i] < 0)
394  {
395  const int *col = D->GetRowColumns(i), n = D->RowSize(i);
396  for (int j = 0; j < n; j++)
397  col_marker[col[j]] = -1;
398  }
399  }
400  else
401  {
402  row_marker.Copy(col_marker);
403  }
404 }
405
406 void FiniteElementSpace::EliminateEssentialBCFromGRM
407 (FiniteElementSpace *cfes, Array<int> &bdr_attr_is_ess, SparseMatrix *R)
408 {
409  int i, j, k, one_vdim;
410  Array<int> dofs;
411
412  one_vdim = (cfes -> GetNDofs() == R -> Height()) ? 1 : 0;
413
414  mesh -> SetState (Mesh::TWO_LEVEL_COARSE);
415  if (bdr_attr_is_ess.Size() != 0)
416  for (i=0; i < cfes -> GetNBE(); i++)
417  if (bdr_attr_is_ess[cfes -> GetBdrAttribute (i)-1])
418  {
419  if (one_vdim == 1)
420  cfes -> GetBdrElementDofs (i, dofs);
421  else
422  cfes -> GetBdrElementVDofs (i, dofs);
423  for (j=0; j < dofs.Size(); j++)
424  if ( (k = dofs[j]) >= 0 )
425  R -> EliminateRow(k);
426  else
427  R -> EliminateRow(-1-k);
428  }
429  R -> Finalize();
430 }
431
432 SparseMatrix * FiniteElementSpace::GlobalRestrictionMatrix
433 (FiniteElementSpace * cfes, Array<int> &bdr_attr_is_ess, int one_vdim)
434 {
435  SparseMatrix * R;
436
437  R = GlobalRestrictionMatrix (cfes, one_vdim);
438  EliminateEssentialBCFromGRM (cfes, bdr_attr_is_ess, R);
439
440  return R;
441 }
442
443 SparseMatrix *
444 FiniteElementSpace::D2C_GlobalRestrictionMatrix (FiniteElementSpace *cfes)
445 {
446  int i, j;
447  Array<int> d_vdofs, c_vdofs;
448  SparseMatrix *R;
449
450  R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
451
452  for (i = 0; i < mesh -> GetNE(); i++)
453  {
454  this -> GetElementVDofs (i, d_vdofs);
455  cfes -> GetElementVDofs (i, c_vdofs);
456
457 #ifdef MFEM_DEBUG
458  if (d_vdofs.Size() != c_vdofs.Size())
459  mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
460 #endif
461
462  for (j = 0; j < d_vdofs.Size(); j++)
463  R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
464  }
465
466  R -> Finalize();
467
468  return R;
469 }
470
471 SparseMatrix *
472 FiniteElementSpace::D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
473 {
474  int i, j;
475  Array<int> d_dofs, c_dofs;
476  SparseMatrix *R;
477
478  R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
479
480  for (i = 0; i < mesh -> GetNE(); i++)
481  {
482  this -> GetElementDofs (i, d_dofs);
483  cfes -> GetElementDofs (i, c_dofs);
484
485 #ifdef MFEM_DEBUG
486  if (c_dofs.Size() != 1)
487  mfem_error ("FiniteElementSpace::"
488  "D2Const_GlobalRestrictionMatrix (...)");
489 #endif
490
491  for (j = 0; j < d_dofs.Size(); j++)
492  R -> Set (c_dofs[0], d_dofs[j], 1.0);
493  }
494
495  R -> Finalize();
496
497  return R;
498 }
499
500 SparseMatrix *
501 FiniteElementSpace::H2L_GlobalRestrictionMatrix (FiniteElementSpace *lfes)
502 {
503  SparseMatrix *R;
504  DenseMatrix loc_restr;
505  Array<int> l_dofs, h_dofs;
506
507  R = new SparseMatrix (lfes -> GetNDofs(), ndofs);
508
509  if (!lfes->GetNE())
510  {
511  R->Finalize();
512  return R;
513  }
514
515  const FiniteElement *h_fe = this -> GetFE (0);
516  const FiniteElement *l_fe = lfes -> GetFE (0);
519  h_fe->Project(*l_fe, T, loc_restr);
520
521  for (int i = 0; i < mesh -> GetNE(); i++)
522  {
523  this -> GetElementDofs (i, h_dofs);
524  lfes -> GetElementDofs (i, l_dofs);
525
526  R -> SetSubMatrix (l_dofs, h_dofs, loc_restr, 1);
527  }
528
529  R -> Finalize();
530
531  return R;
532 }
533
534 FiniteElementSpace::FiniteElementSpace(FiniteElementSpace &fes)
535 {
536  mesh = fes.mesh;
537  vdim = fes.vdim;
538  ndofs = fes.ndofs;
539  ordering = fes.ordering;
540  fec = fes.fec;
541  nvdofs = fes.nvdofs;
542  nedofs = fes.nedofs;
543  nfdofs = fes.nfdofs;
544  nbdofs = fes.nbdofs;
545  fdofs = fes.fdofs;
546  bdofs = fes.bdofs;
547  // keep 'RefData' in 'fes'
548  elem_dof = fes.elem_dof;
549  bdrElem_dof = fes.bdrElem_dof;
550  Swap(dof_elem_array, fes.dof_elem_array);
551  Swap(dof_ldof_array, fes.dof_ldof_array);
552
553  NURBSext = fes.NURBSext;
554  own_ext = 0;
555
556  cP = fes.cP;
557  cR = fes.cR;
558  fes.cP = NULL;
559  fes.cR = NULL;
560
561  fes.bdofs = NULL;
562  fes.fdofs = NULL;
563  fes.elem_dof = NULL;
564  fes.bdrElem_dof = NULL;
565 }
566
567 FiniteElementSpace::FiniteElementSpace(Mesh *m,
568  const FiniteElementCollection *f,
569  int dim, int order)
570 {
571  mesh = m;
572  fec = f;
573  vdim = dim;
574  ordering = order;
575
576  const NURBSFECollection *nurbs_fec =
577  dynamic_cast<const NURBSFECollection *>(fec);
578  if (nurbs_fec)
579  {
580  if (!mesh->NURBSext)
581  {
582  mfem_error("FiniteElementSpace::FiniteElementSpace :\n"
583  " NURBS FE space requires NURBS mesh.");
584  }
585  else
586  {
587  int Order = nurbs_fec->GetOrder();
588  if (mesh->NURBSext->GetOrder() == Order)
589  {
590  NURBSext = mesh->NURBSext;
591  own_ext = 0;
592  }
593  else
594  {
595  NURBSext = new NURBSExtension(mesh->NURBSext, Order);
596  own_ext = 1;
597  }
598  UpdateNURBS();
599  cP = cR = NULL;
600  }
601  }
602  else
603  {
604  NURBSext = NULL;
605  own_ext = 0;
606  Constructor();
607  }
608 }
609
610 NURBSExtension *FiniteElementSpace::StealNURBSext()
611 {
612  if (NURBSext && !own_ext)
613  mfem_error("FiniteElementSpace::StealNURBSext");
614  own_ext = 0;
615
616  return NURBSext;
617 }
618
619 void FiniteElementSpace::UpdateNURBS()
620 {
621  nvdofs = 0;
622  nedofs = 0;
623  nfdofs = 0;
624  nbdofs = 0;
625  fdofs = NULL;
626  bdofs = NULL;
627
628  dynamic_cast<const NURBSFECollection *>(fec)->Reset();
629
630  ndofs = NURBSext->GetNDof();
631
632  elem_dof = NURBSext->GetElementDofTable();
633
634  bdrElem_dof = NURBSext->GetBdrElementDofTable();
635 }
636
637 void FiniteElementSpace::Constructor()
638 {
639  int i;
640
641  elem_dof = NULL;
642  bdrElem_dof = NULL;
643
644  nvdofs = mesh->GetNV() * fec->DofForGeometry(Geometry::POINT);
645
646  if ( mesh->Dimension() > 1 )
647  nedofs = mesh->GetNEdges() * fec->DofForGeometry(Geometry::SEGMENT);
648  else
649  nedofs = 0;
650
651  ndofs = 0;
652  nfdofs = 0;
653  nbdofs = 0;
654  bdofs = NULL;
655  fdofs = NULL;
656  cP = NULL;
657  cR = NULL;
658
659  if (!mesh->GetNE())
660  return;
661
662  if (mesh->Dimension() == 3)
663  {
664  // Here we assume that all faces in the mesh have the same base
665  // geometry -- the base geometry of the 0-th face element.
666  // The class Mesh assumes the same inside GetFaceBaseGeometry(...).
667  // Thus we do not need to generate all the faces in the mesh
668  // if we do not need them.
669  int fdof = fec->DofForGeometry(mesh->GetFaceBaseGeometry(0));
670  if (fdof > 0)
671  {
672  fdofs = new int[mesh->GetNFaces()+1];
673  fdofs[0] = 0;
674  for (i = 0; i < mesh->GetNFaces(); i++)
675  {
676  nfdofs += fdof;
677  // nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(i));
678  fdofs[i+1] = nfdofs;
679  }
680  }
681  }
682
683  bdofs = new int[mesh->GetNE()+1];
684  bdofs[0] = 0;
685  for (i = 0; i < mesh->GetNE(); i++)
686  {
687  nbdofs += fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
688  bdofs[i+1] = nbdofs;
689  }
690
691  ndofs = nvdofs + nedofs + nfdofs + nbdofs;
692
693  if (mesh->ncmesh && ndofs > nbdofs)
694  {
695  cP = mesh->ncmesh->GetInterpolation(this, &cR);
696  if (cP && vdim > 1)
697  {
698  Array<int> cdofs, vcdofs;
699  Vector srow;
700  SparseMatrix *vec_cP =
701  new SparseMatrix(vdim*cP->Height(), vdim*cP->Width());
702  for (int i = 0; i < cP->Height(); i++)
703  {
704  cP->GetRow(i, cdofs, srow);
705  for (int vd = 0; vd < vdim; vd++)
706  {
707  cdofs.Copy(vcdofs);
708  ndofs = cP->Width(); // make DofsToVDofs work on conf. dofs
709  DofsToVDofs(vd, vcdofs);
710  ndofs = cP->Height();
711  vec_cP->SetRow(DofToVDof(i, vd), vcdofs, srow);
712  }
713  }
714  delete cP;
715  vec_cP->Finalize();
716  cP = vec_cP;
717
718  SparseMatrix *vec_cR =
719  new SparseMatrix(vdim*cR->Height(), vdim*cR->Width());
720  for (int i = 0; i < cR->Height(); i++)
721  {
722  cR->GetRow(i, cdofs, srow); // here, cdofs are partially conf. dofs
723  for (int vd = 0; vd < vdim; vd++)
724  {
725  cdofs.Copy(vcdofs); // here, vcdofs are partially conf. dofs
726  DofsToVDofs(vd, vcdofs);
727  ndofs = cR->Height(); // make DofToVDof work on conf. dofs
728  vec_cR->SetRow(DofToVDof(i, vd), vcdofs, srow);
729  ndofs = cR->Width();
730  }
731  }
732  delete cR;
733  vec_cR->Finalize();
734  cR = vec_cR;
735  }
736  }
737 }
738
739 void FiniteElementSpace::GetElementDofs (int i, Array<int> &dofs) const
740 {
741  if (elem_dof)
742  {
743  elem_dof -> GetRow (i, dofs);
744  }
745  else
746  {
747  Array<int> V, E, Eo, F, Fo;
748  int k, j, nv, ne, nf, nb, nfd, nd;
749  int *ind, dim;
750
751  dim = mesh->Dimension();
752  nv = fec->DofForGeometry(Geometry::POINT);
753  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
754  nb = fec->DofForGeometry(mesh->GetElementBaseGeometry(i));
755  if (nv > 0)
756  mesh->GetElementVertices(i, V);
757  if (ne > 0)
758  mesh->GetElementEdges(i, E, Eo);
759  nfd = 0;
760  if (dim == 3)
761  if (fec->HasFaceDofs(mesh->GetElementBaseGeometry(i)))
762  {
763  mesh->GetElementFaces(i, F, Fo);
764  for (k = 0; k < F.Size(); k++)
765  nfd += fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
766  }
767  nd = V.Size() * nv + E.Size() * ne + nfd + nb;
768  dofs.SetSize(nd);
769  if (nv > 0)
770  {
771  for (k = 0; k < V.Size(); k++)
772  for (j = 0; j < nv; j++)
773  dofs[k*nv+j] = V[k]*nv+j;
774  nv *= V.Size();
775  }
776  if (ne > 0)
777  // if (dim > 1)
778  for (k = 0; k < E.Size(); k++)
779  {
780  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
781  for (j = 0; j < ne; j++)
782  if (ind[j] < 0)
783  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
784  else
785  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
786  }
787  ne = nv + ne * E.Size();
788  if (nfd > 0)
789  // if (dim == 3)
790  {
791  for (k = 0; k < F.Size(); k++)
792  {
793  ind = fec->DofOrderForOrientation(mesh->GetFaceBaseGeometry(F[k]),
794  Fo[k]);
795  nf = fec->DofForGeometry(mesh->GetFaceBaseGeometry(F[k]));
796  for (j = 0; j < nf; j++)
797  {
798  if (ind[j] < 0)
799  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[F[k]]+(-1-ind[j]) );
800  else
801  dofs[ne+j] = nvdofs+nedofs+fdofs[F[k]]+ind[j];
802  }
803  ne += nf;
804  }
805  }
806  k = nvdofs + nedofs + nfdofs + bdofs[i];
807  for (j = 0; j < nb; j++)
808  {
809  dofs[ne+j] = k + j;
810  }
811  }
812 }
813
814 const FiniteElement *FiniteElementSpace::GetFE(int i) const
815 {
816  const FiniteElement *FE =
817  fec->FiniteElementForGeometry(mesh->GetElementBaseGeometry(i));
818
819  if (NURBSext)
821
822  return FE;
823 }
824
825 void FiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
826 {
827  if (bdrElem_dof)
828  {
829  bdrElem_dof->GetRow(i, dofs);
830  }
831  else
832  {
833  Array<int> V, E, Eo;
834  int k, j, nv, ne, nf, nd, iF, oF;
835  int *ind, dim;
836
837  dim = mesh->Dimension();
838  nv = fec->DofForGeometry(Geometry::POINT);
839  if (nv > 0)
840  mesh->GetBdrElementVertices(i, V);
841  ne = (dim > 1) ? ( fec->DofForGeometry(Geometry::SEGMENT) ) : ( 0 );
842  if (ne > 0)
843  mesh->GetBdrElementEdges(i, E, Eo);
844  nd = V.Size() * nv + E.Size() * ne;
845  nf = (dim == 3) ? (fec->DofForGeometry(
846  mesh->GetBdrElementBaseGeometry(i))) : (0);
847  if (nf > 0)
848  {
849  nd += nf;
850  mesh->GetBdrElementFace(i, &iF, &oF);
851  }
852  dofs.SetSize(nd);
853  if (nv > 0)
854  {
855  for (k = 0; k < V.Size(); k++)
856  for (j = 0; j < nv; j++)
857  dofs[k*nv+j] = V[k]*nv+j;
858  nv *= V.Size();
859  }
860  if (ne > 0)
861  // if (dim > 1)
862  for (k = 0; k < E.Size(); k++)
863  {
864  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
865  for (j = 0; j < ne; j++)
866  if (ind[j] < 0)
867  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
868  else
869  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
870  }
871  if (nf > 0)
872  // if (dim == 3)
873  {
874  ne = nv + ne * E.Size();
875  ind = (fec->DofOrderForOrientation(
876  mesh->GetBdrElementBaseGeometry(i), oF));
877  for (j = 0; j < nf; j++)
878  {
879  if (ind[j] < 0)
880  dofs[ne+j] = -1 - ( nvdofs+nedofs+fdofs[iF]+(-1-ind[j]) );
881  else
882  dofs[ne+j] = nvdofs+nedofs+fdofs[iF]+ind[j];
883  }
884  }
885  }
886 }
887
888 void FiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs) const
889 {
890  int j, k, nv, ne, nf, nd, dim = mesh->Dimension();
891  Array<int> V, E, Eo;
892  const int *ind;
893
894  // for 1D, 2D and 3D faces
895  nv = fec->DofForGeometry(Geometry::POINT);
896  ne = (dim > 1) ? fec->DofForGeometry(Geometry::SEGMENT) : 0;
897  if (nv > 0)
898  mesh->GetFaceVertices(i, V);
899  if (ne > 0)
900  mesh->GetFaceEdges(i, E, Eo);
901  nf = (fdofs) ? (fdofs[i+1]-fdofs[i]) : (0);
902  nd = V.Size() * nv + E.Size() * ne + nf;
903  dofs.SetSize(nd);
904  if (nv > 0)
905  for (k = 0; k < V.Size(); k++)
906  for (j = 0; j < nv; j++)
907  dofs[k*nv+j] = V[k]*nv+j;
908  nv *= V.Size();
909  if (ne > 0)
910  for (k = 0; k < E.Size(); k++)
911  {
912  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[k]);
913  for (j = 0; j < ne; j++)
914  if (ind[j] < 0)
915  dofs[nv+k*ne+j] = -1 - ( nvdofs+E[k]*ne+(-1-ind[j]) );
916  else
917  dofs[nv+k*ne+j] = nvdofs+E[k]*ne+ind[j];
918  }
919  ne = nv + ne * E.Size();
920  if (nf > 0)
921  for (j = nvdofs+nedofs+fdofs[i], k = 0; k < nf; j++, k++)
922  dofs[ne+k] = j;
923 }
924
925 void FiniteElementSpace::GetEdgeDofs(int i, Array<int> &dofs) const
926 {
927  int j, k, nv, ne;
928  Array<int> V;
929
930  nv = fec->DofForGeometry(Geometry::POINT);
931  if (nv > 0)
932  mesh->GetEdgeVertices(i, V);
933  ne = fec->DofForGeometry(Geometry::SEGMENT);
934  dofs.SetSize(2*nv+ne);
935  if (nv > 0)
936  for (k = 0; k < 2; k++)
937  for (j = 0; j < nv; j++)
938  dofs[k*nv+j] = V[k]*nv+j;
939  nv *= 2;
940  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
941  dofs[nv+j] = k;
942 }
943
944 void FiniteElementSpace::GetVertexDofs(int i, Array<int> &dofs) const
945 {
946  int j, nv;
947
948  nv = fec->DofForGeometry(Geometry::POINT);
949  dofs.SetSize(nv);
950  for (j = 0; j < nv; j++)
951  dofs[j] = i*nv+j;
952 }
953
954 void FiniteElementSpace::GetElementInteriorDofs (int i, Array<int> &dofs) const
955 {
956  int j, k, nb;
957  nb = fec -> DofForGeometry (mesh -> GetElementBaseGeometry (i));
958  dofs.SetSize (nb);
959  k = nvdofs + nedofs + nfdofs + bdofs[i];
960  for (j = 0; j < nb; j++)
961  {
962  dofs[j] = k + j;
963  }
964 }
965
966 void FiniteElementSpace::GetEdgeInteriorDofs (int i, Array<int> &dofs) const
967 {
968  int j, k, ne;
969
970  ne = fec -> DofForGeometry (Geometry::SEGMENT);
971  dofs.SetSize (ne);
972  for (j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
973  dofs[j] = k;
974 }
975
976 const FiniteElement *FiniteElementSpace::GetBE (int i) const
977 {
978  const FiniteElement *BE;
979
980  switch ( mesh->Dimension() )
981  {
982  case 1:
983  BE = fec->FiniteElementForGeometry(Geometry::POINT);
984  case 2:
985  BE = fec->FiniteElementForGeometry(Geometry::SEGMENT);
986  case 3:
987  default:
988  BE = fec->FiniteElementForGeometry(
989  mesh->GetBdrElementBaseGeometry(i));
990  }
991
992  if (NURBSext)
994
995  return BE;
996 }
997
998 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
999 {
1000  const FiniteElement *fe;
1001
1002  switch (mesh->Dimension())
1003  {
1004  case 1:
1005  fe = fec->FiniteElementForGeometry(Geometry::POINT);
1006  case 2:
1007  fe = fec->FiniteElementForGeometry(Geometry::SEGMENT);
1008  case 3:
1009  default:
1010  fe = fec->FiniteElementForGeometry(mesh->GetFaceBaseGeometry(i));
1011  }
1012
1013  // if (NURBSext)
1015
1016  return fe;
1017 }
1018
1019 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i) const
1020 {
1021  return fec->FiniteElementForGeometry(Geometry::SEGMENT);
1022 }
1023
1024 const FiniteElement *FiniteElementSpace::GetTraceElement(
1025  int i, int geom_type) const
1026 {
1027  return fec->TraceFiniteElementForGeometry(geom_type);
1028 }
1029
1030 FiniteElementSpace::~FiniteElementSpace()
1031 {
1032  Destructor();
1033  // delete RefData
1034  for (int i = 0; i < RefData.Size(); i++)
1035  delete RefData[i];
1036 }
1037
1038 void FiniteElementSpace::Destructor()
1039 {
1040  delete cR;
1041  delete cP;
1042
1043  dof_elem_array.DeleteAll();
1044  dof_ldof_array.DeleteAll();
1045
1046  if (NURBSext)
1047  {
1048  if (own_ext) delete NURBSext;
1049  }
1050  else
1051  {
1052  delete elem_dof;
1053  delete bdrElem_dof;
1054
1055  delete [] bdofs;
1056  delete [] fdofs;
1057  }
1058 }
1059
1061 {
1062  if (NURBSext)
1063  {
1064  UpdateNURBS();
1065  }
1066  else
1067  {
1068  Destructor(); // keeps RefData
1069  Constructor();
1070  }
1071 }
1072
1073 FiniteElementSpace *FiniteElementSpace::SaveUpdate()
1074 {
1075  FiniteElementSpace *cfes = new FiniteElementSpace(*this);
1076  Constructor();
1077  return cfes;
1078 }
1079
1080 void FiniteElementSpace::UpdateAndInterpolate(int num_grid_fns, ...)
1081 {
1082  if (mesh->GetState() == Mesh::NORMAL)
1083  {
1084  MFEM_ABORT("Mesh must be in two-level state, please call "
1085  "Mesh::UseTwoLevelState before refining.");
1086  }
1087
1088  FiniteElementSpace *cfes = SaveUpdate();
1089
1090  // obtain the (transpose of) interpolation matrix between mesh levels
1091  SparseMatrix *R = GlobalRestrictionMatrix(cfes, 0);
1092
1093  delete cfes;
1094
1095  // interpolate the grid functions
1096  std::va_list vl;
1097  va_start(vl, num_grid_fns);
1098  for (int i = 0; i < num_grid_fns; i++)
1099  {
1100  GridFunction* gf = va_arg(vl, GridFunction*);
1101  if (gf->FESpace() != this)
1102  {
1103  MFEM_ABORT("Cannot interpolate: grid function is not based "
1104  "on this space.");
1105  }
1106
1107  Vector coarse_gf = *gf;
1108  gf->Update();
1109  R->MultTranspose(coarse_gf, *gf);
1110  }
1111  va_end(vl);
1112
1113  delete R;
1114  mesh->SetState(Mesh::TWO_LEVEL_FINE);
1115 }
1116
1117 void FiniteElementSpace::ConstructRefinementData (int k, int num_c_dofs,
1118  RefinementType type)
1119 {
1120  int i,j,l;
1121  Array<int> dofs, g_dofs;
1122
1123  RefinementData * data = new RefinementData;
1124
1125  data -> type = type;
1126  data -> num_fine_elems = mesh -> GetNumFineElems(k);
1127
1128  // we assume that each fine element has <= dofs than the
1129  // initial coarse element
1130  data -> fl_to_fc = new Table(data -> num_fine_elems, num_c_dofs);
1131  for (i = 0; i < data -> num_fine_elems; i++) {
1132  GetElementDofs(mesh -> GetFineElem(k,i), g_dofs); // TWO_LEVEL_FINE
1133  for (j = 0; j < g_dofs.Size(); j++)
1134  data -> fl_to_fc -> Push (i,dofs.Union(g_dofs[j]));
1135  }
1136  data -> fl_to_fc -> Finalize();
1137  data -> num_fine_dofs = dofs.Size();
1138
1139  // construction of I
1140
1141  // k is a coarse element index but mesh is in TWO_LEVEL_FINE state
1142  int geomtype = mesh->GetElementBaseGeometry(k);
1143  const FiniteElement *fe = fec -> FiniteElementForGeometry(geomtype);
1144  // const IntegrationRule &ir = fe -> GetNodes();
1145  int nedofs = fe -> GetDof(); // number of dofs for each element
1146
1147  ElementTransformation *trans;
1148  // IntegrationPoint ip;
1149  // DenseMatrix tr;
1150  Array<int> row;
1151
1152  // Vector shape(nedofs);
1153  DenseMatrix I (nedofs);
1154  data -> I = new DenseMatrix(data -> num_fine_dofs, nedofs);
1155
1156  for (i=0; i < data -> num_fine_elems; i++) {
1157
1158  trans = mesh -> GetFineElemTrans(k,i);
1159  // trans -> Transform(ir,tr);
1160  fe -> GetLocalInterpolation (*trans, I);
1161  data -> fl_to_fc -> GetRow(i,row);
1162
1163  for (j=0; j < nedofs; j++)
1164  {
1165  /*
1166  ip.x = tr(0,j); ip.y = tr(1,j); if (tr.Height()==3) ip.z = tr(2,j);
1167  fe -> SetCalcElemTrans(trans); shape(0) = j;
1168  fe -> CalcShape(ip, shape);
1169  */
1170  for (l=0; l < nedofs; l++)
1171  (*(data->I))(row[j],l) = I(j,l);
1172  }
1173  }
1174
1175  RefData.Append(data);
1176 }
1177
1178 void FiniteElementSpace::Save(std::ostream &out) const
1179 {
1180  out << "FiniteElementSpace\n"
1181  << "FiniteElementCollection: " << fec->Name() << '\n'
1182  << "VDim: " << vdim << '\n'
1183  << "Ordering: " << ordering << '\n';
1184 }
1185
1186 }
Abstract class for Finite Elements.
Definition: fe.hpp:42
int GetOrder() const
Definition: fe_coll.hpp:195
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:180
int Size() const
Logical size of the array.
Definition: array.hpp:108
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:91
int GetVSize() const
Definition: fespace.hpp:151
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:149
int ndofs
Number of degrees of freedom. Number of unknowns are ndofs*vdim.
Definition: fespace.hpp:71
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
FineTransform * GetFineTransforms()
Definition: ncmesh.cpp:1870
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:132
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:155
T * GetData()
Returns the data.
Definition: array.hpp:90
void SetIdentityTransformation(int GeomType)
Definition: eltrans.cpp:29
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:221
Data type dense matrix.
Definition: densemat.hpp:22
int vdim
Vector dimension (number of unknowns per degree of freedom).
Definition: fespace.hpp:68
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe.cpp:69
int RefinementType
Type of refinement (int, a tree, etc.)
Definition: fespace.hpp:39
const FiniteElementCollection * fec
Definition: fespace.hpp:78
Array< int > dof_elem_array
Definition: fespace.hpp:87
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:169
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Definition: sparsemat.cpp:1661
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
SparseMatrix * cR
Definition: fespace.hpp:97
DenseMatrix point_matrix
for use in IsoparametricTransformation
Definition: ncmesh.hpp:101
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:739
Array< int > dof_ldof_array
Definition: fespace.hpp:87
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:118
SparseMatrix * cP
Definition: fespace.hpp:95
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
Definition: ncmesh.hpp:65
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:246
virtual void Finalize(int skip_zeros=1)
Definition: sparsemat.cpp:529
int Union(const T &el)
Append element when it is not yet in the array, return index.
Definition: array.hpp:375
int GetGeomType() const
Returns the geometry type:
Definition: fe.hpp:83
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:284
Data kept for every type of refinement.
Definition: fespace.hpp:42
Abstract finite element space.
Definition: fespace.hpp:61
int GetDof() const
Returns the degrees of freedom in the FE space.
Definition: fe.hpp:86
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1577
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:293
void Set(const int i, const int j, const double a)
Definition: sparsemat.cpp:1506
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
Definition: solvers.cpp:468
Mesh * mesh
The mesh that FE space lives on.
Definition: fespace.hpp:65
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:67
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:387
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:1700
int DofToVDof(int dof, int vd) const
Definition: fespace.cpp:88
void DofsToVDofs(Array< int > &dofs) const
Definition: fespace.cpp:29
Vector data type.
Definition: vector.hpp:29
NURBSExtension * NURBSext
Definition: fespace.hpp:89
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1667
DenseMatrix * I
Local interpolation matrix.
Definition: fespace.hpp:53