14 #include "../general/text.hpp"
15 #include "../general/forall.hpp"
16 #include "../mesh/mesh_headers.hpp"
18 #include "ceed/util.hpp"
28 template <>
void Ordering::
29 DofsToVDofs<Ordering::byNODES>(
int ndofs,
int vdim,
Array<int> &dofs)
32 int size = dofs.Size();
33 dofs.SetSize(size*vdim);
34 for (
int vd = 1; vd < vdim; vd++)
36 for (
int i = 0; i < size; i++)
38 dofs[i+size*vd] = Map<byNODES>(ndofs, vdim, dofs[i], vd);
43 template <>
void Ordering::
44 DofsToVDofs<Ordering::byVDIM>(
int ndofs,
int vdim,
Array<int> &dofs)
47 int size = dofs.Size();
48 dofs.SetSize(size*vdim);
49 for (
int vd = vdim-1; vd >= 0; vd--)
51 for (
int i = 0; i < size; i++)
53 dofs[i+size*vd] = Map<byVDIM>(ndofs, vdim, dofs[i], vd);
59 FiniteElementSpace::FiniteElementSpace()
61 ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0),
63 elem_dof(NULL), elem_fos(NULL), bdr_elem_dof(NULL), bdr_elem_fos(NULL),
65 NURBSext(NULL), own_ext(false),
66 DoFTrans(0), VDoFTrans(vdim, ordering),
67 cP(NULL), cR(NULL), cR_hp(NULL), cP_is_set(false),
69 sequence(0), mesh_sequence(0), orders_changed(false), relaxed_hp(false)
75 : VDoFTrans(orig.vdim, orig.ordering)
77 mesh_ = mesh_ ? mesh_ : orig.
mesh;
78 fec_ = fec_ ? fec_ : orig.
fec;
103 MFEM_VERIFY(
cP == NULL,
"");
104 MFEM_VERIFY(
cR == NULL,
"");
112 int n = perm->
Size();
114 for (
int i=0; i<n; ++i)
118 perm_mat->
Set(i, j, s);
130 else if (perm != NULL)
141 else if (perm != NULL)
154 "Space has not been Updated() after a Mesh change.");
155 MFEM_VERIFY(i >= 0 && i <
GetNE(),
"Invalid element index");
156 MFEM_VERIFY(p >= 0 && p <=
MaxVarOrder,
"Order out of range");
181 "Space has not been Updated() after a Mesh change.");
182 MFEM_VERIFY(i >= 0 && i <
GetNE(),
"Invalid element index");
197 if (ndofs_ < 0) { ndofs_ = this->
ndofs; }
201 for (
int i = 0; i < dofs.
Size(); i++)
203 dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs_,
vdim, i, vd);
208 for (
int i = 0; i < dofs.
Size(); i++)
210 dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs_,
vdim, i, vd);
217 if (
vdim == 1) {
return; }
218 if (ndofs_ < 0) { ndofs_ = this->
ndofs; }
222 Ordering::DofsToVDofs<Ordering::byNODES>(ndofs_,
vdim, dofs);
226 Ordering::DofsToVDofs<Ordering::byVDIM>(ndofs_,
vdim, dofs);
232 if (
vdim == 1) {
return; }
233 if (ndofs_ < 0) { ndofs_ = this->
ndofs; }
237 for (
int i = 0; i < dofs.
Size(); i++)
239 dofs[i] = Ordering::Map<Ordering::byNODES>(ndofs_,
vdim, dofs[i], vd);
244 for (
int i = 0; i < dofs.
Size(); i++)
246 dofs[i] = Ordering::Map<Ordering::byVDIM>(ndofs_,
vdim, dofs[i], vd);
253 if (
vdim == 1) {
return dof; }
254 if (ndofs_ < 0) { ndofs_ = this->
ndofs; }
258 return Ordering::Map<Ordering::byNODES>(ndofs_,
vdim, dof, vd);
262 return Ordering::Map<Ordering::byVDIM>(ndofs_,
vdim, dof, vd);
269 int n = vdofs.
Size(), *vdof = vdofs;
270 for (
int i = 0; i < n; i++)
273 if ((j = vdof[i]) < 0)
285 if (
vdim == 1 || doftrans == NULL)
301 if (
vdim == 1 || doftrans == NULL)
352 if (el_fos) { el_fos -> MakeI (
mesh ->
GetNE()); }
353 for (
int i = 0; i <
mesh ->
GetNE(); i++)
356 el_dof -> AddColumnsInRow (i, dofs.
Size());
361 el_fos -> AddColumnsInRow (i, Fo.
Size());
365 if (el_fos) { el_fos -> MakeJ(); }
366 for (
int i = 0; i <
mesh ->
GetNE(); i++)
369 el_dof -> AddConnections (i, (
int *)dofs, dofs.
Size());
374 el_fos -> AddConnections (i, (
int *)Fo, Fo.
Size());
377 el_dof -> ShiftUpI();
378 if (el_fos) { el_fos -> ShiftUpI(); }
416 for (
int i = 0; i < fc_dof->
Size(); i++)
422 for (
int i = 0; i < fc_dof->
Size(); i++)
447 for (
int k = 0, dof_counter = 0; k < nnz; k++)
449 const int sdof = J[k];
450 const int dof = (sdof < 0) ? -1-sdof : sdof;
451 int new_dof = dof_marker[dof];
454 dof_marker[dof] = new_dof = dof_counter++;
456 J[k] = (sdof < 0) ? -1-new_dof : new_dof;
469 for (
int i = 0; i <
mesh ->
GetNE(); i++)
471 const int *dofs =
elem_dof -> GetRow(i);
472 const int n =
elem_dof -> RowSize(i);
473 for (
int j = 0; j < n; j++)
487 for (
int i = 0; i < dofs.
Size(); i++)
490 if (k < 0) { k = -1 - k; }
504 for (
int i = 0; i <
GetNBE(); i++)
512 mark_dofs(vdofs, ess_vdofs);
517 for (
int d = 0; d < dofs.
Size(); d++)
518 { dofs[d] =
DofToVDof(dofs[d], component); }
519 mark_dofs(dofs, ess_vdofs);
531 for (
int i = 0; i < bdr_verts.
Size(); i++)
536 mark_dofs(vdofs, ess_vdofs);
541 for (
int d = 0; d < dofs.
Size(); d++)
542 { dofs[d] =
DofToVDof(dofs[d], component); }
543 mark_dofs(dofs, ess_vdofs);
546 for (
int i = 0; i < bdr_edges.
Size(); i++)
551 mark_dofs(vdofs, ess_vdofs);
556 for (
int d = 0; d < dofs.
Size(); d++)
557 { dofs[d] =
DofToVDof(dofs[d], component); }
558 mark_dofs(dofs, ess_vdofs);
603 for (
int i = 0; i < marker.
Size(); i++)
605 if (marker[i]) { num_marked++; }
610 for (
int i = 0; i < marker.
Size(); i++)
612 if (marker[i]) { list.
Append(i); }
624 for (
int i = 0; i < list.
Size(); i++)
626 marker[list[i]] = mark_val;
635 else { dofs.
Copy(cdofs); }
643 else { cdofs.
Copy(dofs); }
661 if (d_vdofs.
Size() != c_vdofs.
Size())
663 mfem_error (
"FiniteElementSpace::D2C_GlobalRestrictionMatrix (...)");
667 for (j = 0; j < d_vdofs.
Size(); j++)
669 R ->
Set (c_vdofs[j], d_vdofs[j], 1.0);
693 if (c_dofs.
Size() != 1)
695 "D2Const_GlobalRestrictionMatrix (...)");
698 for (j = 0; j < d_dofs.
Size(); j++)
700 R ->
Set (c_dofs[0], d_dofs[j], 1.0);
724 for (
int i = 0; i <
mesh ->
GetNE(); i++)
731 if (geom != cached_geom)
733 h_fe =
this ->
GetFE (i);
734 l_fe = lfes ->
GetFE (i);
736 h_fe->
Project(*l_fe, T, loc_restr);
740 for (
int vd = 0; vd < lvdim; vd++)
742 l_dofs.
Copy(l_vdofs);
745 h_dofs.
Copy(h_vdofs);
748 R -> SetSubMatrix (l_vdofs, h_vdofs, loc_restr, 1);
761 for (
int i = skipfirst; i < slave_dofs.
Size(); i++)
763 int sdof = slave_dofs[i];
766 for (
int j = 0; j < master_dofs.
Size(); j++)
768 double coef = I(i, j);
769 if (std::abs(coef) > 1e-12)
771 int mdof = master_dofs[j];
772 if (mdof != sdof && mdof != (-1-sdof))
774 deps.
Add(sdof, mdof, coef);
797 mesh->GetFaceVertices(slave_face, V);
798 mesh->GetFaceEdges(slave_face, E, Eo);
799 MFEM_ASSERT(V.
Size() == E.
Size(),
"");
806 for (
int i = 0; i < E.
Size(); i++)
808 int a = i,
b = (i+1) % V.
Size();
809 if (V[a] > V[b]) { std::swap(a, b); }
816 for (
int j = 0; j < 2; j++)
818 edge_pm(j, 0) = (*pm)(j,
a);
819 edge_pm(j, 1) = (*pm)(j,
b);
820 mid[j] = 0.5*((*pm)(j,
a) + (*pm)(j,
b));
824 const double eps = 1e-14;
825 if (mid[0] > eps && mid[0] < 1-eps &&
826 mid[1] > eps && mid[1] < 1-eps)
828 int order = GetEdgeDofs(E[i], slave_dofs, 0);
831 edge_fe->GetTransferMatrix(*master_fe, edge_T, I);
833 AddDependencies(deps, master_dofs, slave_dofs, I, 0);
845 for (
int i = 0; i < ndep; i++)
847 if (!finalized[dep[i]]) {
return false; }
867 int order =
GetEdgeDofs(-1 - index, edof, variant);
874 if (!dofs.
Size()) {
return 0; }
879 for (
int i = 0; i < nv; i++)
882 dofs[nv+i] = edof[nv+i];
886 for (
int i = 0; i < ne; i++)
888 dofs[face_vert*nv + i] = edof[2*nv + i];
931 MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(
this) == NULL,
932 "This method should not be used with a ParFiniteElementSpace!");
944 Array<int> master_dofs, slave_dofs, highest_dofs;
961 for (
int entity = 2; entity >= 1; entity--)
964 if (!list.
masters.Size()) {
continue; }
972 if (!master_dofs.
Size()) {
continue; }
975 if (!master_fe) {
continue; }
982 default: MFEM_ABORT(
"unsupported geometry");
990 if (!slave_dofs.
Size()) {
break; }
1011 slave_dofs, slave.
index, pm);
1023 master_geom, nvar-1);
1024 const auto *highest_fe =
fec->
GetFE(master_geom, q);
1044 MFEM_ASSERT(ent_dofs.
Size() == num_ent+1,
"");
1048 for (
int i = 0; i < num_ent; i++)
1050 if (ent_dofs.
RowSize(i) <= 1) {
continue; }
1055 if (geom != last_geom)
1063 const auto *master_fe =
fec->
GetFE(geom, p);
1064 if (!master_fe) {
break; }
1067 for (
int variant = 1; ; variant++)
1069 int q =
GetEntityDofs(entity, i, slave_dofs, geom, variant);
1070 if (q < 0) {
break; }
1072 const auto *slave_fe =
fec->
GetFE(geom, q);
1085 int n_true_dofs = 0;
1086 for (
int i = 0; i <
ndofs; i++)
1088 if (!deps.
RowSize(i)) { n_true_dofs++; }
1092 if (n_true_dofs == ndofs)
1107 for (
int i = 0; i < n_true_dofs; i++)
1112 cR_I[n_true_dofs] = n_true_dofs;
1127 for (
int i = 0, true_dof = 0; i <
ndofs; i++)
1131 cP->
Add(i, true_dof, 1.0);
1133 finalized[i] =
true;
1139 inv_deps.
GetRow(i, cols, srow);
1160 int n_finalized = n_true_dofs;
1164 for (
int dof = 0; dof <
ndofs; dof++)
1170 int n_dep = deps.
RowSize(dof);
1172 for (
int j = 0; j < n_dep; j++)
1174 cP->
GetRow(dep_col[j], cols, srow);
1175 srow *= dep_coef[j];
1179 finalized[dof] =
true;
1189 MFEM_VERIFY(n_finalized == ndofs,
1190 "Error creating cP matrix: n_finalized = "
1191 << n_finalized <<
", ndofs = " << ndofs);
1208 if (
vdim == 1) {
return; }
1210 int height = mat.
Height();
1211 int width = mat.
Width();
1217 for (
int i = 0; i < height; i++)
1219 mat.
GetRow(i, dofs, srow);
1220 for (
int vd = 0; vd <
vdim; vd++)
1300 key_face key = std::make_tuple(is_dg_space, e_ordering, type, m);
1301 auto itr =
L2F.find(key);
1302 if (itr !=
L2F.end())
1332 for (
int i = 0; i <
E2Q_array.Size(); i++)
1335 if (qi->
IntRule == &ir) {
return qi; }
1346 for (
int i = 0; i <
E2Q_array.Size(); i++)
1349 if (qi->
qspace == &qs) {
return qi; }
1366 if (qi->
IntRule == &ir) {
return qi; }
1379 if (qi->
IntRule == &ir) {
return qi; }
1390 const int coarse_ndofs,
const Table &coarse_elem_dof,
1403 if (elem_geoms.
Size() == 1)
1405 const int coarse_ldof = localP[elem_geoms[0]].
SizeJ();
1423 const int fine_ldof = localP[geom].
SizeI();
1428 for (
int vd = 0; vd <
vdim; vd++)
1430 coarse_dofs.Copy(coarse_vdofs);
1433 for (
int i = 0; i < fine_ldof; i++)
1436 int m = (r >= 0) ? r : (-1 - r);
1441 P->
SetRow(r, coarse_vdofs, row);
1448 MFEM_ASSERT(mark.
Sum() == P->
Height(),
"Not all rows of P set.");
1461 int nmat = pmats.
SizeK();
1468 localP.
SetSize(ldof, ldof, nmat);
1469 for (
int i = 0; i < nmat; i++)
1477 const Table* old_elem_dof,
1478 const Table* old_elem_fos)
1480 MFEM_VERIFY(
GetNE() >= old_elem_dof->
Size(),
1481 "Previous mesh is not coarser.");
1486 for (
int i = 0; i < elem_geoms.Size(); i++)
1499 , old_elem_dof(old_elem_dof)
1500 , old_elem_fos(old_elem_fos)
1502 MFEM_VERIFY(fespace->
GetNE() >= old_elem_dof->
Size(),
1503 "Previous mesh is not coarser.");
1505 width = old_ndofs * fespace->
GetVDim();
1510 for (
int i = 0; i < elem_geoms.Size(); i++)
1515 ConstructDoFTrans();
1520 :
Operator(fespace->GetVSize(), coarse_fes->GetVSize()),
1521 fespace(fespace), old_elem_dof(NULL), old_elem_fos(NULL)
1525 for (
int i = 0; i < elem_geoms.Size(); i++)
1528 localP[elem_geoms[i]]);
1540 ConstructDoFTrans();
1545 delete old_elem_dof;
1546 delete old_elem_fos;
1549 void FiniteElementSpace::RefinementOperator
1550 ::ConstructDoFTrans()
1553 for (
int i=0; i<old_DoFTrans.Size(); i++)
1555 old_DoFTrans[i] = NULL;
1559 if (dynamic_cast<const ND_FECollection*>(fec_ref))
1569 const FiniteElement * nd_tet =
1574 new ND_TetDofTransformation(nd_tet->GetOrder());
1577 const FiniteElement * nd_pri =
1582 new ND_WedgeDofTransformation(nd_pri->GetOrder());
1594 Array<int> dofs, vdofs, old_dofs, old_vdofs, old_Fo;
1596 int rvdim = fespace->
GetVDim();
1597 int old_ndofs = width / rvdim;
1601 for (
int k = 0; k < mesh_ref->
GetNE(); k++)
1614 for (
int vd = 0; vd < rvdim; vd++)
1618 old_dofs.
Copy(old_vdofs);
1621 lP.
Mult(subX, subY);
1628 old_DoFTrans[geom]->SetFaceOrientations(old_Fo);
1635 new_doftrans = doftrans;
1639 for (
int vd = 0; vd < rvdim; vd++)
1643 old_dofs.
Copy(old_vdofs);
1646 old_DoFTrans[geom]->InvTransformPrimal(subX);
1647 lP.
Mult(subX, subY);
1654 doftrans = new_doftrans;
1672 Array<int> f_dofs, c_dofs, f_vdofs, c_vdofs, old_Fo;
1674 int rvdim = fespace->
GetVDim();
1675 int old_ndofs = width / rvdim;
1677 Vector subY, subX, subYt;
1679 for (
int k = 0; k < mesh_ref->
GetNE(); k++)
1692 for (
int vd = 0; vd < rvdim; vd++)
1694 f_dofs.
Copy(f_vdofs);
1696 c_dofs.
Copy(c_vdofs);
1701 for (
int p = 0;
p < f_dofs.
Size(); ++
p)
1718 old_DoFTrans[geom]->SetFaceOrientations(old_Fo);
1725 new_doftrans = doftrans;
1729 for (
int vd = 0; vd < rvdim; vd++)
1731 f_dofs.
Copy(f_vdofs);
1733 c_dofs.
Copy(c_vdofs);
1737 doftrans->InvTransformDual(subX);
1738 for (
int p = 0;
p < f_dofs.
Size(); ++
p)
1747 old_DoFTrans[geom]->TransformDual(subYt);
1753 doftrans = new_doftrans;
1757 for (
int p = 0;
p < f_dofs.
Size(); ++
p)
1775 : geom(g), num_children(n), children(c) { }
1777 bool operator<(
const RefType &other)
const
1779 if (geom < other.geom) {
return true; }
1780 if (geom > other.geom) {
return false; }
1781 if (num_children < other.num_children) {
return true; }
1782 if (num_children > other.num_children) {
return false; }
1783 for (
int i = 0; i < num_children; i++)
1785 if (children[i].one < other.children[i].one) {
return true; }
1786 if (children[i].one > other.children[i].one) {
return false; }
1792 void GetCoarseToFineMap(
const CoarseFineTransformations &cft,
1794 Table &coarse_to_fine,
1795 Array<int> &coarse_to_ref_type,
1796 Table &ref_type_to_matrix,
1797 Array<Geometry::Type> &ref_type_to_geom)
1799 const int fine_ne = cft.embeddings.Size();
1801 for (
int i = 0; i < fine_ne; i++)
1803 coarse_ne = std::max(coarse_ne, cft.embeddings[i].parent);
1807 coarse_to_ref_type.SetSize(coarse_ne);
1808 coarse_to_fine.SetDims(coarse_ne, fine_ne);
1810 Array<int> cf_i(coarse_to_fine.GetI(), coarse_ne+1);
1811 Array<Pair<int,int> > cf_j(fine_ne);
1813 for (
int i = 0; i < fine_ne; i++)
1815 cf_i[cft.embeddings[i].parent+1]++;
1818 MFEM_ASSERT(cf_i.Last() == cf_j.Size(),
"internal error");
1819 for (
int i = 0; i < fine_ne; i++)
1821 const Embedding &e = cft.embeddings[i];
1822 cf_j[cf_i[e.parent]].one = e.matrix;
1823 cf_j[cf_i[e.parent]].two = i;
1826 std::copy_backward(cf_i.begin(), cf_i.end()-1, cf_i.end());
1828 for (
int i = 0; i < coarse_ne; i++)
1830 std::sort(&cf_j[cf_i[i]], cf_j.GetData() + cf_i[i+1]);
1832 for (
int i = 0; i < fine_ne; i++)
1834 coarse_to_fine.GetJ()[i] = cf_j[i].two;
1840 map<RefType,int> ref_type_map;
1841 for (
int i = 0; i < coarse_ne; i++)
1843 const int num_children = cf_i[i+1]-cf_i[i];
1844 MFEM_ASSERT(num_children > 0,
"");
1845 const int fine_el = cf_j[cf_i[i]].two;
1848 const RefType ref_type(geom, num_children, &cf_j[cf_i[i]]);
1849 pair<map<RefType,int>::iterator,
bool> res =
1850 ref_type_map.insert(
1851 pair<const RefType,int>(ref_type, (
int)ref_type_map.size()));
1852 coarse_to_ref_type[i] = res.first->second;
1855 ref_type_to_matrix.MakeI((
int)ref_type_map.size());
1856 ref_type_to_geom.SetSize((
int)ref_type_map.size());
1857 for (map<RefType,int>::iterator it = ref_type_map.begin();
1858 it != ref_type_map.end(); ++it)
1860 ref_type_to_matrix.AddColumnsInRow(it->second, it->first.num_children);
1861 ref_type_to_geom[it->second] = it->first.geom;
1864 ref_type_to_matrix.MakeJ();
1865 for (map<RefType,int>::iterator it = ref_type_map.begin();
1866 it != ref_type_map.end(); ++it)
1868 const RefType &rt = it->first;
1869 for (
int j = 0; j < rt.num_children; j++)
1871 ref_type_to_matrix.AddConnection(it->second, rt.children[j].one);
1874 ref_type_to_matrix.ShiftUpI();
1889 "incompatible coarse and fine FE spaces");
1897 for (
int gi = 0; gi < elem_geoms.
Size(); gi++)
1900 DenseTensor &lP = localP[geom], &lM = localM[geom];
1909 emb_tr.SetIdentityTransformation(geom);
1910 for (
int i = 0; i < pmats.
SizeK(); i++)
1912 emb_tr.SetPointMat(pmats(i));
1920 Table ref_type_to_matrix;
1921 internal::GetCoarseToFineMap(rtrans, *f_mesh, coarse_to_fine,
1922 coarse_to_ref_type, ref_type_to_matrix,
1924 MFEM_ASSERT(coarse_to_fine.Size() == c_fes->
GetNE(),
"");
1926 const int total_ref_types = ref_type_to_geom.Size();
1928 Array<int> ref_type_to_coarse_elem_offset(total_ref_types);
1929 ref_type_to_fine_elem_offset.
SetSize(total_ref_types);
1932 for (
int i = 0; i < total_ref_types; i++)
1935 ref_type_to_coarse_elem_offset[i] = num_ref_types[g];
1936 ref_type_to_fine_elem_offset[i] = num_fine_elems[g];
1938 num_fine_elems[g] += ref_type_to_matrix.
RowSize(i);
1943 if (num_ref_types[g] == 0) {
continue; }
1944 const int fine_dofs = localP[g].
SizeI();
1945 const int coarse_dofs = localP[g].
SizeJ();
1946 localPtMP[g].
SetSize(coarse_dofs, coarse_dofs, num_ref_types[g]);
1947 localR[g].
SetSize(coarse_dofs, fine_dofs, num_fine_elems[g]);
1949 for (
int i = 0; i < total_ref_types; i++)
1952 DenseMatrix &lPtMP = localPtMP[g](ref_type_to_coarse_elem_offset[i]);
1953 int lR_offset = ref_type_to_fine_elem_offset[i];
1954 const int *mi = ref_type_to_matrix.
GetRow(i);
1955 const int nm = ref_type_to_matrix.
RowSize(i);
1957 for (
int s = 0;
s < nm;
s++)
1966 for (
int s = 0;
s < nm;
s++)
1979 delete coarse_elem_dof;
1988 const int fine_vdim = fine_fes->GetVDim();
1989 const int coarse_ndofs = height/fine_vdim;
1990 for (
int coarse_el = 0; coarse_el < coarse_to_fine.Size(); coarse_el++)
1992 coarse_elem_dof->
GetRow(coarse_el, c_vdofs);
1993 fine_fes->DofsToVDofs(c_vdofs, coarse_ndofs);
1998 const int ref_type = coarse_to_ref_type[coarse_el];
2000 const int *fine_elems = coarse_to_fine.GetRow(coarse_el);
2001 const int num_fine_elems = coarse_to_fine.RowSize(coarse_el);
2002 const int lR_offset = ref_type_to_fine_elem_offset[ref_type];
2003 for (
int s = 0;
s < num_fine_elems;
s++)
2006 fine_fes->GetElementVDofs(fine_elems[
s], f_vdofs);
2010 AddMult(lR, loc_x_mat, loc_y_mat);
2025 const int nmat = pmats.
SizeK();
2026 const int ldof = fe->
GetDof();
2032 localR.
SetSize(ldof, ldof, nmat);
2033 for (
int i = 0; i < nmat; i++)
2041 const Table* old_elem_dof,
2042 const Table* old_elem_fos)
2046 MFEM_VERIFY(
Nonconforming(),
"Not implemented for conforming meshes.");
2047 MFEM_VERIFY(old_ndofs,
"Missing previous (finer) space.");
2048 MFEM_VERIFY(
ndofs <= old_ndofs,
"Previous space is not finer.");
2056 for (
int i = 0; i < elem_geoms.
Size(); i++)
2064 localR[elem_geoms[0]].SizeI());
2072 MFEM_ASSERT(dtrans.
embeddings.Size() == old_elem_dof->
Size(),
"");
2075 for (
int k = 0; k < dtrans.
embeddings.Size(); k++)
2082 old_elem_dof->
GetRow(k, old_dofs);
2084 for (
int vd = 0; vd <
vdim; vd++)
2086 old_dofs.
Copy(old_vdofs);
2089 for (
int i = 0; i < lR.
Height(); i++)
2091 if (!std::isfinite(lR(i, 0))) {
continue; }
2094 int m = (r >= 0) ? r : (-1 - r);
2099 R->
SetRow(r, old_vdofs, row);
2107 MFEM_VERIFY(num_marked == R->
Height(),
2108 "internal error: not all rows of R were set.");
2127 int nmat = pmats.
SizeK();
2134 for (
int i = 0; i < nmat; i++)
2143 int vdim_,
int ordering_)
2164 MFEM_VERIFY(mesh_->
NURBSext,
"NURBS FE space requires a NURBS mesh.");
2166 if (NURBSext_ == NULL)
2198 for (
int i=0; i<
DoFTrans.Size(); i++)
2203 if (dynamic_cast<const ND_FECollection*>(
fec))
2235 mfem_error(
"FiniteElementSpace::StealNURBSext");
2244 MFEM_VERIFY(
NURBSext,
"NURBSExt not defined.");
2288 if (b == -1) {
continue; }
2298 for (
int i = 0; i < nv; i++)
2300 MFEM_VERIFY(fv[i] == bv[i],
2301 "non-matching face and boundary elements detected!");
2306 for (
int i = 0; i < row.
Size(); i++)
2309 face_dof_list.
Append(conn);
2318 MFEM_VERIFY(!
NURBSext,
"internal error");
2323 "Variable order space requires a nonconforming mesh.");
2345 "Mesh was not correctly finalized.");
2356 else if (mixed_faces)
2443 MFEM_ASSERT(bits != 0,
"invalid bit mask");
2444 for (
int order = 0; bits != 0; order++, bits >>= 1)
2446 if (bits & 1) {
return order; }
2473 for (
int j = 0; j < E.
Size(); j++)
2475 edge_orders[E[j]] |= mask;
2481 for (
int j = 0; j < F.
Size(); j++)
2483 face_orders[F[j]] |= mask;
2511 slave_orders |= edge_orders[edge_list.
slaves[i].index];
2514 int min_order =
MinOrder(slave_orders);
2530 if (slave.
index >= 0)
2532 slave_orders |= face_orders[slave.
index];
2535 for (
int j = 0; j < E.
Size(); j++)
2537 slave_orders |= edge_orders[E[j]];
2543 slave_orders |= edge_orders[-1 - slave.
index];
2547 int min_order =
MinOrder(slave_orders);
2559 for (
int j = 0; j < E.
Size(); j++)
2561 edge_orders[E[j]] |= face_orders[i];
2583 int num_ent = entity_orders.
Size();
2592 var_ent_order->
Reserve(num_ent);
2596 for (
int i = 0; i < num_ent; i++)
2601 for (
int order = 0; bits != 0; order++, bits >>= 1)
2609 if (var_ent_order) { var_ent_order->
Append(order); }
2624 int row,
int ndof)
const
2626 const int *beg = var_dof_table.
GetRow(row);
2627 const int *end = var_dof_table.
GetRow(row + 1);
2632 if ((beg[1] - beg[0]) == ndof) {
return beg[0]; }
2636 MFEM_ABORT(
"DOFs not found for ndof = " << ndof);
2646 if (variant >= end - beg) {
return -1; }
2662 if (variant >= end - beg) {
return -1; }
2672 MFEM_ASSERT(index >= 0 && index < dof_table.
Size(),
"");
2673 return dof_table.
GetRow(index + 1) - dof_table.
GetRow(index);
2676 static const char* msg_orders_changed =
2677 "Element orders changed, you need to Update() the space first.";
2714 for (
int i = 0; i < F.
Size(); i++)
2721 -> SetFaceOrientations(Fo);
2730 for (
int i = 0; i < V.
Size(); i++)
2732 for (
int j = 0; j < nv; j++)
2734 dofs.
Append(V[i]*nv + j);
2741 for (
int i = 0; i < E.
Size(); i++)
2746 for (
int j = 0; j < ne; j++)
2755 for (
int i = 0; i < F.
Size(); i++)
2763 for (
int j = 0; j < nf; j++)
2775 for (
int j = 0; j < nb; j++)
2785 if (i < 0 || !mesh->
GetNE()) {
return NULL; }
2786 MFEM_VERIFY(i < mesh->
GetNE(),
2787 "Invalid element id " << i <<
", maximum allowed " <<
mesh->
GetNE()-1);
2804 "internal error: " <<
2827 SetFaceOrientations(Fo);
2860 SetFaceOrientations(Fo);
2869 for (
int i = 0; i < V.
Size(); i++)
2871 for (
int j = 0; j < nv; j++)
2873 dofs.
Append(V[i]*nv + j);
2880 for (
int i = 0; i < E.
Size(); i++)
2885 for (
int j = 0; j < ne; j++)
2897 for (
int j = 0; j < nf; j++)
2920 int order, nf, fbase;
2928 if (variant >= end - beg) {
return -1; }
2930 fbase = beg[variant];
2931 nf = beg[variant+1] - fbase;
2939 if (variant > 0) {
return -1; }
2958 for (
int i = 0; i < V.
Size(); i++)
2960 for (
int j = 0; j < nv; j++)
2962 dofs.
Append(V[i]*nv + j);
2968 for (
int i = 0; i < E.
Size(); i++)
2973 for (
int j = 0; j < ne; j++)
2979 for (
int j = 0; j < nf; j++)
2992 int order, ne, base;
2997 if (variant >= end - beg) {
return -1; }
2999 base = beg[variant];
3000 ne = beg[variant+1] - base;
3007 if (variant > 0) {
return -1; }
3020 for (
int i = 0; i < 2; i++)
3022 for (
int j = 0; j < nv; j++)
3024 dofs.
Append(V[i]*nv + j);
3027 for (
int j = 0; j < ne; j++)
3039 for (
int j = 0; j < nv; j++)
3054 for (
int j = 0; j < nb; j++)
3072 for (
int j = 0, k =
nvdofs+i*ne; j < ne; j++, k++)
3096 for (
int j = 0; j < nf; j++)
3158 "NURBS mesh: only boundary faces are supported!");
3168 MFEM_ASSERT(
mesh->
Dimension() > 1,
"No edges with mesh dimension < 2");
3193 for (
int i = 0; i <
E2Q_array.Size(); i++)
3234 ceed::RemoveBasisAndRestriction(
this);
3239 for (
int i = 0; i <
DoFTrans.Size(); i++)
3256 for (
int i = 0; i < elem_geoms.
Size(); i++)
3259 localP[elem_geoms[i]]);
3295 if (RP_case == 0) {
return; }
3308 cR, T.
Ptr(), coarse_P,
false, owner,
false));
3330 MFEM_ABORT(
"not implemented yet");
3346 MFEM_ABORT(
"Error in update sequence. Space needs to be updated after "
3347 "each mesh modification.");
3354 MFEM_ABORT(
"Updating space after both mesh change and element order "
3355 "change is not supported. Please update separately after "
3366 Table* old_elem_dof = NULL;
3367 Table* old_elem_fos = NULL;
3393 MFEM_VERIFY(!old_orders_changed,
"Interpolation for element order change "
3394 "is not implemented yet, sorry.");
3404 old_elem_fos, old_ndofs));
3407 old_elem_dof = NULL;
3408 old_elem_fos = NULL;
3427 false,
false,
true));
3436 delete old_elem_dof;
3437 delete old_elem_fos;
3448 int fes_format = 90;
3449 bool nurbs_unit_weights =
false;
3460 MFEM_VERIFY(nurbs_fec,
"invalid FE collection");
3462 const double eps = 5e-14;
3473 os << (fes_format == 90 ?
3474 "FiniteElementSpace\n" :
"MFEM FiniteElementSpace v1.0\n")
3475 <<
"FiniteElementCollection: " <<
fec->
Name() <<
'\n'
3476 <<
"VDim: " <<
vdim <<
'\n'
3477 <<
"Ordering: " <<
ordering <<
'\n';
3479 if (fes_format == 100)
3493 os <<
"NURBS_orders\n";
3500 os <<
"NURBS_periodic\n";
3505 if (!nurbs_unit_weights)
3507 os <<
"NURBS_weights\n";
3511 os <<
"End: MFEM FiniteElementSpace v1.0\n";
3518 int fes_format = 0, ord;
3524 getline(input, buff);
3526 if (buff ==
"FiniteElementSpace") { fes_format = 90; }
3527 else if (buff ==
"MFEM FiniteElementSpace v1.0") { fes_format = 100; }
3528 else { MFEM_ABORT(
"input stream is not a FiniteElementSpace!"); }
3529 getline(input, buff,
' ');
3531 getline(input, buff);
3534 getline(input, buff,
' ');
3536 getline(input, buff,
' ');
3541 if (fes_format == 90)
3545 MFEM_VERIFY(m->
NURBSext,
"NURBS FE collection requires a NURBS mesh!");
3546 const int order = nurbs_fec->
GetOrder();
3554 else if (fes_format == 100)
3559 MFEM_VERIFY(input.good(),
"error reading FiniteElementSpace v1.0");
3560 getline(input, buff);
3562 if (buff ==
"NURBS_order" || buff ==
"NURBS_orders")
3564 MFEM_VERIFY(nurbs_fec,
3565 buff <<
": NURBS FE collection is required!");
3566 MFEM_VERIFY(m->
NURBSext, buff <<
": NURBS mesh is required!");
3567 MFEM_VERIFY(!nurbs_ext, buff <<
": order redefinition!");
3568 if (buff ==
"NURBS_order")
3581 else if (buff ==
"NURBS_periodic")
3588 else if (buff ==
"NURBS_weights")
3590 MFEM_VERIFY(nurbs_ext,
"NURBS_weights: NURBS_orders have to be "
3591 "specified before NURBS_weights!");
3594 else if (buff ==
"element_orders")
3596 MFEM_VERIFY(!nurbs_fec,
"section element_orders cannot be used "
3597 "with a NURBS FE collection");
3598 MFEM_ABORT(
"element_orders: not implemented yet!");
3600 else if (buff ==
"End: MFEM FiniteElementSpace v1.0")
3606 MFEM_ABORT(
"unknown section: " << buff);
3622 element_offsets =
new int[num_elem + 1];
3627 for (
int i = 0; i < num_elem; i++)
3629 element_offsets[i] = offset;
3631 if (int_rule[geom] == NULL)
3637 element_offsets[num_elem] = size = offset;
3643 const char *msg =
"invalid input stream";
3646 in >> ident; MFEM_VERIFY(ident ==
"QuadratureSpace", msg);
3647 in >> ident; MFEM_VERIFY(ident ==
"Type:", msg);
3649 if (ident ==
"default_quadrature")
3651 in >> ident; MFEM_VERIFY(ident ==
"Order:", msg);
3656 MFEM_ABORT(
"unknown QuadratureSpace type: " << ident);
3665 os <<
"QuadratureSpace\n"
3666 <<
"Type: default_quadrature\n"
3667 <<
"Order: " <<
order <<
'\n';
int GetNPoints() const
Returns the number of the points in the integration rule.
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof, const Table *old_elem_fos)
Abstract class for all finite elements.
int GetOrder() const
Get the order of the NURBS collection: either a positive number, when using fixed order...
VDofTransformation VDoFTrans
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Field is discontinuous across element interfaces.
Ordering::Type GetOrdering() const
Return the ordering method.
int Size() const
Return the logical size of the array.
int RowSize(const int i) const
Returns the number of elements in row i.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
void GetVertexVDofs(int i, Array< int > &vdofs) const
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
const Vector & GetWeights() const
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
void Load(std::istream &in, int fmt=0)
Read an Array from the stream in using format fmt. The format fmt can be:
int GetNDofs() const
Returns number of degrees of freedom.
const FiniteElement * GetFE(Geometry::Type geom, int p) const
Variable order version of FiniteElementForGeometry().
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
Class for an integration rule - an Array of IntegrationPoint.
int * bdofs
internal DOFs of elements if mixed/var-order; NULL otherwise
const Array< int > & GetMaster() const
void GetElementInteriorVDofs(int i, Array< int > &vdofs) const
int DofToVDof(int dof, int vd, int ndofs=-1) const
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I, int skipfirst=0)
const CoarseFineTransformations & GetDerefinementTransforms()
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
void SetSize(int dim, int connections_per_row)
Set the size and the number of connections for the table.
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
void LoadBE(int i, const FiniteElement *BE) const
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
SparseMatrix * DerefinementMatrix(int old_ndofs, const Table *old_elem_dof, const Table *old_elem_fos)
Calculate GridFunction restriction matrix after mesh derefinement.
Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
void BuildElementToDofTable() const
void AddColumnsInRow(int r, int ncol)
int MakeDofTable(int ent_dim, const Array< int > &entity_orders, Table &entity_dofs, Array< char > *var_ent_order)
void BuildNURBSFaceToDofTable() const
Generates partial face_dof table for a NURBS space.
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
DerefinementOperator(const FiniteElementSpace *f_fes, const FiniteElementSpace *c_fes, BilinearFormIntegrator *mass_integ)
TODO: Implement DofTransformation support.
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
void SetSize(int s)
Resize the vector to size s.
int GetNBE() const
Returns number of boundary elements.
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'.
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
static int DecodeDof(int dof)
Helpers to remove encoded sign from a DOF.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Lists all edges/faces in the nonconforming mesh.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
static int MinOrder(VarOrderBits bits)
Return the minimum order (least significant bit set) in the bit mask.
Pointer to an Operator of a specified type.
Operator::Type Type() const
Get the currently set operator type id.
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
static constexpr int MaxVarOrder
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
virtual ~DerefinementOperator()
unsigned matrix
index into NCList::point_matrices[geom]
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Data type dense matrix using column-major storage.
int vdim
Vector dimension (number of unknowns per degree of freedom).
Geometry::Type Geom() const
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
int GetNumElementInteriorDofs(int i) const
void Save(std::ostream &out, int fmt=0) const
Save the Array to the stream out using the format fmt. The format fmt can be:
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
int GetBdrAttribute(int i) const
Operator that extracts face degrees of freedom for L2 nonconforming spaces.
int GetBdrElementEdgeIndex(int i) const
const SparseMatrix * GetHpConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
int GetNE() const
Returns number of elements.
const Element * GetFace(int i) const
int GetEdgeDofs(int edge, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
T Sum()
Return the sum of all the array entries using the '+'' operator for class 'T'.
MFEM_DEPRECATED void RebuildElementToDofTable()
(
Geometry::Type GetFaceBaseGeometry(int i) const
bool operator<(const Pair< A, B > &p, const Pair< A, B > &q)
Comparison operator for class Pair, based on the first element only.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
virtual void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer true-dof data from coarse_fes, defined on a coarse mesh, to this FE space, defined on a refined mesh.
void AddEdgeFaceDependencies(SparseMatrix &deps, Array< int > &master_dofs, const FiniteElement *master_fe, Array< int > &slave_dofs, int slave_face, const DenseMatrix *pm) const
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
double * GetData() const
Return a pointer to the beginning of the Vector data.
virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes, const Array< int > *perm)
Copies the prolongation and restriction matrices from fes.
const FiniteElementCollection * fec
Associated FE collection (not owned).
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
virtual int GetFaceDofs(int face, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Geometry::Type GetElementBaseGeometry(int i) const
int Size_of_connections() const
void GetLocalDerefinementMatrices(Geometry::Type geom, DenseTensor &localR) const
Geometry::Type GetBdrElementGeometry(int i) const
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
void DeleteAll()
Delete the whole array.
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
void AddConnections(int r, const int *c, int nc)
std::uint64_t VarOrderBits
Bit-mask representing a set of orders needed by an edge/face.
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Array< int > dof_elem_array
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
virtual int DofForGeometry(Geometry::Type GeomType) const =0
virtual void UpdateMeshPointer(Mesh *new_mesh)
void ConvertFromConformingVDofs(const Array< int > &cdofs, Array< int > &dofs)
For a partially conforming FE space, convert a marker array (nonzero entries are true) on the conform...
const IntegrationRule * IntRule
Not owned.
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
int uni_fdof
of single face DOFs if all faces uniform; -1 otherwise
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
ID for class SparseMatrix.
const Array< int > & GetSlave() const
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
int GetNE() const
Returns number of elements in the mesh.
void Set(const int i, const int j, const double val)
void GetVertexDofs(int i, Array< int > &dofs) const
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
General product operator: x -> (A*B)(x) = A(B(x)).
void BuildBdrElementToDofTable() const
int HasFaceDofs(Geometry::Type geom, int p) const
ID for the base class Operator, i.e. any type.
int FindFaceDof(int face, int ndof) const
Similar to FindEdgeDof, but used for mixed meshes too.
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
int GetFaceOrder(int face, int variant=0) const
Returns the polynomial degree of the i'th face finite element.
double f(const Vector &xvec)
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Extract all column indices and values from a given row.
Array< DofTransformation * > DoFTrans
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
int GetNBE() const
Returns number of boundary elements in the mesh.
void MakeFromList(int nrows, const Array< Connection > &list)
bool orders_changed
True if at least one element order changed (variable-order space only).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int FindEdgeDof(int edge, int ndof) const
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void EnsureMultTranspose() const
Ensures that the matrix is capable of performing MultTranspose(), AddMultTranspose(), and AbsMultTranspose().
Mesh * GetMesh() const
Returns the mesh.
Operator that extracts Face degrees of freedom for L2 spaces.
void GetVDofs(int vd, Array< int > &dofs, int ndofs=-1) const
Returns the indices of all of the VDofs for the specified dimension 'vd'.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
SparseMatrix * cR
Conforming restriction matrix such that cR.cP=I.
void GetRow(int r, Vector &row) const
const IntegrationRule * IntRule
Not owned.
SparseMatrix * D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Generate the global restriction matrix from a discontinuous FE space to the piecewise constant FE spa...
static const int NumVerts[NumGeom]
SparseMatrix * D2C_GlobalRestrictionMatrix(FiniteElementSpace *cfes)
Generate the global restriction matrix from a discontinuous FE space to the continuous FE space of th...
void GetBoundaryTrueDofs(Array< int > &boundary_dofs, int component=-1)
Get a list of all boundary true dofs, boundary_dofs. For spaces with 'vdim' > 1, the 'component' param...
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void BuildConformingInterpolation() const
Calculate the cP and cR matrices for a nonconforming mesh.
Geometry::Type GetFaceGeometry(int i) const
NURBSExtension * StealNURBSext()
Array< DenseMatrix * > point_matrices[Geometry::NumGeom]
List of unique point matrices for each slave geometry.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
virtual ~RefinementOperator()
virtual const FiniteElement * TraceFiniteElementForGeometry(Geometry::Type GeomType) const
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i...
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Array< int > dof_ldof_array
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
const Array< int > & GetOrders() const
Read-only access to the orders of all knot vectors.
int GetVDim() const
Returns vector dimension.
int Size() const
Returns the number of TYPE I elements.
const int * GetDofOrdering(Geometry::Type geom, int p, int ori) const
Variable order version of DofOrderForOrientation().
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
static int EncodeDof(int entity_base, int idx)
Helper to encode a sign flip into a DOF index (for Hcurl/Hdiv shapes).
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Array< QuadratureInterpolator * > E2Q_array
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
int GetNumBorderDofs(Geometry::Type geom, int order) const
Operator * Ptr() const
Access the underlying Operator pointer.
void Add(const int i, const int j, const double val)
int GetEntityDofs(int entity, int index, Array< int > &dofs, Geometry::Type master_geom=Geometry::INVALID, int variant=0) const
Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
OperatorHandle L2E_nat
The element restriction operators, see GetElementRestriction().
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
int GetNumDof(Geometry::Type geom, int p) const
Variable order version of DofForGeometry().
int slaves_end
slave faces
void GetTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer GridFunction data from coarse_fes...
void GetElementInteriorDofs(int i, Array< int > &dofs) const
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
List of mesh geometries stored as Array<Geometry::Type>.
double Min() const
Returns the minimal element of the vector.
Type
Enumeration defining IDs for some classes derived from Operator.
void Swap(Array< T > &, Array< T > &)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double p(const Vector &x, double t)
void LoadFE(int i, const FiniteElement *FE) const
void SetOrder(int Order) const
Set the order and the name, based on the given Order: either a positive number for fixed order...
Operator that converts FiniteElementSpace L-vectors to E-vectors.
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Nonconforming edge/face within a bigger edge/face.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
virtual ~FiniteElementSpace()
bool OwnsOperator() const
Return true if the OperatorHandle owns the held Operator.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
SparseMatrix * H2L_GlobalRestrictionMatrix(FiniteElementSpace *lfes)
Construct the restriction matrix from the FE space given by (*this) to the lower degree FE space give...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
RefinementOperator(const FiniteElementSpace *fespace, Table *old_elem_dof, Table *old_elem_fos, int old_ndofs)
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Array< char > var_face_orders
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Helper struct for defining a connectivity table, see Table::MakeFromList.
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Geometry::Type GetElementGeometry(int i) const
Table * GetElementDofTable()
void MakeVDimMatrix(SparseMatrix &mat) const
Replicate 'mat' in the vector dimension, according to vdim ordering mode.
Table * GetBdrElementDofTable()
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
void GetEdgeInteriorVDofs(int i, Array< int > &vdofs) const
Linear2DFiniteElement TriangleFE
int GetNVariants(int entity, int index) const
Return number of possible DOF variants for edge/face (var. order spaces).
const CoarseFineTransformations & GetRefinementTransforms()
Operation GetLastOperation() const
Return type of last modification of the mesh.
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
bool Nonconforming() const
NURBSExtension * NURBSext
Optional NURBS mesh extension.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
int GetDegenerateFaceDofs(int index, Array< int > &dofs, Geometry::Type master_geom, int variant) const
void ConvertToConformingVDofs(const Array< int > &dofs, Array< int > &cdofs)
For a partially conforming FE space, convert a marker array (nonzero entries are true) on the partial...
virtual const char * Name() const
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Mesh * mesh
The mesh that FE space lives on (not owned).
void ReorderElementToDofTable()
Reorder the scalar DOFs based on the element ordering.
GridFunction interpolation operator applicable after mesh refinement.
void CalcEdgeFaceVarOrders(Array< VarOrderBits > &edge_orders, Array< VarOrderBits > &face_orders) const
double Max() const
Returns the maximal element of the vector.
const NCList & GetNCList(int entity)
Return vertex/edge/face list (entity = 0/1/2, respectively).
Array< FaceQuadratureInterpolator * > E2BFQ_array
const Table * GetElementToFaceOrientationTable() const
int GetNConformingDofs() const
Table var_face_dofs
NOTE: also used for spaces with mixed faces.
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
General triple product operator x -> A*B*C*x, with ownership of the factors.
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs. ...
Array< FaceQuadratureInterpolator * > E2IFQ_array
void BuildFaceToDofTable() const
void UpdateElementOrders()
Resize the elem_order array on mesh change.
void OrientedPointMatrix(const Slave &slave, DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int index(int i, int j, int nx, int ny)
Array< char > var_edge_orders
int GetEdgeOrder(int edge, int variant=0) const
Lexicographic ordering for tensor-product FiniteElements.
virtual int GetNVertices() const =0
int parent
Coarse Element index in the coarse mesh.
FiniteElementSpace()
Default constructor: the object is invalid until initialized using the method Load().
QuadratureSpace(Mesh *mesh_, int order_)
Create a QuadratureSpace based on the global rules from IntRules.
void Constructor(Mesh *mesh, NURBSExtension *ext, const FiniteElementCollection *fec, int vdim=1, int ordering=Ordering::byNODES)
Help function for constructors + Load().
NCMesh * ncmesh
Optional nonconforming mesh extension.
Geometry::Type GetBdrElementBaseGeometry(int i) const
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th face element (2D and 3D).
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th edge.
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
int GetNEdges() const
Return the number of edges.
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
MFEM_HOST_DEVICE void Set(const int height, const int width, const double alpha, const TA *Adata, TB *Bdata)
Compute B = alpha*A, where the matrices A and B are of size height x width with data Adata and Bdata...
void MakeRef(T *, int)
Make this Array a reference to a pointer.
int GetNFaces() const
Return the number of faces in a 3D mesh.
BiLinear2DFiniteElement QuadrilateralFE
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
const FiniteElementCollection * FEColl() const
void Mult(const double *x, double *y) const
Matrix vector multiplication.
std::tuple< bool, ElementDofOrdering, FaceType, L2FaceValues > key_face
The face restriction operators, see GetFaceRestriction().
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
SparseMatrix * cR_hp
A version of the conforming restriction matrix for variable-order spaces.
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Class representing the storage layout of a QuadratureFunction.
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
NURBSExtension * NURBSext
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Base class for operators that extracts Face degrees of freedom.
void Swap(SparseMatrix &other)
int GetElementOrderImpl(int i) const
Return element order: internal version of GetElementOrder without checks.
Rank 3 tensor (array of matrices)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
void Save(std::ostream &out) const
Save finite element space to output stream out.
Abstract data type element.
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
Linear1DFiniteElement SegmentFE
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
int FindDofs(const Table &var_dof_table, int row, int ndof) const
Search row of a DOF table for a DOF set of size 'ndof', return first DOF.
static void AdjustVDofs(Array< int > &vdofs)
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
const Element * GetBdrElement(int i) const
Defines the position of a fine element within a coarse element.
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
const QuadratureSpace * qspace
Not owned.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...