1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 // Implementation of FiniteElementSpace
13
14 #include "../general/text.hpp"
15 #include "../general/forall.hpp"
16 #include "../mesh/mesh_headers.hpp"
17 #include "fem.hpp"
18 #include "ceed/util.hpp"
19
20 #include <cmath>
21 #include <cstdarg>
22
23 using namespace std;
24
25 namespace mfem
26 {
27
28 template <> void Ordering::
DofsToVDofs(int ndofs,int vdim,Array<int> & dofs)29 DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
30 {
31 // static method
32 int size = dofs.Size();
33 dofs.SetSize(size*vdim);
34 for (int vd = 1; vd < vdim; vd++)
35 {
36 for (int i = 0; i < size; i++)
37 {
38 dofs[i+size*vd] = Map<byNODES>(ndofs, vdim, dofs[i], vd);
39 }
40 }
41 }
42
43 template <> void Ordering::
DofsToVDofs(int ndofs,int vdim,Array<int> & dofs)44 DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
45 {
46 // static method
47 int size = dofs.Size();
48 dofs.SetSize(size*vdim);
49 for (int vd = vdim-1; vd >= 0; vd--)
50 {
51 for (int i = 0; i < size; i++)
52 {
53 dofs[i+size*vd] = Map<byVDIM>(ndofs, vdim, dofs[i], vd);
54 }
55 }
56 }
57
58
FiniteElementSpace()59 FiniteElementSpace::FiniteElementSpace()
60 : mesh(NULL), fec(NULL), vdim(0), ordering(Ordering::byNODES),
61 ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0), bdofs(NULL),
62 elem_dof(NULL), bdr_elem_dof(NULL), face_dof(NULL),
63 NURBSext(NULL), own_ext(false),
64 cP(NULL), cR(NULL), cR_hp(NULL), cP_is_set(false),
65 Th(Operator::ANY_TYPE),
66 sequence(0), mesh_sequence(0), orders_changed(false), relaxed_hp(false)
67 { }
68
FiniteElementSpace(const FiniteElementSpace & orig,Mesh * mesh,const FiniteElementCollection * fec)69 FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
70 Mesh *mesh,
71 const FiniteElementCollection *fec)
72 {
73 mesh = mesh ? mesh : orig.mesh;
74 fec = fec ? fec : orig.fec;
75
76 NURBSExtension *NURBSext = NULL;
77 if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
78 {
79 #ifdef MFEM_USE_MPI
80 ParNURBSExtension *pNURBSext =
81 dynamic_cast<ParNURBSExtension *>(orig.NURBSext);
82 if (pNURBSext)
83 {
84 NURBSext = new ParNURBSExtension(*pNURBSext);
85 }
86 else
87 #endif
88 {
89 NURBSext = new NURBSExtension(*orig.NURBSext);
90 }
91 }
92
93 Constructor(mesh, NURBSext, fec, orig.vdim, orig.ordering);
94 }
95
CopyProlongationAndRestriction(const FiniteElementSpace & fes,const Array<int> * perm)96 void FiniteElementSpace::CopyProlongationAndRestriction(
97 const FiniteElementSpace &fes, const Array<int> *perm)
98 {
99 MFEM_VERIFY(cP == NULL, "");
100 MFEM_VERIFY(cR == NULL, "");
101
102 SparseMatrix *perm_mat = NULL, *perm_mat_tr = NULL;
103 if (perm)
104 {
105 int n = perm->Size();
106 perm_mat = new SparseMatrix(n, n);
107 for (int i=0; i<n; ++i)
108 {
109 double s;
110 int j = DecodeDof((*perm)[i], s);
111 perm_mat->Set(i, j, s);
112 }
113 perm_mat->Finalize();
114 perm_mat_tr = Transpose(*perm_mat);
115 }
116
117 if (fes.GetConformingProlongation() != NULL)
118 {
119 if (perm) { cP = Mult(*perm_mat, *fes.GetConformingProlongation()); }
120 else { cP = new SparseMatrix(*fes.GetConformingProlongation()); }
121 cP_is_set = true;
122 }
123 if (fes.GetConformingRestriction() != NULL)
124 {
125 if (perm) { cR = Mult(*fes.GetConformingRestriction(), *perm_mat_tr); }
126 else { cR = new SparseMatrix(*fes.GetConformingRestriction()); }
127 }
128
129 delete perm_mat;
130 delete perm_mat_tr;
131 }
132
SetElementOrder(int i,int p)133 void FiniteElementSpace::SetElementOrder(int i, int p)
134 {
135 MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),
136 "Space has not been Updated() after a Mesh change.");
137 MFEM_VERIFY(i >= 0 && i < GetNE(), "Invalid element index");
138 MFEM_VERIFY(p >= 0 && p <= MaxVarOrder, "Order out of range");
139 MFEM_ASSERT(!elem_order.Size() || elem_order.Size() == GetNE(),
140 "Internal error");
141
142 if (elem_order.Size()) // already a variable-order space
143 {
144 if (elem_order[i] != p)
145 {
146 elem_order[i] = p;
147 orders_changed = true;
148 }
149 }
150 else // convert space to variable-order space
151 {
152 elem_order.SetSize(GetNE());
153 elem_order = fec->GetOrder();
154
155 elem_order[i] = p;
156 orders_changed = true;
157 }
158 }
159
GetElementOrder(int i) const160 int FiniteElementSpace::GetElementOrder(int i) const
161 {
162 MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),
163 "Space has not been Updated() after a Mesh change.");
164 MFEM_VERIFY(i >= 0 && i < GetNE(), "Invalid element index");
165 MFEM_ASSERT(!elem_order.Size() || elem_order.Size() == GetNE(),
166 "Internal error");
167
168 return GetElementOrderImpl(i);
169 }
170
GetElementOrderImpl(int i) const171 int FiniteElementSpace::GetElementOrderImpl(int i) const
172 {
173 // (this is an internal version of GetElementOrder without asserts and checks)
174 return elem_order.Size() ? elem_order[i] : fec->GetOrder();
175 }
176
GetVDofs(int vd,Array<int> & dofs,int ndofs) const177 void FiniteElementSpace::GetVDofs(int vd, Array<int>& dofs, int ndofs) const
178 {
179 if (ndofs < 0) { ndofs = this->ndofs; }
180
181 if (ordering == Ordering::byNODES)
182 {
183 for (int i = 0; i < dofs.Size(); i++)
184 {
185 dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs, vdim, i, vd);
186 }
187 }
188 else
189 {
190 for (int i = 0; i < dofs.Size(); i++)
191 {
192 dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs, vdim, i, vd);
193 }
194 }
195 }
196
DofsToVDofs(Array<int> & dofs,int ndofs) const197 void FiniteElementSpace::DofsToVDofs (Array<int> &dofs, int ndofs) const
198 {
199 if (vdim == 1) { return; }
200 if (ndofs < 0) { ndofs = this->ndofs; }
201
202 if (ordering == Ordering::byNODES)
203 {
204 Ordering::DofsToVDofs<Ordering::byNODES>(ndofs, vdim, dofs);
205 }
206 else
207 {
208 Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs, vdim, dofs);
209 }
210 }
211
DofsToVDofs(int vd,Array<int> & dofs,int ndofs) const212 void FiniteElementSpace::DofsToVDofs(int vd, Array<int> &dofs, int ndofs) const
213 {
214 if (vdim == 1) { return; }
215 if (ndofs < 0) { ndofs = this->ndofs; }
216
217 if (ordering == Ordering::byNODES)
218 {
219 for (int i = 0; i < dofs.Size(); i++)
220 {
221 dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs, vdim, dofs[i], vd);
222 }
223 }
224 else
225 {
226 for (int i = 0; i < dofs.Size(); i++)
227 {
228 dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dofs[i], vd);
229 }
230 }
231 }
232
DofToVDof(int dof,int vd,int ndofs) const233 int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs) const
234 {
235 if (vdim == 1) { return dof; }
236 if (ndofs < 0) { ndofs = this->ndofs; }
237
238 if (ordering == Ordering::byNODES)
239 {
240 return Ordering::Map<Ordering::byNODES>(ndofs, vdim, dof, vd);
241 }
242 else
243 {
244 return Ordering::Map<Ordering::byVDIM>(ndofs, vdim, dof, vd);
245 }
246 }
247
248 // static function
AdjustVDofs(Array<int> & vdofs)249 void FiniteElementSpace::AdjustVDofs (Array<int> &vdofs)
250 {
251 int n = vdofs.Size(), *vdof = vdofs;
252 for (int i = 0; i < n; i++)
253 {
254 int j;
255 if ((j = vdof[i]) < 0)
256 {
257 vdof[i] = -1-j;
258 }
259 }
260 }
261
GetElementVDofs(int i,Array<int> & vdofs) const262 void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs) const
263 {
264 GetElementDofs(i, vdofs);
265 DofsToVDofs(vdofs);
266 }
267
GetBdrElementVDofs(int i,Array<int> & vdofs) const268 void FiniteElementSpace::GetBdrElementVDofs(int i, Array<int> &vdofs) const
269 {
270 GetBdrElementDofs(i, vdofs);
271 DofsToVDofs(vdofs);
272 }
273
GetFaceVDofs(int i,Array<int> & vdofs) const274 void FiniteElementSpace::GetFaceVDofs(int i, Array<int> &vdofs) const
275 {
276 GetFaceDofs(i, vdofs);
277 DofsToVDofs(vdofs);
278 }
279
GetEdgeVDofs(int i,Array<int> & vdofs) const280 void FiniteElementSpace::GetEdgeVDofs(int i, Array<int> &vdofs) const
281 {
282 GetEdgeDofs(i, vdofs);
283 DofsToVDofs(vdofs);
284 }
285
GetVertexVDofs(int i,Array<int> & vdofs) const286 void FiniteElementSpace::GetVertexVDofs(int i, Array<int> &vdofs) const
287 {
288 GetVertexDofs(i, vdofs);
289 DofsToVDofs(vdofs);
290 }
291
GetElementInteriorVDofs(int i,Array<int> & vdofs) const292 void FiniteElementSpace::GetElementInteriorVDofs(int i, Array<int> &vdofs) const
293 {
294 GetElementInteriorDofs(i, vdofs);
295 DofsToVDofs(vdofs);
296 }
297
GetEdgeInteriorVDofs(int i,Array<int> & vdofs) const298 void FiniteElementSpace::GetEdgeInteriorVDofs(int i, Array<int> &vdofs) const
299 {
300 GetEdgeInteriorDofs(i, vdofs);
301 DofsToVDofs(vdofs);
302 }
303
BuildElementToDofTable() const304 void FiniteElementSpace::BuildElementToDofTable() const
305 {
306 if (elem_dof) { return; }
307
308 // TODO: can we call GetElementDofs only once per element?
309 Table *el_dof = new Table;
310 Array<int> dofs;
311 el_dof -> MakeI (mesh -> GetNE());
312 for (int i = 0; i < mesh -> GetNE(); i++)
313 {
314 GetElementDofs (i, dofs);
315 el_dof -> AddColumnsInRow (i, dofs.Size());
316 }
317 el_dof -> MakeJ();
318 for (int i = 0; i < mesh -> GetNE(); i++)
319 {
320 GetElementDofs (i, dofs);
321 el_dof -> AddConnections (i, (int *)dofs, dofs.Size());
322 }
323 el_dof -> ShiftUpI();
324 elem_dof = el_dof;
325 }
326
BuildBdrElementToDofTable() const327 void FiniteElementSpace::BuildBdrElementToDofTable() const
328 {
329 if (bdr_elem_dof) { return; }
330
331 Table *bel_dof = new Table;
332 Array<int> dofs;
333 bel_dof->MakeI(mesh->GetNBE());
334 for (int i = 0; i < mesh->GetNBE(); i++)
335 {
336 GetBdrElementDofs(i, dofs);
337 bel_dof->AddColumnsInRow(i, dofs.Size());
338 }
339 bel_dof->MakeJ();
340 for (int i = 0; i < mesh->GetNBE(); i++)
341 {
342 GetBdrElementDofs(i, dofs);
343 bel_dof->AddConnections(i, (int *)dofs, dofs.Size());
344 }
345 bel_dof->ShiftUpI();
346 bdr_elem_dof = bel_dof;
347 }
348
BuildFaceToDofTable() const349 void FiniteElementSpace::BuildFaceToDofTable() const
350 {
351 // Here, "face" == (dim-1)-dimensional mesh entity.
352
353 if (face_dof) { return; }
354
355 if (NURBSext) { BuildNURBSFaceToDofTable(); return; }
356
357 Table *fc_dof = new Table;
358 Array<int> dofs;
359 fc_dof->MakeI(mesh->GetNumFaces());
360 for (int i = 0; i < fc_dof->Size(); i++)
361 {
362 GetFaceDofs(i, dofs, 0);
363 fc_dof->AddColumnsInRow(i, dofs.Size());
364 }
365 fc_dof->MakeJ();
366 for (int i = 0; i < fc_dof->Size(); i++)
367 {
368 GetFaceDofs(i, dofs, 0);
369 fc_dof->AddConnections(i, (int *)dofs, dofs.Size());
370 }
371 fc_dof->ShiftUpI();
372 face_dof = fc_dof;
373 }
374
RebuildElementToDofTable()375 void FiniteElementSpace::RebuildElementToDofTable()
376 {
377 delete elem_dof;
378 elem_dof = NULL;
379 BuildElementToDofTable();
380 }
381
ReorderElementToDofTable()382 void FiniteElementSpace::ReorderElementToDofTable()
383 {
384 Array<int> dof_marker(ndofs);
385
386 dof_marker = -1;
387
388 int *J = elem_dof->GetJ(), nnz = elem_dof->Size_of_connections();
389 for (int k = 0, dof_counter = 0; k < nnz; k++)
390 {
391 const int sdof = J[k]; // signed dof
392 const int dof = (sdof < 0) ? -1-sdof : sdof;
393 int new_dof = dof_marker[dof];
394 if (new_dof < 0)
395 {
396 dof_marker[dof] = new_dof = dof_counter++;
397 }
398 J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
399 }
400 }
401
BuildDofToArrays()402 void FiniteElementSpace::BuildDofToArrays()
403 {
404 if (dof_elem_array.Size()) { return; }
405
406 BuildElementToDofTable();
407
408 dof_elem_array.SetSize (ndofs);
409 dof_ldof_array.SetSize (ndofs);
410 dof_elem_array = -1;
411 for (int i = 0; i < mesh -> GetNE(); i++)
412 {
413 const int *dofs = elem_dof -> GetRow(i);
414 const int n = elem_dof -> RowSize(i);
415 for (int j = 0; j < n; j++)
416 {
417 int dof = DecodeDof(dofs[j]);
418 if (dof_elem_array[dof] < 0)
419 {
420 dof_elem_array[dof] = i;
421 dof_ldof_array[dof] = j;
422 }
423 }
424 }
425 }
426
mark_dofs(const Array<int> & dofs,Array<int> & mark_array)427 static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
428 {
429 for (int i = 0; i < dofs.Size(); i++)
430 {
431 int k = dofs[i];
432 if (k < 0) { k = -1 - k; }
433 mark_array[k] = -1;
434 }
435 }
436
GetEssentialVDofs(const Array<int> & bdr_attr_is_ess,Array<int> & ess_vdofs,int component) const437 void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
438 Array<int> &ess_vdofs,
439 int component) const
440 {
441 Array<int> vdofs, dofs;
442
443 ess_vdofs.SetSize(GetVSize());
444 ess_vdofs = 0;
445
446 for (int i = 0; i < GetNBE(); i++)
447 {
448 if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
449 {
450 if (component < 0)
451 {
452 // Mark all components.
453 GetBdrElementVDofs(i, vdofs);
454 mark_dofs(vdofs, ess_vdofs);
455 }
456 else
457 {
458 GetBdrElementDofs(i, dofs);
459 for (int d = 0; d < dofs.Size(); d++)
460 { dofs[d] = DofToVDof(dofs[d], component); }
461 mark_dofs(dofs, ess_vdofs);
462 }
463 }
464 }
465
466 // mark possible hidden boundary edges in a non-conforming mesh, also
467 // local DOFs affected by boundary elements on other processors
468 if (Nonconforming())
469 {
470 Array<int> bdr_verts, bdr_edges;
471 mesh->ncmesh->GetBoundaryClosure(bdr_attr_is_ess, bdr_verts, bdr_edges);
472
473 for (int i = 0; i < bdr_verts.Size(); i++)
474 {
475 if (component < 0)
476 {
477 GetVertexVDofs(bdr_verts[i], vdofs);
478 mark_dofs(vdofs, ess_vdofs);
479 }
480 else
481 {
482 GetVertexDofs(bdr_verts[i], dofs);
483 for (int d = 0; d < dofs.Size(); d++)
484 { dofs[d] = DofToVDof(dofs[d], component); }
485 mark_dofs(dofs, ess_vdofs);
486 }
487 }
488 for (int i = 0; i < bdr_edges.Size(); i++)
489 {
490 if (component < 0)
491 {
492 GetEdgeVDofs(bdr_edges[i], vdofs);
493 mark_dofs(vdofs, ess_vdofs);
494 }
495 else
496 {
497 GetEdgeDofs(bdr_edges[i], dofs);
498 for (int d = 0; d < dofs.Size(); d++)
499 { dofs[d] = DofToVDof(dofs[d], component); }
500 mark_dofs(dofs, ess_vdofs);
501 }
502 }
503 }
504 }
505
GetEssentialTrueDofs(const Array<int> & bdr_attr_is_ess,Array<int> & ess_tdof_list,int component)506 void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
507 Array<int> &ess_tdof_list,
508 int component)
509 {
510 Array<int> ess_vdofs, ess_tdofs;
511 GetEssentialVDofs(bdr_attr_is_ess, ess_vdofs, component);
512 const SparseMatrix *R = GetConformingRestriction();
513 if (!R)
514 {
515 ess_tdofs.MakeRef(ess_vdofs);
516 }
517 else
518 {
519 R->BooleanMult(ess_vdofs, ess_tdofs);
520 }
521 MarkerToList(ess_tdofs, ess_tdof_list);
522 }
523
GetBoundaryTrueDofs(Array<int> & boundary_dofs,int component)524 void FiniteElementSpace::GetBoundaryTrueDofs(Array<int> &boundary_dofs,
525 int component)
526 {
527 if (mesh->bdr_attributes.Size())
528 {
529 Array<int> ess_bdr(mesh->bdr_attributes.Max());
530 ess_bdr = 1;
531 GetEssentialTrueDofs(ess_bdr, boundary_dofs, component);
532 }
533 else
534 {
535 boundary_dofs.DeleteAll();
536 }
537 }
538
539 // static method
MarkerToList(const Array<int> & marker,Array<int> & list)540 void FiniteElementSpace::MarkerToList(const Array<int> &marker,
541 Array<int> &list)
542 {
543 int num_marked = 0;
544 marker.HostRead(); // make sure we can read the array on host
545 for (int i = 0; i < marker.Size(); i++)
546 {
547 if (marker[i]) { num_marked++; }
548 }
549 list.SetSize(0);
550 list.HostWrite();
551 list.Reserve(num_marked);
552 for (int i = 0; i < marker.Size(); i++)
553 {
554 if (marker[i]) { list.Append(i); }
555 }
556 }
557
558 // static method
ListToMarker(const Array<int> & list,int marker_size,Array<int> & marker,int mark_val)559 void FiniteElementSpace::ListToMarker(const Array<int> &list, int marker_size,
560 Array<int> &marker, int mark_val)
561 {
562 list.HostRead(); // make sure we can read the array on host
563 marker.SetSize(marker_size);
564 marker.HostWrite();
565 marker = 0;
566 for (int i = 0; i < list.Size(); i++)
567 {
568 marker[list[i]] = mark_val;
569 }
570 }
571
ConvertToConformingVDofs(const Array<int> & dofs,Array<int> & cdofs)572 void FiniteElementSpace::ConvertToConformingVDofs(const Array<int> &dofs,
573 Array<int> &cdofs)
574 {
575 GetConformingProlongation();
576 if (cP) { cP->BooleanMultTranspose(dofs, cdofs); }
577 else { dofs.Copy(cdofs); }
578 }
579
ConvertFromConformingVDofs(const Array<int> & cdofs,Array<int> & dofs)580 void FiniteElementSpace::ConvertFromConformingVDofs(const Array<int> &cdofs,
581 Array<int> &dofs)
582 {
583 GetConformingRestriction();
584 if (cR) { cR->BooleanMultTranspose(cdofs, dofs); }
585 else { cdofs.Copy(dofs); }
586 }
587
588 SparseMatrix *
D2C_GlobalRestrictionMatrix(FiniteElementSpace * cfes)589 FiniteElementSpace::D2C_GlobalRestrictionMatrix (FiniteElementSpace *cfes)
590 {
591 int i, j;
592 Array<int> d_vdofs, c_vdofs;
593 SparseMatrix *R;
594
595 R = new SparseMatrix (cfes -> GetVSize(), GetVSize());
596
597 for (i = 0; i < mesh -> GetNE(); i++)
598 {
599 this -> GetElementVDofs (i, d_vdofs);
600 cfes -> GetElementVDofs (i, c_vdofs);
601
602 #ifdef MFEM_DEBUG
603 if (d_vdofs.Size() != c_vdofs.Size())
604 {
605 mfem_error ("FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
606 }
607 #endif
608
609 for (j = 0; j < d_vdofs.Size(); j++)
610 {
611 R -> Set (c_vdofs[j], d_vdofs[j], 1.0);
612 }
613 }
614
615 R -> Finalize();
616
617 return R;
618 }
619
620 SparseMatrix *
D2Const_GlobalRestrictionMatrix(FiniteElementSpace * cfes)621 FiniteElementSpace::D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
622 {
623 int i, j;
624 Array<int> d_dofs, c_dofs;
625 SparseMatrix *R;
626
627 R = new SparseMatrix (cfes -> GetNDofs(), ndofs);
628
629 for (i = 0; i < mesh -> GetNE(); i++)
630 {
631 this -> GetElementDofs (i, d_dofs);
632 cfes -> GetElementDofs (i, c_dofs);
633
634 #ifdef MFEM_DEBUG
635 if (c_dofs.Size() != 1)
636 mfem_error ("FiniteElementSpace::"
637 "D2Const_GlobalRestrictionMatrix (...)");
638 #endif
639
640 for (j = 0; j < d_dofs.Size(); j++)
641 {
642 R -> Set (c_dofs[0], d_dofs[j], 1.0);
643 }
644 }
645
646 R -> Finalize();
647
648 return R;
649 }
650
651 SparseMatrix *
H2L_GlobalRestrictionMatrix(FiniteElementSpace * lfes)652 FiniteElementSpace::H2L_GlobalRestrictionMatrix (FiniteElementSpace *lfes)
653 {
654 SparseMatrix *R;
655 DenseMatrix loc_restr;
656 Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
657
658 int vdim = lfes->GetVDim();
659 R = new SparseMatrix (vdim * lfes -> GetNDofs(), vdim * ndofs);
660
661 Geometry::Type cached_geom = Geometry::INVALID;
662 const FiniteElement *h_fe = NULL;
663 const FiniteElement *l_fe = NULL;
664 IsoparametricTransformation T;
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 // Assuming 'loc_restr' depends only on the Geometry::Type.
672 const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
673 if (geom != cached_geom)
674 {
675 h_fe = this -> GetFE (i);
676 l_fe = lfes -> GetFE (i);
677 T.SetIdentityTransformation(h_fe->GetGeomType());
678 h_fe->Project(*l_fe, T, loc_restr);
679 cached_geom = geom;
680 }
681
682 for (int vd = 0; vd < vdim; vd++)
683 {
684 l_dofs.Copy(l_vdofs);
685 lfes->DofsToVDofs(vd, l_vdofs);
686
687 h_dofs.Copy(h_vdofs);
688 this->DofsToVDofs(vd, h_vdofs);
689
690 R -> SetSubMatrix (l_vdofs, h_vdofs, loc_restr, 1);
691 }
692 }
693
694 R -> Finalize();
695
696 return R;
697 }
698
699 void FiniteElementSpace
AddDependencies(SparseMatrix & deps,Array<int> & master_dofs,Array<int> & slave_dofs,DenseMatrix & I,int skipfirst)700 ::AddDependencies(SparseMatrix& deps, Array<int>& master_dofs,
701 Array<int>& slave_dofs, DenseMatrix& I, int skipfirst)
702 {
703 for (int i = skipfirst; i < slave_dofs.Size(); i++)
704 {
705 int sdof = slave_dofs[i];
706 if (!deps.RowSize(sdof)) // not processed yet
707 {
708 for (int j = 0; j < master_dofs.Size(); j++)
709 {
710 double coef = I(i, j);
711 if (std::abs(coef) > 1e-12)
712 {
713 int mdof = master_dofs[j];
714 if (mdof != sdof && mdof != (-1-sdof))
715 {
716 deps.Add(sdof, mdof, coef);
717 }
718 }
719 }
720 }
721 }
722 }
723
724 void FiniteElementSpace
AddEdgeFaceDependencies(SparseMatrix & deps,Array<int> & master_dofs,const FiniteElement * master_fe,Array<int> & slave_dofs,int slave_face,const DenseMatrix * pm) const725 ::AddEdgeFaceDependencies(SparseMatrix &deps, Array<int> &master_dofs,
726 const FiniteElement *master_fe,
727 Array<int> &slave_dofs, int slave_face,
728 const DenseMatrix *pm) const
729 {
730 // In variable-order spaces in 3D, we need to only constrain interior face
731 // DOFs (this is done one level up), since edge dependencies can be more
732 // complex and are primarily handled by edge-edge dependencies. The one
733 // exception is edges of slave faces that lie in the interior of the master
734 // face, which are not covered by edge-edge relations. This function finds
735 // such edges and makes them constrained by the master face.
736 // See also https://github.com/mfem/mfem/pull/1423#issuecomment-633916643
737
738 Array<int> V, E, Eo; // TODO: LocalArray
739 mesh->GetFaceVertices(slave_face, V);
740 mesh->GetFaceEdges(slave_face, E, Eo);
741 MFEM_ASSERT(V.Size() == E.Size(), "");
742
743 DenseMatrix I;
744 IsoparametricTransformation edge_T;
745 edge_T.SetFE(&SegmentFE);
746
747 // constrain each edge of the slave face
748 for (int i = 0; i < E.Size(); i++)
749 {
750 int a = i, b = (i+1) % V.Size();
751 if (V[a] > V[b]) { std::swap(a, b); }
752
753 DenseMatrix &edge_pm = edge_T.GetPointMat();
754 edge_pm.SetSize(2, 2);
755
756 // copy two points from the face point matrix
757 double mid[2];
758 for (int j = 0; j < 2; j++)
759 {
760 edge_pm(j, 0) = (*pm)(j, a);
761 edge_pm(j, 1) = (*pm)(j, b);
762 mid[j] = 0.5*((*pm)(j, a) + (*pm)(j, b));
763 }
764
765 // check that the edge does not coincide with the master face's edge
766 const double eps = 1e-14;
767 if (mid[0] > eps && mid[0] < 1-eps &&
768 mid[1] > eps && mid[1] < 1-eps)
769 {
770 int order = GetEdgeDofs(E[i], slave_dofs, 0);
771
772 const auto *edge_fe = fec->GetFE(Geometry::SEGMENT, order);
773 edge_fe->GetTransferMatrix(*master_fe, edge_T, I);
774
775 AddDependencies(deps, master_dofs, slave_dofs, I, 0);
776 }
777 }
778 }
779
DofFinalizable(int dof,const Array<bool> & finalized,const SparseMatrix & deps)780 bool FiniteElementSpace::DofFinalizable(int dof, const Array<bool>& finalized,
781 const SparseMatrix& deps)
782 {
783 const int* dep = deps.GetRowColumns(dof);
784 int ndep = deps.RowSize(dof);
785
786 // are all constraining DOFs finalized?
787 for (int i = 0; i < ndep; i++)
788 {
789 if (!finalized[dep[i]]) { return false; }
790 }
791 return true;
792 }
793
GetDegenerateFaceDofs(int index,Array<int> & dofs,Geometry::Type master_geom,int variant) const794 int FiniteElementSpace::GetDegenerateFaceDofs(int index, Array<int> &dofs,
795 Geometry::Type master_geom,
796 int variant) const
797 {
798 // In NC meshes with prisms/tets, a special constraint occurs where a
799 // prism/tet edge is slave to another element's face (see illustration
800 // here: https://github.com/mfem/mfem/pull/713#issuecomment-495786362)
801 // Rather than introduce a new edge-face constraint type, we handle such
802 // cases as degenerate face-face constraints, where the point-matrix
803 // rectangle has zero height. This method returns DOFs for the first edge
804 // of the rectangle, duplicated in the orthogonal direction, to resemble
805 // DOFs for a quadrilateral face. The extra DOFs are ignored by
806 // FiniteElementSpace::AddDependencies.
807
808 Array<int> edof;
809 int order = GetEdgeDofs(-1 - index, edof, variant);
810
811 int nv = fec->DofForGeometry(Geometry::POINT);
812 int ne = fec->DofForGeometry(Geometry::SEGMENT);
813 int nn = 2*nv + ne;
814
815 dofs.SetSize(nn*nn);
816 if (!dofs.Size()) { return 0; }
817
818 dofs = edof[0];
819
820 // copy first two vertex DOFs
821 for (int i = 0; i < nv; i++)
822 {
823 dofs[i] = edof[i];
824 dofs[nv+i] = edof[nv+i];
825 }
826 // copy first edge DOFs
827 int face_vert = Geometry::NumVerts[master_geom];
828 for (int i = 0; i < ne; i++)
829 {
830 dofs[face_vert*nv + i] = edof[2*nv + i];
831 }
832
833 return order;
834 }
835
GetNumBorderDofs(Geometry::Type geom,int order) const836 int FiniteElementSpace::GetNumBorderDofs(Geometry::Type geom, int order) const
837 {
838 // return the number of vertex and edge DOFs that precede inner DOFs
839 int nv = fec->GetNumDof(Geometry::POINT, order);
840 int ne = fec->GetNumDof(Geometry::SEGMENT, order);
841 return Geometry::NumVerts[geom] * (nv + ne);
842 }
843
GetEntityDofs(int entity,int index,Array<int> & dofs,Geometry::Type master_geom,int variant) const844 int FiniteElementSpace::GetEntityDofs(int entity, int index, Array<int> &dofs,
845 Geometry::Type master_geom,
846 int variant) const
847 {
848 switch (entity)
849 {
850 case 0:
851 GetVertexDofs(index, dofs);
852 return 0;
853
854 case 1:
855 return GetEdgeDofs(index, dofs, variant);
856
857 default:
858 if (index >= 0)
859 {
860 return GetFaceDofs(index, dofs, variant);
861 }
862 else
863 {
864 return GetDegenerateFaceDofs(index, dofs, master_geom, variant);
865 }
866 }
867 }
868
BuildConformingInterpolation() const869 void FiniteElementSpace::BuildConformingInterpolation() const
870 {
871 #ifdef MFEM_USE_MPI
872 MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
873 "This method should not be used with a ParFiniteElementSpace!");
874 #endif
875
876 if (cP_is_set) { return; }
877 cP_is_set = true;
878
879 Array<int> master_dofs, slave_dofs, highest_dofs;
880
881 IsoparametricTransformation T;
882 DenseMatrix I;
883
884 // For each slave DOF, the dependency matrix will contain a row that
885 // expresses the slave DOF as a linear combination of its immediate master
886 // DOFs. Rows of independent DOFs will remain empty.
887 SparseMatrix deps(ndofs);
888
889 // Inverse dependencies for the cR_hp matrix in variable order spaces:
890 // For each master edge/face with more DOF sets, the inverse dependency
891 // matrix contains a row that expresses the master true DOF (lowest order)
892 // as a linear combination of the highest order set of DOFs.
893 SparseMatrix inv_deps(ndofs);
894
895 // collect local face/edge dependencies
896 for (int entity = 2; entity >= 1; entity--)
897 {
898 const NCMesh::NCList &list = mesh->ncmesh->GetNCList(entity);
899 if (!list.masters.Size()) { continue; }
900
901 // loop through all master edges/faces, constrain their slave edges/faces
902 for (const NCMesh::Master &master : list.masters)
903 {
904 Geometry::Type master_geom = master.Geom();
905
906 int p = GetEntityDofs(entity, master.index, master_dofs, master_geom);
907 if (!master_dofs.Size()) { continue; }
908
909 const FiniteElement *master_fe = fec->GetFE(master_geom, p);
910 if (!master_fe) { continue; }
911
912 switch (master_geom)
913 {
914 case Geometry::SQUARE: T.SetFE(&QuadrilateralFE); break;
915 case Geometry::TRIANGLE: T.SetFE(&TriangleFE); break;
916 case Geometry::SEGMENT: T.SetFE(&SegmentFE); break;
917 default: MFEM_ABORT("unsupported geometry");
918 }
919
920 for (int si = master.slaves_begin; si < master.slaves_end; si++)
921 {
922 const NCMesh::Slave &slave = list.slaves[si];
923
924 int q = GetEntityDofs(entity, slave.index, slave_dofs, master_geom);
925 if (!slave_dofs.Size()) { break; }
926
927 const FiniteElement *slave_fe = fec->GetFE(slave.Geom(), q);
928 list.OrientedPointMatrix(slave, T.GetPointMat());
929 slave_fe->GetTransferMatrix(*master_fe, T, I);
930
931 // variable order spaces: face edges need to be handled separately
932 int skipfirst = 0;
933 if (IsVariableOrder() && entity == 2 && slave.index >= 0)
934 {
935 skipfirst = GetNumBorderDofs(master_geom, q);
936 }
937
938 // make each slave DOF dependent on all master DOFs
939 AddDependencies(deps, master_dofs, slave_dofs, I, skipfirst);
940
941 if (skipfirst)
942 {
943 // constrain internal edge DOFs if they were skipped
944 const auto *pm = list.point_matrices[master_geom][slave.matrix];
945 AddEdgeFaceDependencies(deps, master_dofs, master_fe,
946 slave_dofs, slave.index, pm);
947 }
948 }
949
950 // Add inverse dependencies for the cR_hp matrix; if a master has
951 // more DOF sets, the lowest order set interpolates the highest one.
952 if (IsVariableOrder())
953 {
954 int nvar = GetNVariants(entity, master.index);
955 if (nvar > 1)
956 {
957 int q = GetEntityDofs(entity, master.index, highest_dofs,
958 master_geom, nvar-1);
959 const auto *highest_fe = fec->GetFE(master_geom, q);
960
961 T.SetIdentityTransformation(master_geom);
962 master_fe->GetTransferMatrix(*highest_fe, T, I);
963
964 // add dependencies only for the inner dofs
965 int skip = GetNumBorderDofs(master_geom, p);
966 AddDependencies(inv_deps, highest_dofs, master_dofs, I, skip);
967 }
968 }
969 }
970 }
971
972 // variable order spaces: enforce minimum rule on conforming edges/faces
973 if (IsVariableOrder())
974 {
975 for (int entity = 1; entity < mesh->Dimension(); entity++)
976 {
977 const Table &ent_dofs = (entity == 1) ? var_edge_dofs : var_face_dofs;
978 int num_ent = (entity == 1) ? mesh->GetNEdges() : mesh->GetNFaces();
979 MFEM_ASSERT(ent_dofs.Size() == num_ent+1, "");
980
981 // add constraints within edges/faces holding multiple DOF sets
982 Geometry::Type last_geom = Geometry::INVALID;
983 for (int i = 0; i < num_ent; i++)
984 {
985 if (ent_dofs.RowSize(i) <= 1) { continue; }
986
987 Geometry::Type geom =
988 (entity == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
989
990 if (geom != last_geom)
991 {
992 T.SetIdentityTransformation(geom);
993 last_geom = geom;
994 }
995
996 // get lowest order variant DOFs and FE
997 int p = GetEntityDofs(entity, i, master_dofs, geom, 0);
998 const auto *master_fe = fec->GetFE(geom, p);
999
1000 // constrain all higher order DOFs: interpolate lowest order function
1001 for (int variant = 1; ; variant++)
1002 {
1003 int q = GetEntityDofs(entity, i, slave_dofs, geom, variant);
1004 if (q < 0) { break; }
1005
1006 const auto *slave_fe = fec->GetFE(geom, q);
1007 slave_fe->GetTransferMatrix(*master_fe, T, I);
1008
1009 AddDependencies(deps, master_dofs, slave_dofs, I);
1010 }
1011 }
1012 }
1013 }
1014
1015 deps.Finalize();
1016 inv_deps.Finalize();
1017
1018 // DOFs that stayed independent are true DOFs
1019 int n_true_dofs = 0;
1020 for (int i = 0; i < ndofs; i++)
1021 {
1022 if (!deps.RowSize(i)) { n_true_dofs++; }
1023 }
1024
1025 // if all dofs are true dofs leave cP and cR NULL
1026 if (n_true_dofs == ndofs)
1027 {
1028 cP = cR = cR_hp = NULL; // will be treated as identities
1029 return;
1030 }
1031
1032 // create the conforming prolongation matrix cP
1033 cP = new SparseMatrix(ndofs, n_true_dofs);
1034
1035 // create the conforming restriction matrix cR
1036 int *cR_J;
1037 {
1038 int *cR_I = Memory<int>(n_true_dofs+1);
1039 double *cR_A = Memory<double>(n_true_dofs);
1040 cR_J = Memory<int>(n_true_dofs);
1041 for (int i = 0; i < n_true_dofs; i++)
1042 {
1043 cR_I[i] = i;
1044 cR_A[i] = 1.0;
1045 }
1046 cR_I[n_true_dofs] = n_true_dofs;
1047 cR = new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, ndofs);
1048 }
1049
1050 // In var. order spaces, create the restriction matrix cR_hp which is similar
1051 // to cR, but has interpolation in the extra master edge/face DOFs.
1052 cR_hp = IsVariableOrder() ? new SparseMatrix(n_true_dofs, ndofs) : NULL;
1053
1054 Array<bool> finalized(ndofs);
1055 finalized = false;
1056
1057 Array<int> cols;
1058 Vector srow;
1059
1060 // put identity in the prolongation matrix for true DOFs, initialize cR_hp
1061 for (int i = 0, true_dof = 0; i < ndofs; i++)
1062 {
1063 if (!deps.RowSize(i)) // true dof
1064 {
1065 cP->Add(i, true_dof, 1.0);
1066 cR_J[true_dof] = i;
1067 finalized[i] = true;
1068
1069 if (cR_hp)
1070 {
1071 if (inv_deps.RowSize(i))
1072 {
1073 inv_deps.GetRow(i, cols, srow);
1074 cR_hp->AddRow(true_dof, cols, srow);
1075 }
1076 else
1077 {
1078 cR_hp->Add(true_dof, i, 1.0);
1079 }
1080 }
1081
1082 true_dof++;
1083 }
1084 }
1085
1086 // Now calculate cP rows of slave DOFs as combinations of cP rows of their
1087 // master DOFs. It is possible that some slave DOFs depend on DOFs that are
1088 // themselves slaves. Here we resolve such indirect constraints by first
1089 // calculating rows of the cP matrix for DOFs whose master DOF cP rows are
1090 // already known (in the first iteration these are the true DOFs). In the
1091 // second iteration, slaves of slaves can be 'finalized' (given a row in the
1092 // cP matrix), in the third iteration slaves of slaves of slaves, etc.
1093 bool finished;
1094 int n_finalized = n_true_dofs;
1095 do
1096 {
1097 finished = true;
1098 for (int dof = 0; dof < ndofs; dof++)
1099 {
1100 if (!finalized[dof] && DofFinalizable(dof, finalized, deps))
1101 {
1102 const int* dep_col = deps.GetRowColumns(dof);
1103 const double* dep_coef = deps.GetRowEntries(dof);
1104 int n_dep = deps.RowSize(dof);
1105
1106 for (int j = 0; j < n_dep; j++)
1107 {
1108 cP->GetRow(dep_col[j], cols, srow);
1109 srow *= dep_coef[j];
1110 cP->AddRow(dof, cols, srow);
1111 }
1112
1113 finalized[dof] = true;
1114 n_finalized++;
1115 finished = false;
1116 }
1117 }
1118 }
1119 while (!finished);
1120
1121 // If everything is consistent (mesh, face orientations, etc.), we should
1122 // be able to finalize all slave DOFs, otherwise it's a serious error.
1123 MFEM_VERIFY(n_finalized == ndofs,
1124 "Error creating cP matrix: n_finalized = "
1125 << n_finalized << ", ndofs = " << ndofs);
1126
1127 cP->Finalize();
1128 if (cR_hp) { cR_hp->Finalize(); }
1129
1130 if (vdim > 1)
1131 {
1132 MakeVDimMatrix(*cP);
1133 MakeVDimMatrix(*cR);
1134 if (cR_hp) { MakeVDimMatrix(*cR_hp); }
1135 }
1136
1137 if (Device::IsEnabled()) { cP->BuildTranspose(); }
1138 }
1139
MakeVDimMatrix(SparseMatrix & mat) const1140 void FiniteElementSpace::MakeVDimMatrix(SparseMatrix &mat) const
1141 {
1142 if (vdim == 1) { return; }
1143
1144 int height = mat.Height();
1145 int width = mat.Width();
1146
1147 SparseMatrix *vmat = new SparseMatrix(vdim*height, vdim*width);
1148
1149 Array<int> dofs, vdofs;
1150 Vector srow;
1151 for (int i = 0; i < height; i++)
1152 {
1153 mat.GetRow(i, dofs, srow);
1154 for (int vd = 0; vd < vdim; vd++)
1155 {
1156 dofs.Copy(vdofs);
1157 DofsToVDofs(vd, vdofs, width);
1158 vmat->SetRow(DofToVDof(i, vd, height), vdofs, srow);
1159 }
1160 }
1161 vmat->Finalize();
1162
1163 mat.Swap(*vmat);
1164 delete vmat;
1165 }
1166
1167
GetConformingProlongation() const1168 const SparseMatrix* FiniteElementSpace::GetConformingProlongation() const
1169 {
1170 if (Conforming()) { return NULL; }
1171 if (!cP_is_set) { BuildConformingInterpolation(); }
1172 return cP;
1173 }
1174
GetConformingRestriction() const1175 const SparseMatrix* FiniteElementSpace::GetConformingRestriction() const
1176 {
1177 if (Conforming()) { return NULL; }
1178 if (!cP_is_set) { BuildConformingInterpolation(); }
1179 return cR;
1180 }
1181
GetHpConformingRestriction() const1182 const SparseMatrix* FiniteElementSpace::GetHpConformingRestriction() const
1183 {
1184 if (Conforming()) { return NULL; }
1185 if (!cP_is_set) { BuildConformingInterpolation(); }
1186 return IsVariableOrder() ? cR_hp : cR;
1187 }
1188
GetNConformingDofs() const1189 int FiniteElementSpace::GetNConformingDofs() const
1190 {
1191 const SparseMatrix* P = GetConformingProlongation();
1192 return P ? (P->Width() / vdim) : ndofs;
1193 }
1194
GetElementRestriction(ElementDofOrdering e_ordering) const1195 const Operator *FiniteElementSpace::GetElementRestriction(
1196 ElementDofOrdering e_ordering) const
1197 {
1198 // Check if we have a discontinuous space using the FE collection:
1199 if (IsDGSpace())
1200 {
1201 // TODO: when VDIM is 1, we can return IdentityOperator.
1202 if (L2E_nat.Ptr() == NULL)
1203 {
1204 // The input L-vector layout is:
1205 // * ND x NE x VDIM, for Ordering::byNODES, or
1206 // * VDIM x ND x NE, for Ordering::byVDIM.
1207 // The output E-vector layout is: ND x VDIM x NE.
1208 L2E_nat.Reset(new L2ElementRestriction(*this));
1209 }
1210 return L2E_nat.Ptr();
1211 }
1212 if (e_ordering == ElementDofOrdering::LEXICOGRAPHIC)
1213 {
1214 if (L2E_lex.Ptr() == NULL)
1215 {
1216 L2E_lex.Reset(new ElementRestriction(*this, e_ordering));
1217 }
1218 return L2E_lex.Ptr();
1219 }
1220 // e_ordering == ElementDofOrdering::NATIVE
1221 if (L2E_nat.Ptr() == NULL)
1222 {
1223 L2E_nat.Reset(new ElementRestriction(*this, e_ordering));
1224 }
1225 return L2E_nat.Ptr();
1226 }
1227
GetFaceRestriction(ElementDofOrdering e_ordering,FaceType type,L2FaceValues mul) const1228 const FaceRestriction *FiniteElementSpace::GetFaceRestriction(
1229 ElementDofOrdering e_ordering, FaceType type, L2FaceValues mul) const
1230 {
1231 const bool is_dg_space = IsDGSpace();
1232 const L2FaceValues m = (is_dg_space && mul==L2FaceValues::DoubleValued) ?
1233 L2FaceValues::DoubleValued : L2FaceValues::SingleValued;
1234 key_face key = std::make_tuple(is_dg_space, e_ordering, type, m);
1235 auto itr = L2F.find(key);
1236 if (itr != L2F.end())
1237 {
1238 return itr->second;
1239 }
1240 else
1241 {
1242 FaceRestriction *res;
1243 if (is_dg_space)
1244 {
1245 res = new L2FaceRestriction(*this, e_ordering, type, m);
1246 }
1247 else
1248 {
1249 res = new H1FaceRestriction(*this, e_ordering, type);
1250 }
1251 L2F[key] = res;
1252 return res;
1253 }
1254 }
1255
GetQuadratureInterpolator(const IntegrationRule & ir) const1256 const QuadratureInterpolator *FiniteElementSpace::GetQuadratureInterpolator(
1257 const IntegrationRule &ir) const
1258 {
1259 for (int i = 0; i < E2Q_array.Size(); i++)
1260 {
1261 const QuadratureInterpolator *qi = E2Q_array[i];
1262 if (qi->IntRule == &ir) { return qi; }
1263 }
1264
1265 QuadratureInterpolator *qi = new QuadratureInterpolator(*this, ir);
1266 E2Q_array.Append(qi);
1267 return qi;
1268 }
1269
GetQuadratureInterpolator(const QuadratureSpace & qs) const1270 const QuadratureInterpolator *FiniteElementSpace::GetQuadratureInterpolator(
1271 const QuadratureSpace &qs) const
1272 {
1273 for (int i = 0; i < E2Q_array.Size(); i++)
1274 {
1275 const QuadratureInterpolator *qi = E2Q_array[i];
1276 if (qi->qspace == &qs) { return qi; }
1277 }
1278
1279 QuadratureInterpolator *qi = new QuadratureInterpolator(*this, qs);
1280 E2Q_array.Append(qi);
1281 return qi;
1282 }
1283
1284 const FaceQuadratureInterpolator
GetFaceQuadratureInterpolator(const IntegrationRule & ir,FaceType type) const1285 *FiniteElementSpace::GetFaceQuadratureInterpolator(
1286 const IntegrationRule &ir, FaceType type) const
1287 {
1288 if (type==FaceType::Interior)
1289 {
1290 for (int i = 0; i < E2IFQ_array.Size(); i++)
1291 {
1292 const FaceQuadratureInterpolator *qi = E2IFQ_array[i];
1293 if (qi->IntRule == &ir) { return qi; }
1294 }
1295
1296 FaceQuadratureInterpolator *qi = new FaceQuadratureInterpolator(*this, ir,
1297 type);
1298 E2IFQ_array.Append(qi);
1299 return qi;
1300 }
1301 else //Boundary
1302 {
1303 for (int i = 0; i < E2BFQ_array.Size(); i++)
1304 {
1305 const FaceQuadratureInterpolator *qi = E2BFQ_array[i];
1306 if (qi->IntRule == &ir) { return qi; }
1307 }
1308
1309 FaceQuadratureInterpolator *qi = new FaceQuadratureInterpolator(*this, ir,
1310 type);
1311 E2BFQ_array.Append(qi);
1312 return qi;
1313 }
1314 }
1315
RefinementMatrix_main(const int coarse_ndofs,const Table & coarse_elem_dof,const DenseTensor localP[]) const1316 SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
1317 const int coarse_ndofs, const Table &coarse_elem_dof,
1318 const DenseTensor localP[]) const
1319 {
1320 MFEM_VERIFY(mesh->GetLastOperation() == Mesh::REFINE, "");
1321
1322 Array<int> dofs, coarse_dofs, coarse_vdofs;
1323 Vector row;
1324
1325 Mesh::GeometryList elem_geoms(*mesh);
1326
1327 SparseMatrix *P;
1328 if (elem_geoms.Size() == 1)
1329 {
1330 const int coarse_ldof = localP[elem_geoms[0]].SizeJ();
1331 P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim, coarse_ldof);
1332 }
1333 else
1334 {
1335 P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim);
1336 }
1337
1338 Array<int> mark(P->Height());
1339 mark = 0;
1340
1341 const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
1342
1343 for (int k = 0; k < mesh->GetNE(); k++)
1344 {
1345 const Embedding &emb = rtrans.embeddings[k];
1346 const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
1347 const DenseMatrix &lP = localP[geom](emb.matrix);
1348 const int fine_ldof = localP[geom].SizeI();
1349
1350 elem_dof->GetRow(k, dofs);
1351 coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
1352
1353 for (int vd = 0; vd < vdim; vd++)
1354 {
1355 coarse_dofs.Copy(coarse_vdofs);
1356 DofsToVDofs(vd, coarse_vdofs, coarse_ndofs);
1357
1358 for (int i = 0; i < fine_ldof; i++)
1359 {
1360 int r = DofToVDof(dofs[i], vd);
1361 int m = (r >= 0) ? r : (-1 - r);
1362
1363 if (!mark[m])
1364 {
1365 lP.GetRow(i, row);
1366 P->SetRow(r, coarse_vdofs, row);
1367 mark[m] = 1;
1368 }
1369 }
1370 }
1371 }
1372
1373 MFEM_ASSERT(mark.Sum() == P->Height(), "Not all rows of P set.");
1374 if (elem_geoms.Size() != 1) { P->Finalize(); }
1375 return P;
1376 }
1377
GetLocalRefinementMatrices(Geometry::Type geom,DenseTensor & localP) const1378 void FiniteElementSpace::GetLocalRefinementMatrices(
1379 Geometry::Type geom, DenseTensor &localP) const
1380 {
1381 const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
1382
1383 const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
1384 const DenseTensor &pmats = rtrans.point_matrices[geom];
1385
1386 int nmat = pmats.SizeK();
1387 int ldof = fe->GetDof();
1388
1389 IsoparametricTransformation isotr;
1390 isotr.SetIdentityTransformation(geom);
1391
1392 // calculate local interpolation matrices for all refinement types
1393 localP.SetSize(ldof, ldof, nmat);
1394 for (int i = 0; i < nmat; i++)
1395 {
1396 isotr.SetPointMat(pmats(i));
1397 fe->GetLocalInterpolation(isotr, localP(i));
1398 }
1399 }
1400
RefinementMatrix(int old_ndofs,const Table * old_elem_dof)1401 SparseMatrix* FiniteElementSpace::RefinementMatrix(int old_ndofs,
1402 const Table* old_elem_dof)
1403 {
1404 MFEM_VERIFY(GetNE() >= old_elem_dof->Size(),
1405 "Previous mesh is not coarser.");
1406
1407 Mesh::GeometryList elem_geoms(*mesh);
1408
1409 DenseTensor localP[Geometry::NumGeom];
1410 for (int i = 0; i < elem_geoms.Size(); i++)
1411 {
1412 GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1413 }
1414
1415 return RefinementMatrix_main(old_ndofs, *old_elem_dof, localP);
1416 }
1417
RefinementOperator(const FiniteElementSpace * fespace,Table * old_elem_dof,int old_ndofs)1418 FiniteElementSpace::RefinementOperator::RefinementOperator
1419 (const FiniteElementSpace* fespace, Table* old_elem_dof, int old_ndofs)
1420 : fespace(fespace)
1421 , old_elem_dof(old_elem_dof)
1422 {
1423 MFEM_VERIFY(fespace->GetNE() >= old_elem_dof->Size(),
1424 "Previous mesh is not coarser.");
1425
1426 width = old_ndofs * fespace->GetVDim();
1427 height = fespace->GetVSize();
1428
1429 Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1430
1431 for (int i = 0; i < elem_geoms.Size(); i++)
1432 {
1433 fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
1434 }
1435 }
1436
RefinementOperator(const FiniteElementSpace * fespace,const FiniteElementSpace * coarse_fes)1437 FiniteElementSpace::RefinementOperator::RefinementOperator(
1438 const FiniteElementSpace *fespace, const FiniteElementSpace *coarse_fes)
1439 : Operator(fespace->GetVSize(), coarse_fes->GetVSize()),
1440 fespace(fespace), old_elem_dof(NULL)
1441 {
1442 Mesh::GeometryList elem_geoms(*fespace->GetMesh());
1443
1444 for (int i = 0; i < elem_geoms.Size(); i++)
1445 {
1446 fespace->GetLocalRefinementMatrices(*coarse_fes, elem_geoms[i],
1447 localP[elem_geoms[i]]);
1448 }
1449
1450 // Make a copy of the coarse elem_dof Table.
1451 old_elem_dof = new Table(coarse_fes->GetElementToDofTable());
1452 }
1453
~RefinementOperator()1454 FiniteElementSpace::RefinementOperator::~RefinementOperator()
1455 {
1456 delete old_elem_dof;
1457 }
1458
1459 void FiniteElementSpace::RefinementOperator
Mult(const Vector & x,Vector & y) const1460 ::Mult(const Vector &x, Vector &y) const
1461 {
1462 Mesh* mesh = fespace->GetMesh();
1463 const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
1464
1465 Array<int> dofs, vdofs, old_dofs, old_vdofs;
1466
1467 int vdim = fespace->GetVDim();
1468 int old_ndofs = width / vdim;
1469
1470 Vector subY, subX;
1471
1472 for (int k = 0; k < mesh->GetNE(); k++)
1473 {
1474 const Embedding &emb = rtrans.embeddings[k];
1475 const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
1476 const DenseMatrix &lP = localP[geom](emb.matrix);
1477
1478 subY.SetSize(lP.Height());
1479
1480 fespace->GetElementDofs(k, dofs);
1481 old_elem_dof->GetRow(emb.parent, old_dofs);
1482
1483 for (int vd = 0; vd < vdim; vd++)
1484 {
1485 dofs.Copy(vdofs);
1486 fespace->DofsToVDofs(vd, vdofs);
1487 old_dofs.Copy(old_vdofs);
1488 fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
1489 x.GetSubVector(old_vdofs, subX);
1490 lP.Mult(subX, subY);
1491 y.SetSubVector(vdofs, subY);
1492 }
1493 }
1494 }
1495
1496 void FiniteElementSpace::RefinementOperator
MultTranspose(const Vector & x,Vector & y) const1497 ::MultTranspose(const Vector &x, Vector &y) const
1498 {
1499 y = 0.0;
1500
1501 Mesh* mesh = fespace->GetMesh();
1502 const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
1503
1504 Array<char> processed(fespace->GetVSize());
1505 processed = 0;
1506
1507 Array<int> f_dofs, c_dofs, f_vdofs, c_vdofs;
1508
1509 int vdim = fespace->GetVDim();
1510 int old_ndofs = width / vdim;
1511
1512 Vector subY, subX;
1513
1514 for (int k = 0; k < mesh->GetNE(); k++)
1515 {
1516 const Embedding &emb = rtrans.embeddings[k];
1517 const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
1518 const DenseMatrix &lP = localP[geom](emb.matrix);
1519
1520 fespace->GetElementDofs(k, f_dofs);
1521 old_elem_dof->GetRow(emb.parent, c_dofs);
1522
1523 subY.SetSize(lP.Width());
1524
1525 for (int vd = 0; vd < vdim; vd++)
1526 {
1527 f_dofs.Copy(f_vdofs);
1528 fespace->DofsToVDofs(vd, f_vdofs);
1529 c_dofs.Copy(c_vdofs);
1530 fespace->DofsToVDofs(vd, c_vdofs, old_ndofs);
1531
1532 x.GetSubVector(f_vdofs, subX);
1533
1534 for (int p = 0; p < f_dofs.Size(); ++p)
1535 {
1536 if (processed[DecodeDof(f_dofs[p])])
1537 {
1538 subX[p] = 0.0;
1539 }
1540 }
1541
1542 lP.MultTranspose(subX, subY);
1543 y.AddElementVector(c_vdofs, subY);
1544 }
1545
1546 for (int p = 0; p < f_dofs.Size(); ++p)
1547 {
1548 processed[DecodeDof(f_dofs[p])] = 1;
1549 }
1550 }
1551 }
1552
DerefinementOperator(const FiniteElementSpace * f_fes,const FiniteElementSpace * c_fes,BilinearFormIntegrator * mass_integ)1553 FiniteElementSpace::DerefinementOperator::DerefinementOperator(
1554 const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes,
1555 BilinearFormIntegrator *mass_integ)
1556 : Operator(c_fes->GetVSize(), f_fes->GetVSize()),
1557 fine_fes(f_fes)
1558 {
1559 MFEM_VERIFY(c_fes->GetOrdering() == f_fes->GetOrdering() &&
1560 c_fes->GetVDim() == f_fes->GetVDim(),
1561 "incompatible coarse and fine FE spaces");
1562
1563 IsoparametricTransformation emb_tr;
1564 Mesh *f_mesh = f_fes->GetMesh();
1565 const CoarseFineTransformations &rtrans = f_mesh->GetRefinementTransforms();
1566
1567 Mesh::GeometryList elem_geoms(*f_mesh);
1568 DenseTensor localP[Geometry::NumGeom], localM[Geometry::NumGeom];
1569 for (int gi = 0; gi < elem_geoms.Size(); gi++)
1570 {
1571 const Geometry::Type geom = elem_geoms[gi];
1572 DenseTensor &lP = localP[geom], &lM = localM[geom];
1573 const FiniteElement *fine_fe =
1574 f_fes->fec->FiniteElementForGeometry(geom);
1575 const FiniteElement *coarse_fe =
1576 c_fes->fec->FiniteElementForGeometry(geom);
1577 const DenseTensor &pmats = rtrans.point_matrices[geom];
1578
1579 lP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), pmats.SizeK());
1580 lM.SetSize(fine_fe->GetDof(), fine_fe->GetDof(), pmats.SizeK());
1581 emb_tr.SetIdentityTransformation(geom);
1582 for (int i = 0; i < pmats.SizeK(); i++)
1583 {
1584 emb_tr.SetPointMat(pmats(i));
1585 // Get the local interpolation matrix for this refinement type
1586 fine_fe->GetTransferMatrix(*coarse_fe, emb_tr, lP(i));
1587 // Get the local mass matrix for this refinement type
1588 mass_integ->AssembleElementMatrix(*fine_fe, emb_tr, lM(i));
1589 }
1590 }
1591
1592 Table ref_type_to_matrix;
1593 rtrans.GetCoarseToFineMap(*f_mesh, coarse_to_fine, coarse_to_ref_type,
1594 ref_type_to_matrix, ref_type_to_geom);
1595 MFEM_ASSERT(coarse_to_fine.Size() == c_fes->GetNE(), "");
1596
1597 const int total_ref_types = ref_type_to_geom.Size();
1598 int num_ref_types[Geometry::NumGeom], num_fine_elems[Geometry::NumGeom];
1599 Array<int> ref_type_to_coarse_elem_offset(total_ref_types);
1600 ref_type_to_fine_elem_offset.SetSize(total_ref_types);
1601 std::fill(num_ref_types, num_ref_types+Geometry::NumGeom, 0);
1602 std::fill(num_fine_elems, num_fine_elems+Geometry::NumGeom, 0);
1603 for (int i = 0; i < total_ref_types; i++)
1604 {
1605 Geometry::Type g = ref_type_to_geom[i];
1606 ref_type_to_coarse_elem_offset[i] = num_ref_types[g];
1607 ref_type_to_fine_elem_offset[i] = num_fine_elems[g];
1608 num_ref_types[g]++;
1609 num_fine_elems[g] += ref_type_to_matrix.RowSize(i);
1610 }
1611 DenseTensor localPtMP[Geometry::NumGeom];
1612 for (int g = 0; g < Geometry::NumGeom; g++)
1613 {
1614 if (num_ref_types[g] == 0) { continue; }
1615 const int fine_dofs = localP[g].SizeI();
1616 const int coarse_dofs = localP[g].SizeJ();
1617 localPtMP[g].SetSize(coarse_dofs, coarse_dofs, num_ref_types[g]);
1618 localR[g].SetSize(coarse_dofs, fine_dofs, num_fine_elems[g]);
1619 }
1620 for (int i = 0; i < total_ref_types; i++)
1621 {
1622 Geometry::Type g = ref_type_to_geom[i];
1623 DenseMatrix &lPtMP = localPtMP[g](ref_type_to_coarse_elem_offset[i]);
1624 int lR_offset = ref_type_to_fine_elem_offset[i]; // offset in localR[g]
1625 const int *mi = ref_type_to_matrix.GetRow(i);
1626 const int nm = ref_type_to_matrix.RowSize(i);
1627 lPtMP = 0.0;
1628 for (int s = 0; s < nm; s++)
1629 {
1630 DenseMatrix &lP = localP[g](mi[s]);
1631 DenseMatrix &lM = localM[g](mi[s]);
1632 DenseMatrix &lR = localR[g](lR_offset+s);
1633 MultAtB(lP, lM, lR); // lR = lP^T lM
1634 AddMult(lR, lP, lPtMP); // lPtMP += lP^T lM lP
1635 }
1636 DenseMatrixInverse lPtMP_inv(lPtMP);
1637 for (int s = 0; s < nm; s++)
1638 {
1639 DenseMatrix &lR = localR[g](lR_offset+s);
1640 lPtMP_inv.Mult(lR); // lR <- (P^T M P)^{-1} P^T M
1641 }
1642 }
1643
1644 // Make a copy of the coarse element-to-dof Table.
1645 coarse_elem_dof = new Table(c_fes->GetElementToDofTable());
1646 }
1647
~DerefinementOperator()1648 FiniteElementSpace::DerefinementOperator::~DerefinementOperator()
1649 {
1650 delete coarse_elem_dof;
1651 }
1652
1653 void FiniteElementSpace::DerefinementOperator
Mult(const Vector & x,Vector & y) const1654 ::Mult(const Vector &x, Vector &y) const
1655 {
1656 Array<int> c_vdofs, f_vdofs;
1657 Vector loc_x, loc_y;
1658 DenseMatrix loc_x_mat, loc_y_mat;
1659 const int vdim = fine_fes->GetVDim();
1660 const int coarse_ndofs = height/vdim;
1661 for (int coarse_el = 0; coarse_el < coarse_to_fine.Size(); coarse_el++)
1662 {
1663 coarse_elem_dof->GetRow(coarse_el, c_vdofs);
1664 fine_fes->DofsToVDofs(c_vdofs, coarse_ndofs);
1665 loc_y.SetSize(c_vdofs.Size());
1666 loc_y = 0.0;
1667 loc_y_mat.UseExternalData(loc_y.GetData(), c_vdofs.Size()/vdim, vdim);
1668 const int ref_type = coarse_to_ref_type[coarse_el];
1669 const Geometry::Type geom = ref_type_to_geom[ref_type];
1670 const int *fine_elems = coarse_to_fine.GetRow(coarse_el);
1671 const int num_fine_elems = coarse_to_fine.RowSize(coarse_el);
1672 const int lR_offset = ref_type_to_fine_elem_offset[ref_type];
1673 for (int s = 0; s < num_fine_elems; s++)
1674 {
1675 const DenseMatrix &lR = localR[geom](lR_offset+s);
1676 fine_fes->GetElementVDofs(fine_elems[s], f_vdofs);
1677 x.GetSubVector(f_vdofs, loc_x);
1678 loc_x_mat.UseExternalData(loc_x.GetData(), f_vdofs.Size()/vdim, vdim);
1679 AddMult(lR, loc_x_mat, loc_y_mat);
1680 }
1681 y.SetSubVector(c_vdofs, loc_y);
1682 }
1683 }
1684
GetLocalDerefinementMatrices(Geometry::Type geom,DenseTensor & localR) const1685 void FiniteElementSpace::GetLocalDerefinementMatrices(Geometry::Type geom,
1686 DenseTensor &localR) const
1687 {
1688 const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
1689
1690 const CoarseFineTransformations &dtrans =
1691 mesh->ncmesh->GetDerefinementTransforms();
1692 const DenseTensor &pmats = dtrans.point_matrices[geom];
1693
1694 const int nmat = pmats.SizeK();
1695 const int ldof = fe->GetDof();
1696
1697 IsoparametricTransformation isotr;
1698 isotr.SetIdentityTransformation(geom);
1699
1700 // calculate local restriction matrices for all refinement types
1701 localR.SetSize(ldof, ldof, nmat);
1702 for (int i = 0; i < nmat; i++)
1703 {
1704 isotr.SetPointMat(pmats(i));
1705 fe->GetLocalRestriction(isotr, localR(i));
1706 }
1707 }
1708
DerefinementMatrix(int old_ndofs,const Table * old_elem_dof)1709 SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
1710 const Table* old_elem_dof)
1711 {
1712 MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
1713 MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
1714 MFEM_VERIFY(ndofs <= old_ndofs, "Previous space is not finer.");
1715
1716 Array<int> dofs, old_dofs, old_vdofs;
1717 Vector row;
1718
1719 Mesh::GeometryList elem_geoms(*mesh);
1720
1721 DenseTensor localR[Geometry::NumGeom];
1722 for (int i = 0; i < elem_geoms.Size(); i++)
1723 {
1724 GetLocalDerefinementMatrices(elem_geoms[i], localR[elem_geoms[i]]);
1725 }
1726
1727 SparseMatrix *R = (elem_geoms.Size() != 1)
1728 ? new SparseMatrix(ndofs*vdim, old_ndofs*vdim) // variable row size
1729 : new SparseMatrix(ndofs*vdim, old_ndofs*vdim,
1730 localR[elem_geoms[0]].SizeI());
1731
1732 Array<int> mark(R->Height());
1733 mark = 0;
1734
1735 const CoarseFineTransformations &dtrans =
1736 mesh->ncmesh->GetDerefinementTransforms();
1737
1738 MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
1739
1740 int num_marked = 0;
1741 for (int k = 0; k < dtrans.embeddings.Size(); k++)
1742 {
1743 const Embedding &emb = dtrans.embeddings[k];
1744 Geometry::Type geom = mesh->GetElementBaseGeometry(emb.parent);
1745 DenseMatrix &lR = localR[geom](emb.matrix);
1746
1747 elem_dof->GetRow(emb.parent, dofs);
1748 old_elem_dof->GetRow(k, old_dofs);
1749
1750 for (int vd = 0; vd < vdim; vd++)
1751 {
1752 old_dofs.Copy(old_vdofs);
1753 DofsToVDofs(vd, old_vdofs, old_ndofs);
1754
1755 for (int i = 0; i < lR.Height(); i++)
1756 {
1757 if (!std::isfinite(lR(i, 0))) { continue; }
1758
1759 int r = DofToVDof(dofs[i], vd);
1760 int m = (r >= 0) ? r : (-1 - r);
1761
1762 if (!mark[m])
1763 {
1764 lR.GetRow(i, row);
1765 R->SetRow(r, old_vdofs, row);
1766 mark[m] = 1;
1767 num_marked++;
1768 }
1769 }
1770 }
1771 }
1772
1773 MFEM_VERIFY(num_marked == R->Height(),
1774 "internal error: not all rows of R were set.");
1775
1776 R->Finalize(); // no-op if fixed width
1777 return R;
1778 }
1779
GetLocalRefinementMatrices(const FiniteElementSpace & coarse_fes,Geometry::Type geom,DenseTensor & localP) const1780 void FiniteElementSpace::GetLocalRefinementMatrices(
1781 const FiniteElementSpace &coarse_fes, Geometry::Type geom,
1782 DenseTensor &localP) const
1783 {
1784 // Assumptions: see the declaration of the method.
1785
1786 const FiniteElement *fine_fe = fec->FiniteElementForGeometry(geom);
1787 const FiniteElement *coarse_fe =
1788 coarse_fes.fec->FiniteElementForGeometry(geom);
1789
1790 const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
1791 const DenseTensor &pmats = rtrans.point_matrices[geom];
1792
1793 int nmat = pmats.SizeK();
1794
1795 IsoparametricTransformation isotr;
1796 isotr.SetIdentityTransformation(geom);
1797
1798 // Calculate the local interpolation matrices for all refinement types
1799 localP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), nmat);
1800 for (int i = 0; i < nmat; i++)
1801 {
1802 isotr.SetPointMat(pmats(i));
1803 fine_fe->GetTransferMatrix(*coarse_fe, isotr, localP(i));
1804 }
1805 }
1806
Constructor(Mesh * mesh,NURBSExtension * NURBSext,const FiniteElementCollection * fec,int vdim,int ordering)1807 void FiniteElementSpace::Constructor(Mesh *mesh, NURBSExtension *NURBSext,
1808 const FiniteElementCollection *fec,
1809 int vdim, int ordering)
1810 {
1811 this->mesh = mesh;
1812 this->fec = fec;
1813 this->vdim = vdim;
1814 this->ordering = (Ordering::Type) ordering;
1815
1816 elem_dof = NULL;
1817 face_dof = NULL;
1818
1819 sequence = 0;
1820 orders_changed = false;
1821 relaxed_hp = false;
1822
1823 Th.SetType(Operator::ANY_TYPE);
1824
1825 const NURBSFECollection *nurbs_fec =
1826 dynamic_cast<const NURBSFECollection *>(fec);
1827 if (nurbs_fec)
1828 {
1829 MFEM_VERIFY(mesh->NURBSext, "NURBS FE space requires a NURBS mesh.");
1830
1831 if (NURBSext == NULL)
1832 {
1833 this->NURBSext = mesh->NURBSext;
1834 own_ext = 0;
1835 }
1836 else
1837 {
1838 this->NURBSext = NURBSext;
1839 own_ext = 1;
1840 }
1841 UpdateNURBS();
1842 cP = cR = cR_hp = NULL;
1843 cP_is_set = false;
1844 }
1845 else
1846 {
1847 this->NURBSext = NULL;
1848 own_ext = 0;
1849 Construct();
1850 }
1851 BuildElementToDofTable();
1852 }
1853
StealNURBSext()1854 NURBSExtension *FiniteElementSpace::StealNURBSext()
1855 {
1856 if (NURBSext && !own_ext)
1857 {
1858 mfem_error("FiniteElementSpace::StealNURBSext");
1859 }
1860 own_ext = 0;
1861
1862 return NURBSext;
1863 }
1864
UpdateNURBS()1865 void FiniteElementSpace::UpdateNURBS()
1866 {
1867 MFEM_VERIFY(NURBSext, "NURBSExt not defined.");
1868
1869 nvdofs = 0;
1870 nedofs = 0;
1871 nfdofs = 0;
1872 nbdofs = 0;
1873 bdofs = NULL;
1874
1875 delete face_dof;
1876 face_dof = NULL;
1877 face_to_be.DeleteAll();
1878
1879 dynamic_cast<const NURBSFECollection *>(fec)->Reset();
1880
1881 ndofs = NURBSext->GetNDof();
1882 elem_dof = NURBSext->GetElementDofTable();
1883 bdr_elem_dof = NURBSext->GetBdrElementDofTable();
1884
1885 mesh_sequence = mesh->GetSequence();
1886 sequence++;
1887 }
1888
BuildNURBSFaceToDofTable() const1889 void FiniteElementSpace::BuildNURBSFaceToDofTable() const
1890 {
1891 if (face_dof) { return; }
1892
1893 const int dim = mesh->Dimension();
1894
1895 // Find bdr to face mapping
1896 face_to_be.SetSize(GetNF());
1897 face_to_be = -1;
1898 for (int b = 0; b < GetNBE(); b++)
1899 {
1900 int f = mesh->GetBdrElementEdgeIndex(b);
1901 face_to_be[f] = b;
1902 }
1903
1904 // Loop over faces in correct order, to prevent a sort
1905 // Sort will destroy orientation info in ordering of dofs
1906 Array<Connection> face_dof_list;
1907 Array<int> row;
1908 for (int f = 0; f < GetNF(); f++)
1909 {
1910 int b = face_to_be[f];
1911 if (b == -1) { continue; }
1912 // FIXME: this assumes the boundary element and the face element have the
1913 // same orientation.
1914 if (dim > 1)
1915 {
1916 const Element *fe = mesh->GetFace(f);
1917 const Element *be = mesh->GetBdrElement(b);
1918 const int nv = be->GetNVertices();
1919 const int *fv = fe->GetVertices();
1920 const int *bv = be->GetVertices();
1921 for (int i = 0; i < nv; i++)
1922 {
1923 MFEM_VERIFY(fv[i] == bv[i],
1924 "non-matching face and boundary elements detected!");
1925 }
1926 }
1927 GetBdrElementDofs(b, row);
1928 Connection conn(f,0);
1929 for (int i = 0; i < row.Size(); i++)
1930 {
1931 conn.to = row[i];
1932 face_dof_list.Append(conn);
1933 }
1934 }
1935 face_dof = new Table(GetNF(), face_dof_list);
1936 }
1937
Construct()1938 void FiniteElementSpace::Construct()
1939 {
1940 // This method should be used only for non-NURBS spaces.
1941 MFEM_VERIFY(!NURBSext, "internal error");
1942
1943 // Variable order space needs a nontrivial P matrix + also ghost elements
1944 // in parallel, we thus require the mesh to be NC.
1945 MFEM_VERIFY(!IsVariableOrder() || Nonconforming(),
1946 "Variable order space requires a nonconforming mesh.");
1947
1948 elem_dof = NULL;
1949 bdr_elem_dof = NULL;
1950 face_dof = NULL;
1951
1952 ndofs = 0;
1953 nvdofs = nedofs = nfdofs = nbdofs = 0;
1954 bdofs = NULL;
1955
1956 cP = NULL;
1957 cR = NULL;
1958 cR_hp = NULL;
1959 cP_is_set = false;
1960 // 'Th' is initialized/destroyed before this method is called.
1961
1962 int dim = mesh->Dimension();
1963 int order = fec->GetOrder();
1964
1965 MFEM_VERIFY((mesh->GetNumGeometries(dim) > 0) || (mesh->GetNE() == 0),
1966 "Mesh was not correctly finalized.");
1967
1968 bool mixed_elements = (mesh->GetNumGeometries(dim) > 1);
1969 bool mixed_faces = (dim > 2 && mesh->GetNumGeometries(2) > 1);
1970
1971 Array<VarOrderBits> edge_orders, face_orders;
1972 if (IsVariableOrder())
1973 {
1974 // for variable order spaces, calculate orders of edges and faces
1975 CalcEdgeFaceVarOrders(edge_orders, face_orders);
1976 }
1977 else if (mixed_faces)
1978 {
1979 // for mixed faces we also create the var_face_dofs table, see below
1980 face_orders.SetSize(mesh->GetNFaces());
1981 face_orders = (VarOrderBits(1) << order);
1982 }
1983
1984 // assign vertex DOFs
1985 if (mesh->GetNV())
1986 {
1987 nvdofs = mesh->GetNV() * fec->GetNumDof(Geometry::POINT, order);
1988 }
1989
1990 // assign edge DOFs
1991 if (mesh->GetNEdges())
1992 {
1993 if (IsVariableOrder())
1994 {
1995 nedofs = MakeDofTable(1, edge_orders, var_edge_dofs, &var_edge_orders);
1996 }
1997 else
1998 {
1999 // the simple case: all edges are of the same order
2000 nedofs = mesh->GetNEdges() * fec->GetNumDof(Geometry::SEGMENT, order);
2001 }
2002 }
2003
2004 // assign face DOFs
2005 if (mesh->GetNFaces())
2006 {
2007 if (IsVariableOrder() || mixed_faces)
2008 {
2009 // NOTE: for simplicity, we also use Table var_face_dofs for mixed faces
2010 nfdofs = MakeDofTable(2, face_orders, var_face_dofs,
2011 IsVariableOrder() ? &var_face_orders : NULL);
2012 uni_fdof = -1;
2013 }
2014 else
2015 {
2016 // the simple case: all faces are of the same geometry and order
2017 uni_fdof = fec->GetNumDof(mesh->GetFaceGeometry(0), order);
2018 nfdofs = mesh->GetNFaces() * uni_fdof;
2019 }
2020 }
2021
2022 // assign internal ("bubble") DOFs
2023 if (mesh->GetNE() && dim > 0)
2024 {
2025 if (IsVariableOrder() || mixed_elements)
2026 {
2027 bdofs = new int[mesh->GetNE()+1];
2028 bdofs[0] = 0;
2029 for (int i = 0; i < mesh->GetNE(); i++)
2030 {
2031 int p = GetElementOrderImpl(i);
2032 nbdofs += fec->GetNumDof(mesh->GetElementGeometry(i), p);
2033 bdofs[i+1] = nbdofs;
2034 }
2035 }
2036 else
2037 {
2038 // the simple case: all elements are the same
2039 bdofs = NULL;
2040 Geometry::Type geom = mesh->GetElementGeometry(0);
2041 nbdofs = mesh->GetNE() * fec->GetNumDof(geom, order);
2042 }
2043 }
2044
2045 ndofs = nvdofs + nedofs + nfdofs + nbdofs;
2046
2047 // record the current mesh sequence number to detect refinement etc.
2048 mesh_sequence = mesh->GetSequence();
2049
2050 // increment our sequence number to let GridFunctions know they need updating
2051 sequence++;
2052
2053 // DOFs are now assigned according to current element orders
2054 orders_changed = false;
2055
2056 // Do not build elem_dof Table here: in parallel it has to be constructed
2057 // later.
2058 }
2059
MinOrder(VarOrderBits bits)2060 int FiniteElementSpace::MinOrder(VarOrderBits bits)
2061 {
2062 MFEM_ASSERT(bits != 0, "invalid bit mask");
2063 for (int order = 0; bits != 0; order++, bits >>= 1)
2064 {
2065 if (bits & 1) { return order; }
2066 }
2067 return 0;
2068 }
2069
2070 void FiniteElementSpace
CalcEdgeFaceVarOrders(Array<VarOrderBits> & edge_orders,Array<VarOrderBits> & face_orders) const2071 ::CalcEdgeFaceVarOrders(Array<VarOrderBits> &edge_orders,
2072 Array<VarOrderBits> &face_orders) const
2073 {
2074 MFEM_ASSERT(IsVariableOrder(), "");
2075 MFEM_ASSERT(Nonconforming(), "");
2076 MFEM_ASSERT(elem_order.Size() == mesh->GetNE(), "");
2077
2078 edge_orders.SetSize(mesh->GetNEdges()); edge_orders = 0;
2079 face_orders.SetSize(mesh->GetNFaces()); face_orders = 0;
2080
2081 // Calculate initial edge/face orders, as required by incident elements.
2082 // For each edge/face we accumulate in a bit-mask the orders of elements
2083 // sharing the edge/face.
2084 Array<int> E, F, ori;
2085 for (int i = 0; i < mesh->GetNE(); i++)
2086 {
2087 int order = elem_order[i];
2088 MFEM_ASSERT(order <= MaxVarOrder, "");
2089 VarOrderBits mask = (VarOrderBits(1) << order);
2090
2091 mesh->GetElementEdges(i, E, ori);
2092 for (int j = 0; j < E.Size(); j++)
2093 {
2094 edge_orders[E[j]] |= mask;
2095 }
2096
2097 if (mesh->Dimension() > 2)
2098 {
2099 mesh->GetElementFaces(i, F, ori);
2100 for (int j = 0; j < F.Size(); j++)
2101 {
2102 face_orders[F[j]] |= mask;
2103 }
2104 }
2105 }
2106
2107 if (relaxed_hp)
2108 {
2109 // for relaxed conformity we don't need the masters to match the minimum
2110 // orders of the slaves, we can stop now
2111 return;
2112 }
2113
2114 // Iterate while minimum orders propagate by master/slave relations
2115 // (and new orders also propagate from faces to incident edges).
2116 // See https://github.com/mfem/mfem/pull/1423#issuecomment-638930559
2117 // for an illustration of why this is necessary in hp meshes.
2118 bool done;
2119 do
2120 {
2121 done = true;
2122
2123 // propagate from slave edges to master edges
2124 const NCMesh::NCList &edge_list = mesh->ncmesh->GetEdgeList();
2125 for (const NCMesh::Master &master : edge_list.masters)
2126 {
2127 VarOrderBits slave_orders = 0;
2128 for (int i = master.slaves_begin; i < master.slaves_end; i++)
2129 {
2130 slave_orders |= edge_orders[edge_list.slaves[i].index];
2131 }
2132
2133 int min_order = MinOrder(slave_orders);
2134 if (min_order < MinOrder(edge_orders[master.index]))
2135 {
2136 edge_orders[master.index] |= (VarOrderBits(1) << min_order);
2137 done = false;
2138 }
2139 }
2140
2141 // propagate from slave faces(+edges) to master faces(+edges)
2142 const NCMesh::NCList &face_list = mesh->ncmesh->GetFaceList();
2143 for (const NCMesh::Master &master : face_list.masters)
2144 {
2145 VarOrderBits slave_orders = 0;
2146 for (int i = master.slaves_begin; i < master.slaves_end; i++)
2147 {
2148 const NCMesh::Slave &slave = face_list.slaves[i];
2149 if (slave.index >= 0)
2150 {
2151 slave_orders |= face_orders[slave.index];
2152
2153 mesh->GetFaceEdges(slave.index, E, ori);
2154 for (int j = 0; j < E.Size(); j++)
2155 {
2156 slave_orders |= edge_orders[E[j]];
2157 }
2158 }
2159 else
2160 {
2161 // degenerate face (i.e., edge-face constraint)
2162 slave_orders |= edge_orders[-1 - slave.index];
2163 }
2164 }
2165
2166 int min_order = MinOrder(slave_orders);
2167 if (min_order < MinOrder(face_orders[master.index]))
2168 {
2169 face_orders[master.index] |= (VarOrderBits(1) << min_order);
2170 done = false;
2171 }
2172 }
2173
2174 // make sure edges support (new) orders required by incident faces
2175 for (int i = 0; i < mesh->GetNFaces(); i++)
2176 {
2177 mesh->GetFaceEdges(i, E, ori);
2178 for (int j = 0; j < E.Size(); j++)
2179 {
2180 edge_orders[E[j]] |= face_orders[i];
2181 }
2182 }
2183 }
2184 while (!done);
2185 }
2186
MakeDofTable(int ent_dim,const Array<int> & entity_orders,Table & entity_dofs,Array<char> * var_ent_order)2187 int FiniteElementSpace::MakeDofTable(int ent_dim,
2188 const Array<int> &entity_orders,
2189 Table &entity_dofs,
2190 Array<char> *var_ent_order)
2191 {
2192 // The tables var_edge_dofs and var_face_dofs hold DOF assignments for edges
2193 // and faces of a variable order space, in which each edge/face may host
2194 // several DOF sets, called DOF set variants. Example: an edge 'i' shared by
2195 // 4 hexes of orders 2, 3, 4, 5 will hold four DOF sets, each starting at
2196 // indices e.g. 100, 101, 103, 106, respectively. These numbers are stored
2197 // in row 'i' of var_edge_dofs. Variant zero is always the lowest order DOF
2198 // set, followed by consecutive ranges of higher order DOFs. Variable order
2199 // faces are handled similarly by var_face_dofs, which holds at most two DOF
2200 // set variants per face. The tables are empty for constant-order spaces.
2201
2202 int num_ent = entity_orders.Size();
2203 int total_dofs = 0;
2204
2205 Array<Connection> list;
2206 list.Reserve(2*num_ent);
2207
2208 if (var_ent_order)
2209 {
2210 var_ent_order->SetSize(0);
2211 var_ent_order->Reserve(num_ent);
2212 }
2213
2214 // assign DOFs according to order bit masks
2215 for (int i = 0; i < num_ent; i++)
2216 {
2217 auto geom = (ent_dim == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
2218
2219 VarOrderBits bits = entity_orders[i];
2220 for (int order = 0; bits != 0; order++, bits >>= 1)
2221 {
2222 if (bits & 1)
2223 {
2224 int dofs = fec->GetNumDof(geom, order);
2225 list.Append(Connection(i, total_dofs));
2226 total_dofs += dofs;
2227
2228 if (var_ent_order) { var_ent_order->Append(order); }
2229 }
2230 }
2231 }
2232
2233 // append a dummy row as terminator
2234 list.Append(Connection(num_ent, total_dofs));
2235
2236 // build the table
2237 entity_dofs.MakeFromList(num_ent+1, list);
2238
2239 return total_dofs;
2240 }
2241
FindDofs(const Table & var_dof_table,int row,int ndof) const2242 int FiniteElementSpace::FindDofs(const Table &var_dof_table,
2243 int row, int ndof) const
2244 {
2245 const int *beg = var_dof_table.GetRow(row);
2246 const int *end = var_dof_table.GetRow(row + 1); // terminator, see above
2247
2248 while (beg < end)
2249 {
2250 // return the appropriate range of DOFs
2251 if ((beg[1] - beg[0]) == ndof) { return beg[0]; }
2252 beg++;
2253 }
2254
2255 MFEM_ABORT("DOFs not found for ndof = " << ndof);
2256 return 0;
2257 }
2258
GetEdgeOrder(int edge,int variant) const2259 int FiniteElementSpace::GetEdgeOrder(int edge, int variant) const
2260 {
2261 if (!IsVariableOrder()) { return fec->GetOrder(); }
2262
2263 const int* beg = var_edge_dofs.GetRow(edge);
2264 const int* end = var_edge_dofs.GetRow(edge + 1);
2265 if (variant >= end - beg) { return -1; } // past last variant
2266
2267 return var_edge_orders[var_edge_dofs.GetI()[edge] + variant];
2268 }
2269
GetFaceOrder(int face,int variant) const2270 int FiniteElementSpace::GetFaceOrder(int face, int variant) const
2271 {
2272 if (!IsVariableOrder())
2273 {
2274 // face order can be different from fec->GetOrder()
2275 Geometry::Type geom = mesh->GetFaceGeometry(face);
2276 return fec->FiniteElementForGeometry(geom)->GetOrder();
2277 }
2278
2279 const int* beg = var_face_dofs.GetRow(face);
2280 const int* end = var_face_dofs.GetRow(face + 1);
2281 if (variant >= end - beg) { return -1; } // past last variant
2282
2283 return var_face_orders[var_face_dofs.GetI()[face] + variant];
2284 }
2285
GetNVariants(int entity,int index) const2286 int FiniteElementSpace::GetNVariants(int entity, int index) const
2287 {
2288 MFEM_ASSERT(IsVariableOrder(), "");
2289 const Table &dof_table = (entity == 1) ? var_edge_dofs : var_face_dofs;
2290
2291 MFEM_ASSERT(index >= 0 && index < dof_table.Size(), "");
2292 return dof_table.GetRow(index + 1) - dof_table.GetRow(index);
2293 }
2294
2295 static const char* msg_orders_changed =
2296 "Element orders changed, you need to Update() the space first.";
2297
GetElementDofs(int elem,Array<int> & dofs) const2298 void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
2299 {
2300 MFEM_VERIFY(!orders_changed, msg_orders_changed);
2301
2302 if (elem_dof)
2303 {
2304 elem_dof->GetRow(elem, dofs);
2305 return;
2306 }
2307
2308 Array<int> V, E, Eo, F, Fo; // TODO: LocalArray
2309
2310 int dim = mesh->Dimension();
2311 auto geom = mesh->GetElementGeometry(elem);
2312 int order = GetElementOrderImpl(elem);
2313
2314 int nv = fec->GetNumDof(Geometry::POINT, order);
2315 int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
2316 int nb = (dim > 0) ? fec->GetNumDof(geom, order) : 0;
2317
2318 if (nv) { mesh->GetElementVertices(elem, V); }
2319 if (ne) { mesh->GetElementEdges(elem, E, Eo); }
2320
2321 int nfd = 0;
2322 if (dim > 2 && fec->HasFaceDofs(geom, order))
2323 {
2324 mesh->GetElementFaces(elem, F, Fo);
2325 for (int i = 0; i < F.Size(); i++)
2326 {
2327 nfd += fec->GetNumDof(mesh->GetFaceGeometry(F[i]), order);
2328 }
2329 }
2330
2331 dofs.SetSize(0);
2332 dofs.Reserve(nv*V.Size() + ne*E.Size() + nfd + nb);
2333
2334 if (nv) // vertex DOFs
2335 {
2336 for (int i = 0; i < V.Size(); i++)
2337 {
2338 for (int j = 0; j < nv; j++)
2339 {
2340 dofs.Append(V[i]*nv + j);
2341 }
2342 }
2343 }
2344
2345 if (ne) // edge DOFs
2346 {
2347 for (int i = 0; i < E.Size(); i++)
2348 {
2349 int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
2350 const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
2351
2352 for (int j = 0; j < ne; j++)
2353 {
2354 dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
2355 }
2356 }
2357 }
2358
2359 if (nfd) // face DOFs
2360 {
2361 for (int i = 0; i < F.Size(); i++)
2362 {
2363 auto fgeom = mesh->GetFaceGeometry(F[i]);
2364 int nf = fec->GetNumDof(fgeom, order);
2365
2366 int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F[i], nf) : F[i]*nf;
2367 const int *ind = fec->GetDofOrdering(fgeom, order, Fo[i]);
2368
2369 for (int j = 0; j < nf; j++)
2370 {
2371 dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
2372 }
2373 }
2374 }
2375
2376 if (nb) // interior ("bubble") DOFs
2377 {
2378 int bbase = bdofs ? bdofs[elem] : elem*nb;
2379 bbase += nvdofs + nedofs + nfdofs;
2380
2381 for (int j = 0; j < nb; j++)
2382 {
2383 dofs.Append(bbase + j);
2384 }
2385 }
2386 }
2387
GetFE(int i) const2388 const FiniteElement *FiniteElementSpace::GetFE(int i) const
2389 {
2390 if (i < 0 || !mesh->GetNE()) { return NULL; }
2391 MFEM_VERIFY(i < mesh->GetNE(),
2392 "Invalid element id " << i << ", maximum allowed " << mesh->GetNE()-1);
2393
2394 const FiniteElement *FE =
2395 fec->GetFE(mesh->GetElementGeometry(i), GetElementOrderImpl(i));
2396
2397 if (NURBSext)
2398 {
2399 NURBSext->LoadFE(i, FE);
2400 }
2401 else
2402 {
2403 #ifdef MFEM_DEBUG
2404 // consistency check: fec->GetOrder() and FE->GetOrder() should return
2405 // the same value (for standard, constant-order spaces)
2406 if (!IsVariableOrder() && FE->GetDim() > 0)
2407 {
2408 MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
2409 "internal error: " <<
2410 FE->GetOrder() << " != " << fec->GetOrder());
2411 }
2412 #endif
2413 }
2414
2415 return FE;
2416 }
2417
GetBdrElementDofs(int bel,Array<int> & dofs) const2418 void FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs) const
2419 {
2420 MFEM_VERIFY(!orders_changed, msg_orders_changed);
2421
2422 if (bdr_elem_dof)
2423 {
2424 bdr_elem_dof->GetRow(bel, dofs);
2425 return;
2426 }
2427
2428 Array<int> V, E, Eo; // TODO: LocalArray
2429 int F, Fo;
2430
2431 int dim = mesh->Dimension();
2432 auto geom = mesh->GetBdrElementGeometry(bel);
2433 int order = fec->GetOrder();
2434
2435 if (IsVariableOrder()) // determine order from adjacent element
2436 {
2437 int elem, info;
2438 mesh->GetBdrElementAdjacentElement(bel, elem, info);
2439 order = elem_order[elem];
2440 }
2441
2442 int nv = fec->GetNumDof(Geometry::POINT, order);
2443 int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
2444 int nf = (dim > 2) ? fec->GetNumDof(geom, order) : 0;
2445
2446 if (nv) { mesh->GetBdrElementVertices(bel, V); }
2447 if (ne) { mesh->GetBdrElementEdges(bel, E, Eo); }
2448 if (nf) { mesh->GetBdrElementFace(bel, &F, &Fo); }
2449
2450 dofs.SetSize(0);
2451 dofs.Reserve(nv*V.Size() + ne*E.Size() + nf);
2452
2453 if (nv) // vertex DOFs
2454 {
2455 for (int i = 0; i < V.Size(); i++)
2456 {
2457 for (int j = 0; j < nv; j++)
2458 {
2459 dofs.Append(V[i]*nv + j);
2460 }
2461 }
2462 }
2463
2464 if (ne) // edge DOFs
2465 {
2466 for (int i = 0; i < E.Size(); i++)
2467 {
2468 int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
2469 const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
2470
2471 for (int j = 0; j < ne; j++)
2472 {
2473 dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
2474 }
2475 }
2476 }
2477
2478 if (nf) // face DOFs
2479 {
2480 int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F, nf) : F*nf;
2481 const int *ind = fec->GetDofOrdering(geom, order, Fo);
2482
2483 for (int j = 0; j < nf; j++)
2484 {
2485 dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
2486 }
2487 }
2488 }
2489
GetFaceDofs(int face,Array<int> & dofs,int variant) const2490 int FiniteElementSpace::GetFaceDofs(int face, Array<int> &dofs,
2491 int variant) const
2492 {
2493 MFEM_VERIFY(!orders_changed, msg_orders_changed);
2494
2495 // If face_dof is already built, use it.
2496 // If it is not and we have a NURBS space, build the face_dof and use it.
2497 if ((face_dof && variant == 0) ||
2498 (NURBSext && (BuildNURBSFaceToDofTable(), true)))
2499 {
2500 face_dof->GetRow(face, dofs);
2501 return fec->GetOrder();
2502 }
2503
2504 int order, nf, fbase;
2505 int dim = mesh->Dimension();
2506 auto fgeom = (dim > 2) ? mesh->GetFaceGeometry(face) : Geometry::INVALID;
2507
2508 if (var_face_dofs.Size() > 0) // variable orders or *mixed* faces
2509 {
2510 const int* beg = var_face_dofs.GetRow(face);
2511 const int* end = var_face_dofs.GetRow(face + 1);
2512 if (variant >= end - beg) { return -1; } // past last face DOFs
2513
2514 fbase = beg[variant];
2515 nf = beg[variant+1] - fbase;
2516
2517 order = !IsVariableOrder() ? fec->GetOrder() :
2518 var_face_orders[var_face_dofs.GetI()[face] + variant];
2519 MFEM_ASSERT(fec->GetNumDof(fgeom, order) == nf, "");
2520 }
2521 else
2522 {
2523 if (variant > 0) { return -1; }
2524 order = fec->GetOrder();
2525 nf = (dim > 2) ? fec->GetNumDof(fgeom, order) : 0;
2526 fbase = face*nf;
2527 }
2528
2529 // for 1D, 2D and 3D faces
2530 int nv = fec->GetNumDof(Geometry::POINT, order);
2531 int ne = (dim > 1) ? fec->GetNumDof(Geometry::SEGMENT, order) : 0;
2532
2533 Array<int> V, E, Eo;
2534 if (nv) { mesh->GetFaceVertices(face, V); }
2535 if (ne) { mesh->GetFaceEdges(face, E, Eo); }
2536
2537 dofs.SetSize(0);
2538 dofs.Reserve(V.Size() * nv + E.Size() * ne + nf);
2539
2540 if (nv) // vertex DOFs
2541 {
2542 for (int i = 0; i < V.Size(); i++)
2543 {
2544 for (int j = 0; j < nv; j++)
2545 {
2546 dofs.Append(V[i]*nv + j);
2547 }
2548 }
2549 }
2550 if (ne) // edge DOFs
2551 {
2552 for (int i = 0; i < E.Size(); i++)
2553 {
2554 int ebase = IsVariableOrder() ? FindEdgeDof(E[i], ne) : E[i]*ne;
2555 const int *ind = fec->GetDofOrdering(Geometry::SEGMENT, order, Eo[i]);
2556
2557 for (int j = 0; j < ne; j++)
2558 {
2559 dofs.Append(EncodeDof(nvdofs + ebase, ind[j]));
2560 }
2561 }
2562 }
2563 for (int j = 0; j < nf; j++)
2564 {
2565 dofs.Append(nvdofs + nedofs + fbase + j);
2566 }
2567
2568 return order;
2569 }
2570
GetEdgeDofs(int edge,Array<int> & dofs,int variant) const2571 int FiniteElementSpace::GetEdgeDofs(int edge, Array<int> &dofs,
2572 int variant) const
2573 {
2574 MFEM_VERIFY(!orders_changed, msg_orders_changed);
2575
2576 int order, ne, base;
2577 if (IsVariableOrder())
2578 {
2579 const int* beg = var_edge_dofs.GetRow(edge);
2580 const int* end = var_edge_dofs.GetRow(edge + 1);
2581 if (variant >= end - beg) { return -1; } // past last edge DOFs
2582
2583 base = beg[variant];
2584 ne = beg[variant+1] - base;
2585
2586 order = var_edge_orders[var_edge_dofs.GetI()[edge] + variant];
2587 MFEM_ASSERT(fec->GetNumDof(Geometry::SEGMENT, order) == ne, "");
2588 }
2589 else
2590 {
2591 if (variant > 0) { return -1; }
2592 order = fec->GetOrder();
2593 ne = fec->GetNumDof(Geometry::SEGMENT, order);
2594 base = edge*ne;
2595 }
2596
2597 Array<int> V; // TODO: LocalArray
2598 int nv = fec->GetNumDof(Geometry::POINT, order);
2599 if (nv) { mesh->GetEdgeVertices(edge, V); }
2600
2601 dofs.SetSize(0);
2602 dofs.Reserve(2*nv + ne);
2603
2604 for (int i = 0; i < 2; i++)
2605 {
2606 for (int j = 0; j < nv; j++)
2607 {
2608 dofs.Append(V[i]*nv + j);
2609 }
2610 }
2611 for (int j = 0; j < ne; j++)
2612 {
2613 dofs.Append(nvdofs + base + j);
2614 }
2615
2616 return order;
2617 }
2618
GetVertexDofs(int i,Array<int> & dofs) const2619 void FiniteElementSpace::GetVertexDofs(int i, Array<int> &dofs) const
2620 {
2621 int nv = fec->DofForGeometry(Geometry::POINT);
2622 dofs.SetSize(nv);
2623 for (int j = 0; j < nv; j++)
2624 {
2625 dofs[j] = i*nv+j;
2626 }
2627 }
2628
GetElementInteriorDofs(int i,Array<int> & dofs) const2629 void FiniteElementSpace::GetElementInteriorDofs(int i, Array<int> &dofs) const
2630 {
2631 MFEM_VERIFY(!orders_changed, msg_orders_changed);
2632
2633 int nb = fec->GetNumDof(mesh->GetElementGeometry(i), GetElementOrderImpl(i));
2634 int base = bdofs ? bdofs[i] : i*nb;
2635
2636 dofs.SetSize(nb);
2637 base += nvdofs + nedofs + nfdofs;
2638 for (int j = 0; j < nb; j++)
2639 {
2640 dofs[j] = base + j;
2641 }
2642 }
2643
GetNumElementInteriorDofs(int i) const2644 int FiniteElementSpace::GetNumElementInteriorDofs(int i) const
2645 {
2646 return fec->GetNumDof(mesh->GetElementGeometry(i),
2647 GetElementOrderImpl(i));
2648 }
2649
GetEdgeInteriorDofs(int i,Array<int> & dofs) const2650 void FiniteElementSpace::GetEdgeInteriorDofs(int i, Array<int> &dofs) const
2651 {
2652 MFEM_VERIFY(!IsVariableOrder(), "not implemented");
2653
2654 int ne = fec->DofForGeometry(Geometry::SEGMENT);
2655 dofs.SetSize (ne);
2656 for (int j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
2657 {
2658 dofs[j] = k;
2659 }
2660 }
2661
GetFaceInteriorDofs(int i,Array<int> & dofs) const2662 void FiniteElementSpace::GetFaceInteriorDofs(int i, Array<int> &dofs) const
2663 {
2664 MFEM_VERIFY(!IsVariableOrder(), "not implemented");
2665
2666 int nf, base;
2667 if (var_face_dofs.Size() > 0) // mixed faces
2668 {
2669 base = var_face_dofs.GetRow(i)[0];
2670 nf = var_face_dofs.GetRow(i)[1] - base;
2671 }
2672 else
2673 {
2674 auto geom = mesh->GetFaceGeometry(0);
2675 nf = fec->GetNumDof(geom, fec->GetOrder());
2676 base = i*nf;
2677 }
2678
2679 dofs.SetSize(nf);
2680 for (int j = 0; j < nf; j++)
2681 {
2682 dofs[j] = nvdofs + nedofs + base + j;
2683 }
2684 }
2685
GetBE(int i) const2686 const FiniteElement *FiniteElementSpace::GetBE(int i) const
2687 {
2688 int order = fec->GetOrder();
2689
2690 if (IsVariableOrder()) // determine order from adjacent element
2691 {
2692 int elem, info;
2693 mesh->GetBdrElementAdjacentElement(i, elem, info);
2694 order = elem_order[elem];
2695 }
2696
2697 const FiniteElement *BE;
2698 switch (mesh->Dimension())
2699 {
2700 case 1:
2701 BE = fec->GetFE(Geometry::POINT, order);
2702 break;
2703 case 2:
2704 BE = fec->GetFE(Geometry::SEGMENT, order);
2705 break;
2706 case 3:
2707 default:
2708 BE = fec->GetFE(mesh->GetBdrElementBaseGeometry(i), order);
2709 }
2710
2711 if (NURBSext)
2712 {
2713 NURBSext->LoadBE(i, BE);
2714 }
2715
2716 return BE;
2717 }
2718
GetFaceElement(int i) const2719 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
2720 {
2721 MFEM_VERIFY(!IsVariableOrder(), "not implemented");
2722
2723 const FiniteElement *fe;
2724 switch (mesh->Dimension())
2725 {
2726 case 1:
2727 fe = fec->FiniteElementForGeometry(Geometry::POINT);
2728 break;
2729 case 2:
2730 fe = fec->FiniteElementForGeometry(Geometry::SEGMENT);
2731 break;
2732 case 3:
2733 default:
2734 fe = fec->FiniteElementForGeometry(mesh->GetFaceBaseGeometry(i));
2735 }
2736
2737 if (NURBSext)
2738 {
2739 // Ensure 'face_to_be' is built:
2740 if (!face_dof) { BuildNURBSFaceToDofTable(); }
2741 MFEM_ASSERT(face_to_be[i] >= 0,
2742 "NURBS mesh: only boundary faces are supported!");
2743 NURBSext->LoadBE(face_to_be[i], fe);
2744 }
2745
2746 return fe;
2747 }
2748
GetEdgeElement(int i,int variant) const2749 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
2750 int variant) const
2751 {
2752 MFEM_ASSERT(mesh->Dimension() > 1, "No edges with mesh dimension < 2");
2753
2754 int eo = IsVariableOrder() ? GetEdgeOrder(i, variant) : fec->GetOrder();
2755 return fec->GetFE(Geometry::SEGMENT, eo);
2756 }
2757
2758 const FiniteElement *FiniteElementSpace
GetTraceElement(int i,Geometry::Type geom_type) const2759 ::GetTraceElement(int i, Geometry::Type geom_type) const
2760 {
2761 return fec->TraceFiniteElementForGeometry(geom_type);
2762 }
2763
~FiniteElementSpace()2764 FiniteElementSpace::~FiniteElementSpace()
2765 {
2766 Destroy();
2767 }
2768
Destroy()2769 void FiniteElementSpace::Destroy()
2770 {
2771 delete cR;
2772 delete cR_hp;
2773 delete cP;
2774 Th.Clear();
2775 L2E_nat.Clear();
2776 L2E_lex.Clear();
2777 for (int i = 0; i < E2Q_array.Size(); i++)
2778 {
2779 delete E2Q_array[i];
2780 }
2781 E2Q_array.SetSize(0);
2782 for (auto &x : L2F)
2783 {
2784 delete x.second;
2785 }
2786 for (int i = 0; i < E2IFQ_array.Size(); i++)
2787 {
2788 delete E2IFQ_array[i];
2789 }
2790 E2IFQ_array.SetSize(0);
2791 for (int i = 0; i < E2BFQ_array.Size(); i++)
2792 {
2793 delete E2BFQ_array[i];
2794 }
2795 E2BFQ_array.SetSize(0);
2796
2797 dof_elem_array.DeleteAll();
2798 dof_ldof_array.DeleteAll();
2799
2800 if (NURBSext)
2801 {
2802 if (own_ext) { delete NURBSext; }
2803 delete face_dof;
2804 face_to_be.DeleteAll();
2805 }
2806 else
2807 {
2808 delete elem_dof;
2809 delete bdr_elem_dof;
2810 delete face_dof;
2811
2812 delete [] bdofs;
2813 }
2814 ceed::RemoveBasisAndRestriction(this);
2815 }
2816
GetTransferOperator(const FiniteElementSpace & coarse_fes,OperatorHandle & T) const2817 void FiniteElementSpace::GetTransferOperator(
2818 const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
2819 {
2820 // Assumptions: see the declaration of the method.
2821
2822 if (T.Type() == Operator::MFEM_SPARSEMAT)
2823 {
2824 Mesh::GeometryList elem_geoms(*mesh);
2825
2826 DenseTensor localP[Geometry::NumGeom];
2827 for (int i = 0; i < elem_geoms.Size(); i++)
2828 {
2829 GetLocalRefinementMatrices(coarse_fes, elem_geoms[i],
2830 localP[elem_geoms[i]]);
2831 }
2832 T.Reset(RefinementMatrix_main(coarse_fes.GetNDofs(),
2833 coarse_fes.GetElementToDofTable(),
2834 localP));
2835 }
2836 else
2837 {
2838 T.Reset(new RefinementOperator(this, &coarse_fes));
2839 }
2840 }
2841
GetTrueTransferOperator(const FiniteElementSpace & coarse_fes,OperatorHandle & T) const2842 void FiniteElementSpace::GetTrueTransferOperator(
2843 const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
2844 {
2845 const SparseMatrix *coarse_P = coarse_fes.GetConformingProlongation();
2846
2847 Operator::Type req_type = T.Type();
2848 GetTransferOperator(coarse_fes, T);
2849
2850 if (req_type == Operator::MFEM_SPARSEMAT)
2851 {
2852 if (GetConformingRestriction())
2853 {
2854 T.Reset(mfem::Mult(*cR, *T.As<SparseMatrix>()));
2855 }
2856 if (coarse_P)
2857 {
2858 T.Reset(mfem::Mult(*T.As<SparseMatrix>(), *coarse_P));
2859 }
2860 }
2861 else
2862 {
2863 const int RP_case = bool(GetConformingRestriction()) + 2*bool(coarse_P);
2864 if (RP_case == 0) { return; }
2865 const bool owner = T.OwnsOperator();
2866 T.SetOperatorOwner(false);
2867 switch (RP_case)
2868 {
2869 case 1:
2870 T.Reset(new ProductOperator(cR, T.Ptr(), false, owner));
2871 break;
2872 case 2:
2873 T.Reset(new ProductOperator(T.Ptr(), coarse_P, owner, false));
2874 break;
2875 case 3:
2876 T.Reset(new TripleProductOperator(
2877 cR, T.Ptr(), coarse_P, false, owner, false));
2878 break;
2879 }
2880 }
2881 }
2882
UpdateElementOrders()2883 void FiniteElementSpace::UpdateElementOrders()
2884 {
2885 const CoarseFineTransformations &cf_tr = mesh->GetRefinementTransforms();
2886
2887 Array<char> new_order(mesh->GetNE());
2888 switch (mesh->GetLastOperation())
2889 {
2890 case Mesh::REFINE:
2891 {
2892 for (int i = 0; i < mesh->GetNE(); i++)
2893 {
2894 new_order[i] = elem_order[cf_tr.embeddings[i].parent];
2895 }
2896 break;
2897 }
2898 default:
2899 MFEM_ABORT("not implemented yet");
2900 }
2901
2902 mfem::Swap(elem_order, new_order);
2903 }
2904
Update(bool want_transform)2905 void FiniteElementSpace::Update(bool want_transform)
2906 {
2907 if (!orders_changed)
2908 {
2909 if (mesh->GetSequence() == mesh_sequence)
2910 {
2911 return; // mesh and space are in sync, no-op
2912 }
2913 if (want_transform && mesh->GetSequence() != mesh_sequence + 1)
2914 {
2915 MFEM_ABORT("Error in update sequence. Space needs to be updated after "
2916 "each mesh modification.");
2917 }
2918 }
2919 else
2920 {
2921 if (mesh->GetSequence() != mesh_sequence)
2922 {
2923 MFEM_ABORT("Updating space after both mesh change and element order "
2924 "change is not supported. Please update separately after "
2925 "each change.");
2926 }
2927 }
2928
2929 if (NURBSext)
2930 {
2931 UpdateNURBS();
2932 return;
2933 }
2934
2935 Table* old_elem_dof = NULL;
2936 int old_ndofs;
2937 bool old_orders_changed = orders_changed;
2938
2939 // save old DOF table
2940 if (want_transform)
2941 {
2942 old_elem_dof = elem_dof;
2943 elem_dof = NULL;
2944 old_ndofs = ndofs;
2945 }
2946
2947 // update the 'elem_order' array if the mesh has changed
2948 if (IsVariableOrder() && mesh->GetSequence() != mesh_sequence)
2949 {
2950 UpdateElementOrders();
2951 }
2952
2953 Destroy(); // calls Th.Clear()
2954 Construct();
2955 BuildElementToDofTable();
2956
2957 if (want_transform)
2958 {
2959 MFEM_VERIFY(!old_orders_changed, "Interpolation for element order change "
2960 "is not implemented yet, sorry.");
2961
2962 // calculate appropriate GridFunction transformation
2963 switch (mesh->GetLastOperation())
2964 {
2965 case Mesh::REFINE:
2966 {
2967 if (Th.Type() != Operator::MFEM_SPARSEMAT)
2968 {
2969 Th.Reset(new RefinementOperator(this, old_elem_dof, old_ndofs));
2970 // The RefinementOperator takes ownership of 'old_elem_dof', so
2971 // we no longer own it:
2972 old_elem_dof = NULL;
2973 }
2974 else
2975 {
2976 // calculate fully assembled matrix
2977 Th.Reset(RefinementMatrix(old_ndofs, old_elem_dof));
2978 }
2979 break;
2980 }
2981
2982 case Mesh::DEREFINE:
2983 {
2984 BuildConformingInterpolation();
2985 Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof));
2986 if (cP && cR)
2987 {
2988 Th.SetOperatorOwner(false);
2989 Th.Reset(new TripleProductOperator(cP, cR, Th.Ptr(),
2990 false, false, true));
2991 }
2992 break;
2993 }
2994
2995 default:
2996 break;
2997 }
2998
2999 delete old_elem_dof;
3000 }
3001 }
3002
UpdateMeshPointer(Mesh * new_mesh)3003 void FiniteElementSpace::UpdateMeshPointer(Mesh *new_mesh)
3004 {
3005 mesh = new_mesh;
3006 }
3007
Save(std::ostream & out) const3008 void FiniteElementSpace::Save(std::ostream &out) const
3009 {
3010 int fes_format = 90; // the original format, v0.9
3011 bool nurbs_unit_weights = false;
3012
3013 // Determine the format that should be used.
3014 if (!NURBSext)
3015 {
3016 // TODO: if this is a variable-order FE space, use fes_format = 100.
3017 }
3018 else
3019 {
3020 const NURBSFECollection *nurbs_fec =
3021 dynamic_cast<const NURBSFECollection *>(fec);
3022 MFEM_VERIFY(nurbs_fec, "invalid FE collection");
3023 nurbs_fec->SetOrder(NURBSext->GetOrder());
3024 const double eps = 5e-14;
3025 nurbs_unit_weights = (NURBSext->GetWeights().Min() >= 1.0-eps &&
3026 NURBSext->GetWeights().Max() <= 1.0+eps);
3027 if ((NURBSext->GetOrder() == NURBSFECollection::VariableOrder) ||
3028 (NURBSext != mesh->NURBSext && !nurbs_unit_weights) ||
3029 (NURBSext->GetMaster().Size() != 0 ))
3030 {
3031 fes_format = 100; // v1.0 format
3032 }
3033 }
3034
3035 out << (fes_format == 90 ?
3036 "FiniteElementSpace\n" : "MFEM FiniteElementSpace v1.0\n")
3037 << "FiniteElementCollection: " << fec->Name() << '\n'
3038 << "VDim: " << vdim << '\n'
3039 << "Ordering: " << ordering << '\n';
3040
3041 if (fes_format == 100) // v1.0
3042 {
3043 if (!NURBSext)
3044 {
3045 // TODO: this is a variable-order FE space --> write 'element_orders'.
3046 }
3047 else if (NURBSext != mesh->NURBSext)
3048 {
3049 if (NURBSext->GetOrder() != NURBSFECollection::VariableOrder)
3050 {
3051 out << "NURBS_order\n" << NURBSext->GetOrder() << '\n';
3052 }
3053 else
3054 {
3055 out << "NURBS_orders\n";
3056 // 1 = do not write the size, just the entries:
3057 NURBSext->GetOrders().Save(out, 1);
3058 }
3059 // If periodic BCs are given, write connectivity
3060 if (NURBSext->GetMaster().Size() != 0 )
3061 {
3062 out <<"NURBS_periodic\n";
3063 NURBSext->GetMaster().Save(out);
3064 NURBSext->GetSlave().Save(out);
3065 }
3066 // If the weights are not unit, write them to the output:
3067 if (!nurbs_unit_weights)
3068 {
3069 out << "NURBS_weights\n";
3070 NURBSext->GetWeights().Print(out, 1);
3071 }
3072 }
3073 out << "End: MFEM FiniteElementSpace v1.0\n";
3074 }
3075 }
3076
Load(Mesh * m,std::istream & input)3077 FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)
3078 {
3079 string buff;
3080 int fes_format = 0, ord;
3081 FiniteElementCollection *r_fec;
3082
3083 Destroy();
3084
3085 input >> std::ws;
3086 getline(input, buff); // 'FiniteElementSpace'
3087 filter_dos(buff);
3088 if (buff == "FiniteElementSpace") { fes_format = 90; /* v0.9 */ }
3089 else if (buff == "MFEM FiniteElementSpace v1.0") { fes_format = 100; }
3090 else { MFEM_ABORT("input stream is not a FiniteElementSpace!"); }
3091 getline(input, buff, ' '); // 'FiniteElementCollection:'
3092 input >> std::ws;
3093 getline(input, buff);
3094 filter_dos(buff);
3095 r_fec = FiniteElementCollection::New(buff.c_str());
3096 getline(input, buff, ' '); // 'VDim:'
3097 input >> vdim;
3098 getline(input, buff, ' '); // 'Ordering:'
3099 input >> ord;
3100
3101 NURBSFECollection *nurbs_fec = dynamic_cast<NURBSFECollection*>(r_fec);
3102 NURBSExtension *NURBSext = NULL;
3103 if (fes_format == 90) // original format, v0.9
3104 {
3105 if (nurbs_fec)
3106 {
3107 MFEM_VERIFY(m->NURBSext, "NURBS FE collection requires a NURBS mesh!");
3108 const int order = nurbs_fec->GetOrder();
3109 if (order != m->NURBSext->GetOrder() &&
3110 order != NURBSFECollection::VariableOrder)
3111 {
3112 NURBSext = new NURBSExtension(m->NURBSext, order);
3113 }
3114 }
3115 }
3116 else if (fes_format == 100) // v1.0
3117 {
3118 while (1)
3119 {
3120 skip_comment_lines(input, '#');
3121 MFEM_VERIFY(input.good(), "error reading FiniteElementSpace v1.0");
3122 getline(input, buff);
3123 filter_dos(buff);
3124 if (buff == "NURBS_order" || buff == "NURBS_orders")
3125 {
3126 MFEM_VERIFY(nurbs_fec,
3127 buff << ": NURBS FE collection is required!");
3128 MFEM_VERIFY(m->NURBSext, buff << ": NURBS mesh is required!");
3129 MFEM_VERIFY(!NURBSext, buff << ": order redefinition!");
3130 if (buff == "NURBS_order")
3131 {
3132 int order;
3133 input >> order;
3134 NURBSext = new NURBSExtension(m->NURBSext, order);
3135 }
3136 else
3137 {
3138 Array<int> orders;
3139 orders.Load(m->NURBSext->GetNKV(), input);
3140 NURBSext = new NURBSExtension(m->NURBSext, orders);
3141 }
3142 }
3143 else if (buff == "NURBS_periodic")
3144 {
3145 Array<int> master, slave;
3146 master.Load(input);
3147 slave.Load(input);
3148 NURBSext->ConnectBoundaries(master,slave);
3149 }
3150 else if (buff == "NURBS_weights")
3151 {
3152 MFEM_VERIFY(NURBSext, "NURBS_weights: NURBS_orders have to be "
3153 "specified before NURBS_weights!");
3154 NURBSext->GetWeights().Load(input, NURBSext->GetNDof());
3155 }
3156 else if (buff == "element_orders")
3157 {
3158 MFEM_VERIFY(!nurbs_fec, "section element_orders cannot be used "
3159 "with a NURBS FE collection");
3160 MFEM_ABORT("element_orders: not implemented yet!");
3161 }
3162 else if (buff == "End: MFEM FiniteElementSpace v1.0")
3163 {
3164 break;
3165 }
3166 else
3167 {
3168 MFEM_ABORT("unknown section: " << buff);
3169 }
3170 }
3171 }
3172
3173 Constructor(m, NURBSext, r_fec, vdim, ord);
3174
3175 return r_fec;
3176 }
3177
3178
Construct()3179 void QuadratureSpace::Construct()
3180 {
3181 // protected method
3182 int offset = 0;
3183 const int num_elem = mesh->GetNE();
3184 element_offsets = new int[num_elem + 1];
3185 for (int g = 0; g < Geometry::NumGeom; g++)
3186 {
3187 int_rule[g] = NULL;
3188 }
3189 for (int i = 0; i < num_elem; i++)
3190 {
3191 element_offsets[i] = offset;
3192 int geom = mesh->GetElementBaseGeometry(i);
3193 if (int_rule[geom] == NULL)
3194 {
3195 int_rule[geom] = &IntRules.Get(geom, order);
3196 }
3197 offset += int_rule[geom]->GetNPoints();
3198 }
3199 element_offsets[num_elem] = size = offset;
3200 }
3201
QuadratureSpace(Mesh * mesh_,std::istream & in)3202 QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
3203 : mesh(mesh_)
3204 {
3205 const char *msg = "invalid input stream";
3206 string ident;
3207
3208 in >> ident; MFEM_VERIFY(ident == "QuadratureSpace", msg);
3209 in >> ident; MFEM_VERIFY(ident == "Type:", msg);
3210 in >> ident;
3211 if (ident == "default_quadrature")
3212 {
3213 in >> ident; MFEM_VERIFY(ident == "Order:", msg);
3214 in >> order;
3215 }
3216 else
3217 {
3218 MFEM_ABORT("unknown QuadratureSpace type: " << ident);
3219 return;
3220 }
3221
3222 Construct();
3223 }
3224
Save(std::ostream & out) const3225 void QuadratureSpace::Save(std::ostream &out) const
3226 {
3227 out << "QuadratureSpace\n"
3228 << "Type: default_quadrature\n"
3229 << "Order: " << order << '\n';
3230 }
3231
3232 } // namespace mfem
3233