45 istream::int_type next_char = input.peek();
51 if (buff ==
"NURBS_patches")
54 "NURBS_patches requires NURBS FE space");
59 MFEM_ABORT(
"unknown section: " << buff);
100 int g_nddofs = g_ndofs - (g_nvdofs + g_nedofs + g_nfdofs);
102 vi = ei = fi = di = 0;
103 for (
int i = 0; i < num_pieces; i++)
110 int l_nddofs = l_ndofs - (l_nvdofs + l_nedofs + l_nfdofs);
115 for (
int d = 0; d < vdim; d++)
117 memcpy(g_data+vi, l_data, l_nvdofs*
sizeof(
real_t));
120 memcpy(g_data+ei, l_data, l_nedofs*
sizeof(
real_t));
123 memcpy(g_data+fi, l_data, l_nfdofs*
sizeof(
real_t));
126 memcpy(g_data+di, l_data, l_nddofs*
sizeof(
real_t));
133 memcpy(g_data+vdim*vi, l_data, l_nvdofs*
sizeof(
real_t)*vdim);
134 l_data += vdim*l_nvdofs;
135 g_data += vdim*g_nvdofs;
136 memcpy(g_data+vdim*ei, l_data, l_nedofs*
sizeof(
real_t)*vdim);
137 l_data += vdim*l_nedofs;
138 g_data += vdim*g_nedofs;
139 memcpy(g_data+vdim*fi, l_data, l_nfdofs*
sizeof(
real_t)*vdim);
140 l_data += vdim*l_nfdofs;
141 g_data += vdim*g_nfdofs;
142 memcpy(g_data+vdim*di, l_data, l_nddofs*
sizeof(
real_t)*vdim);
143 l_data += vdim*l_nddofs;
144 g_data += vdim*g_nddofs;
182 old_data.
Swap(*
this);
185 T->
Mult(old_data, *
this);
213 MFEM_ASSERT(v.
Size() >= v_offset +
f->GetVSize(),
"");
245 MFEM_ASSERT(tv.
Size() >= tv_offset +
f->GetTrueVSize(),
"");
266 int nfe = ufes->
GetNE();
274 for (
int i = 0; i < nfe; i++)
276 if (subdomain >= 0 && ufes->
GetAttribute(i) != subdomain)
284 u.GetSubVector(udofs, ul);
292 *ffes->
GetFE(i), fl, wcoef);
301 for (
int j = 0; j < fdofs.
Size(); j++)
317 for (
int i = 0; i < count.
Size(); i++)
319 if (count[i] != 0) { flux(i) /= count[i]; }
403 int dof = FElem->
GetDof();
419 for (
int k = 0; k < n; k++)
422 nval[k] = shape * (&loc_data[dof * vdim]);
428 for (
int k = 0; k < n; k++)
432 nval[k] = shape * (&loc_data[dof * vdim]);
440 for (
int k = 0; k < n; k++)
444 nval[k] = loc_data * (&vshape(0,vdim));
473 return (DofVal * LocVec);
480 int dof = FElem->
GetDof();
504 for (
int k = 0; k < vdim; k++)
506 val(k) = shape * (&loc_data[dof * k]);
531 int dof = FElem->
GetDof();
532 Vector DofVal(dof), loc_data(dof);
540 for (
int k = 0; k < n; k++)
543 vals(k) = DofVal * loc_data;
549 for (
int k = 0; k < n; k++)
553 vals(k) = DofVal * loc_data;
582 "invalid FE map type");
584 int dof = FElem->
GetDof();
585 Vector DofLap(dof), loc_data(dof);
587 for (
int k = 0; k < n; k++)
592 laps(k) = DofLap * loc_data;
625 "invalid FE map type");
627 int dof = FElem->
GetDof();
635 for (
int k = 0; k < n; k++)
641 for (
int j = 0; j <
size; j++)
643 for (
int d = 0; d < dof; d++)
645 hess(k,j) += DofHes(d,j) * loc_data[d];
725 int comp,
Vector *tr)
const
751 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \""
753 <<
"on mesh edges.");
766 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \""
768 <<
"on mesh faces.");
787 MFEM_ASSERT(FET !=
nullptr,
788 "FaceElementTransformation must be valid for a boundary element");
817 MFEM_ABORT(
"GridFunction::GetValue: Unsupported element type \""
835 return (DofVal * LocVec);
850 for (
int j = 0; j < nip; j++)
887 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \""
889 <<
"on mesh edges.");
902 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \""
904 <<
"on mesh faces.");
923 MFEM_ASSERT(FET !=
nullptr,
924 "FaceElementTransformation must be valid for a boundary element");
944 MFEM_ASSERT(FET !=
nullptr,
945 "FaceElementTransformation must be valid for a boundary element");
955 MFEM_ABORT(
"GridFunction::GetVectorValue: Unsupported element type \""
957 if (val.
Size() > 0) { val = NAN; }
982 for (
int k = 0; k < vdim; k++)
984 val(k) = shape * (&loc_data[dof * k]);
1009 int dof = FElem->
GetDof();
1027 for (
int j = 0; j < nip; j++)
1033 for (
int k = 0; k < vdim; k++)
1035 vals(k,j) = shape * (&loc_data[dof * k]);
1042 int vdim = std::max(spaceDim, FElem->
GetRangeDim());
1048 for (
int j = 0; j < nip; j++)
1089 MFEM_ASSERT(Transf !=
nullptr,
"FaceElementTransformation cannot be null!");
1096 MFEM_ASSERT(Transf !=
nullptr,
"FaceElementTransformation cannot be null!");
1112 Vector shape, loc_values, orig_loc_values;
1113 int i, j, d, ne, dof, odof, vdim;
1117 for (i = 0; i < ne; i++)
1129 odof = orig_fe->
GetDof();
1130 loc_values.
SetSize(dof * vdim);
1133 for (j = 0; j < dof; j++)
1137 for (d = 0; d < vdim; d++)
1139 loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
1158 Vector shape, loc_values, loc_values_t, orig_loc_values, orig_loc_values_t;
1159 int i, j, d, nbe, dof, odof, vdim;
1163 for (i = 0; i < nbe; i++)
1171 odof = orig_fe->
GetDof();
1172 loc_values.
SetSize(dof * vdim);
1175 for (j = 0; j < dof; j++)
1179 for (d = 0; d < vdim; d++)
1181 loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
1195 int d, k, n, sdim, dof;
1207 Vector loc_data, val(sdim);
1213 for (k = 0; k < n; k++)
1219 for (d = 0; d < sdim; d++)
1239 for (j = 0; j < ndofs; j++)
1240 for (i = 0; i < vdim; i++)
1242 temp[j+i*ndofs] =
data[k++];
1245 for (i = 0; i <
size; i++)
1273 val(vertices[k]) += vals(k, comp);
1274 overlap[vertices[k]]++;
1278 for (i = 0; i < overlap.
Size(); i++)
1280 val(i) /= overlap[i];
1288 int d, i, k, ind, dof, sdim;
1297 for (i = 0; i < new_fes->
GetNE(); i++)
1304 for (d = 0; d < sdim; d++)
1306 for (k = 0; k < dof; k++)
1308 if ( (ind=new_vdofs[dof*d+k]) < 0 )
1310 ind = -1-ind, vals(k, d) = - vals(k, d);
1312 vec_field(ind) += vals(k, d);
1318 for (i = 0; i < overlap.
Size(); i++)
1320 vec_field(i) /= overlap[i];
1333 Vector pt_grad, loc_func;
1334 int i, j, k,
dim, dof, der_dof, ind;
1341 for (i = 0; i < der_fes->
GetNE(); i++)
1350 der_dof = der_fe->
GetDof();
1356 for (j = 0; j < dof; j++)
1357 loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
1359 for (k = 0; k < der_dof; k++)
1367 for (j = 0; j <
dim; j++)
1369 a += inv_jac(j, der_comp) * pt_grad(j);
1371 der(der_dofs[k]) +=
a;
1372 zones_per_dof[der_dofs[k]]++;
1382 for (
int i = 0; i < overlap.
Size(); i++)
1384 der(i) /= overlap[i];
1401 MultAtB(loc_data_mat, dshape, gh);
1415 "invalid FE map type");
1420 for (
int i = 0; i < Jinv.
Width(); i++)
1422 for (
int j = 0; j < Jinv.
Height(); j++)
1424 div_v += grad_hat(i, j) * Jinv(j, i);
1441 return (loc_data * divshape) / T.
Weight();
1483 MFEM_ABORT(
"GridFunction::GetDivergence: Unsupported element type \""
1501 "invalid FE map type");
1507 Mult(grad_hat, Jinv, grad);
1512 curl(0) = grad(2,1) - grad(1,2);
1513 curl(1) = grad(0,2) - grad(2,0);
1514 curl(2) = grad(1,0) - grad(0,1);
1516 else if (grad.
Height() == 2)
1519 curl(0) = grad(1,0) - grad(0,1);
1576 MFEM_ABORT(
"GridFunction::GetCurl: Unsupported element type \""
1590 "invalid FE map type");
1591 MFEM_ASSERT(
fes->
GetVDim() == 1,
"Defined for scalar functions.");
1641 MFEM_ABORT(
"GridFunction::GetGradient: Unsupported element type \""
1651 int elNo = tr.ElementNo;
1664 tr.SetIntPoint(&ip);
1684 Mult(grad_hat, Jinv, grad);
1723 MFEM_ABORT(
"GridFunction::GetVectorGradient: "
1724 "Unsupported element type \"" << T.
ElementType <<
"\"");
1736 Vector loc_avgs, loc_this;
1741 for (
int i = 0; i <
fes->
GetNE(); i++)
1753 loc_mass.
Mult(loc_this, loc_avgs);
1760 loc_mass.
Mult(loc_this, loc_avgs);
1763 for (
int i = 0; i < avgs.
Size(); i++)
1765 avgs(i) /= int_psi(i);
1786 if (!mesh->
GetNE()) {
return; }
1797 MFEM_VERIFY(vdim == src.
fes->
GetVDim(),
"incompatible vector dimensions!");
1802 for (
int i = 0; i < mesh->
GetNE(); i++)
1819 for (
int vd = 0; vd < vdim; vd++)
1846 MFEM_ASSERT(weights.
Size() ==
size,
"Different # of weights and dofs.");
1847 MFEM_ASSERT(lo_.
Size() ==
size,
"Different # of lower bounds and dofs.");
1848 MFEM_ASSERT(hi_.
Size() ==
size,
"Different # of upper bounds and dofs.");
1859 slbqp.
Mult(vals, new_vals);
1884 if (max_val <= min_)
1895 if (min_ <= min_val && max_val <= max_)
1901 minv = (min_ > min_val) ? min_ : min_val;
1902 maxv = (max_ < max_val) ? max_ : max_val;
1915 R->
Mult(*
this, tmp);
1916 P->
Mult(tmp, *
this);
1934 for (j = 0; j < vertices.
Size(); j++)
1936 nval(vertices[j]) += values[j];
1937 overlap[vertices[j]]++;
1940 for (i = 0; i < overlap.
Size(); i++)
1942 nval(i) /= overlap[i];
1953 for (
int i = 0; i <
fes->
GetNE(); i++)
1957 for (
int j = 0; j < vdofs.
Size(); j++)
1959 elem_per_vdof[vdofs[j]]++;
1978 for (
int i = 0; i <
fes->
GetNE(); i++)
1986 for (
int j = 0; j < vdofs.
Size(); j++)
1990 MFEM_VERIFY(vals[j] != 0.0,
1991 "Coefficient has zeros, harmonic avg is undefined!");
1992 (*this)(vdofs[j]) += 1.0 / vals[j];
1996 (*this)(vdofs[j]) += vals[j];
1998 else { MFEM_ABORT(
"Not implemented"); }
2000 zones_per_vdof[vdofs[j]]++;
2019 for (
int i = 0; i <
fes->
GetNE(); i++)
2027 for (
int j = 0; j < vdofs.
Size(); j++)
2044 MFEM_VERIFY(vals[j] != 0.0,
2045 "Coefficient has zeros, harmonic avg is undefined!");
2046 (*this)(ldof) += isign / vals[j];
2050 (*this)(ldof) += isign*vals[j];
2053 else { MFEM_ABORT(
"Not implemented"); }
2055 zones_per_vdof[ldof]++;
2064 int i, j, fdof, d, ind, vdim;
2088 for (j = 0; j < fdof; j++)
2092 if (vcoeff) { vcoeff->
Eval(vc, *transf, ip); }
2093 for (d = 0; d < vdim; d++)
2095 if (!vcoeff && !coeff[d]) {
continue; }
2097 val = vcoeff ? vc(d) : coeff[d]->
Eval(*transf, ip);
2098 if ( (ind = vdofs[fdof*d+j]) < 0 )
2100 val = -val, ind = -1-ind;
2102 if (++values_counter[ind] == 1)
2108 (*this)(ind) += val;
2128 Array<int> bdr_edges, bdr_vertices, bdr_faces;
2131 for (i = 0; i < bdr_edges.
Size(); i++)
2133 int edge = bdr_edges[i];
2135 if (vdofs.
Size() == 0) {
continue; }
2143 for (d = 0; d < vdim; d++)
2145 if (!coeff[d]) {
continue; }
2147 fe->
Project(*coeff[d], *transf, vals);
2148 for (
int k = 0; k < vals.
Size(); k++)
2150 ind = vdofs[d*vals.
Size()+k];
2151 if (++values_counter[ind] == 1)
2153 (*this)(ind) = vals(k);
2157 (*this)(ind) += vals(k);
2165 fe->
Project(*vcoeff, *transf, vals);
2166 for (
int k = 0; k < vals.
Size(); k++)
2169 if (++values_counter[ind] == 1)
2171 (*this)(ind) = vals(k);
2175 (*this)(ind) += vals(k);
2186 for (
int i = 0; i < dofs.
Size(); i++)
2190 if (k < 0) { k = -1 - k; val = -val; }
2191 if (++values_counter[k] == 1)
2226 fe->
Project(vcoeff, *T, lvec);
2228 accumulate_dofs(dofs, lvec, *
this, values_counter);
2235 Array<int> bdr_edges, bdr_vertices, bdr_faces;
2238 for (
int i = 0; i < bdr_edges.
Size(); i++)
2240 int edge = bdr_edges[i];
2242 if (dofs.
Size() == 0) {
continue; }
2248 fe->
Project(vcoeff, *T, lvec);
2249 accumulate_dofs(dofs, lvec, *
this, values_counter);
2259 for (
int i = 0; i <
size; i++)
2261 const int nz = zones_per_vdof[i];
2262 if (nz) { (*this)(i) /= nz; }
2267 for (
int i = 0; i <
size; i++)
2269 const int nz = zones_per_vdof[i];
2270 if (nz) { (*this)(i) = nz/(*
this)(i); }
2275 MFEM_ABORT(
"invalid AvgType");
2297 for (
int i = 0; i < mesh->
GetNV(); i++)
2301 if (dist < min_dist)
2311 if (min_dist >= delta_coeff.
Tol())
2320 Vector vals, loc_mass_vals;
2321 for (
int i = 0; i < mesh->
GetNE(); i++)
2324 for (
int j = 0; j < vertices.
Size(); j++)
2325 if (vertices[j] == v_idx)
2339 loc_mass.
Mult(vals, loc_mass_vals);
2340 integral += loc_mass_vals.
Sum();
2351 if (delta_c == NULL)
2356 for (
int i = 0; i <
fes->
GetNE(); i++)
2374 (*this) *= (delta_c->
Scale() / integral);
2387 for (
int i = 0; i < dofs.
Size(); i++)
2400 (*this)(vdof) = coeff.
Eval(*T, ip);
2436 for (
int i = 0; i < dofs.
Size(); i++)
2448 vcoeff.
Eval(val, *T, ip);
2452 (*this)(vdof) = val(vd);
2485 int i, j, fdof, d, ind, vdim;
2501 for (j = 0; j < fdof; j++)
2505 for (d = 0; d < vdim; d++)
2507 if (!coeff[d]) {
continue; }
2509 val = coeff[d]->
Eval(*transf, ip);
2510 if ( (ind = vdofs[fdof*d+j]) < 0 )
2512 val = -val, ind = -1-ind;
2532 for (
int i = 0; i <
fes->
GetNE(); i++)
2541 for (
int j = 0; j < vdofs.
Size(); j++)
2543 if (attr > dof_attr[vdofs[j]])
2545 (*this)(vdofs[j]) = vals[j];
2546 dof_attr[vdofs[j]] = attr;
2588 for (
int i = 0; i < values_counter.
Size(); i++)
2590 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2606 ess_vdofs_marker = 0;
2610 if (!coeff[i]) {
continue; }
2612 for (
int j = 0; j<
Size(); j++)
2614 ess_vdofs_marker[j] = bool(ess_vdofs_marker[j]) ||
2615 bool(component_dof_marker[j]);
2618 for (
int i = 0; i < values_counter.
Size(); i++)
2620 MFEM_ASSERT(
bool(values_counter[i]) == ess_vdofs_marker[i],
2656 vcoeff.
Eval(vc, *T, ip);
2659 lvec.
Add(ip.
weight * (vc * nor), shape);
2687 vcoeff.
Eval(vc, *T, ip);
2689 lvec(j) = (vc * nor);
2710 for (
int i = 0; i < values_counter.
Size(); i++)
2712 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2727 int fdof, d, i, intorder, j, k;
2731 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2755 for (k = 0; k < fdof; k++)
2756 if (vdofs[fdof*d+k] >= 0)
2758 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2762 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2764 a -= exsol[d]->
Eval(*transf, ip);
2770 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2783 for (
int i = 0; i <
fes->
GetNE(); i++)
2785 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2787 int intorder = 2*fe->
GetOrder() + 3;
2799 exsol.
Eval(exact_vals, *T, *ir);
2802 vals.
Norm2(loc_errs);
2807 error += ip.
weight * T->
Weight() * (loc_errs(j) * loc_errs(j));
2811 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2845 exgrad->
Eval(vec,*Tr,ip);
2849 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2864 for (
int i = 0; i <
fes->
GetNE(); i++)
2884 exgrad->
Eval(vec,*Tr,ip);
2889 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2904 for (
int i = 0; i <
fes->
GetNE(); i++)
2924 excurl->
Eval(vec,*Tr,ip);
2930 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2942 for (
int i = 0; i <
fes->
GetNE(); i++)
2966 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2974 int fdof, intorder, k;
2979 Vector shape, el_dofs, err_val, ell_coeff_val;
3001 intorder = 2 * intorder;
3015 transf = face_elem_transf->
Elem1;
3021 for (k = 0; k < fdof; k++)
3024 el_dofs(k) = (*this)(vdofs[k]);
3028 el_dofs(k) = - (*this)(-1-vdofs[k]);
3035 ell_coeff_val(j) = ell_coeff->
Eval(*transf, eip);
3036 err_val(j) = exsol->
Eval(*transf, eip) - (shape * el_dofs);
3042 transf = face_elem_transf->
Elem2;
3048 for (k = 0; k < fdof; k++)
3051 el_dofs(k) = (*this)(vdofs[k]);
3055 el_dofs(k) = - (*this)(-1-vdofs[k]);
3062 ell_coeff_val(j) += ell_coeff->
Eval(*transf, eip);
3063 ell_coeff_val(j) *= 0.5;
3064 err_val(j) -= (exsol->
Eval(*transf, eip) - (shape * el_dofs));
3068 transf = face_elem_transf;
3074 error += (ip.
weight * nu * ell_coeff_val(j) *
3076 err_val(j) * err_val(j));
3080 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
3095 int norm_type)
const
3106 return sqrt(error1 * error1 + error2 * error2);
3115 return sqrt(L2error*L2error + GradError*GradError);
3124 return sqrt(L2error*L2error + DivError*DivError);
3133 return sqrt(L2error*L2error + CurlError*CurlError);
3144 int fdof, d, i, intorder, j, k;
3171 for (k = 0; k < fdof; k++)
3172 if (vdofs[fdof*d+k] >= 0)
3174 a += (*this)(vdofs[fdof*d+k]) * shape(k);
3178 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
3180 a -= exsol[d]->
Eval(*transf, ip);
3197 int i, fdof,
dim, intorder, j, k;
3201 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
3213 for (i = 0; i < mesh->
GetNE(); i++)
3215 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3232 for (k = 0; k < fdof; k++)
3235 el_dofs(k) = (*this)(vdofs[k]);
3239 el_dofs(k) = -(*this)(-1-vdofs[k]);
3246 a = (el_dofs * shape) - (exsol->
Eval(*transf, ip));
3252 for (i = 0; i < mesh->
GetNE(); i++)
3254 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3272 for (k = 0; k < fdof; k++)
3275 el_dofs(k) = (*this)(vdofs[k]);
3279 el_dofs(k) = -(*this)(-1-vdofs[k]);
3286 exgrad->
Eval(e_grad, *transf, ip);
3288 Mult(dshape, Jinv, dshapet);
3308 for (
int i = 0; i <
fes->
GetNE(); i++)
3310 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3319 int intorder = 2*fe->
GetOrder() + 3;
3328 real_t diff = fabs(vals(j) - exsol.
Eval(*T, ip));
3331 diff = pow(diff,
p);
3334 diff *= weight->
Eval(*T, ip);
3342 diff *= weight->
Eval(*T, ip);
3344 error = std::max(error, diff);
3354 error = -pow(-error, 1./
p);
3358 error = pow(error, 1./
p);
3371 "Incorrect size for result vector");
3378 for (
int i = 0; i <
fes->
GetNE(); i++)
3388 int intorder = 2*fe->
GetOrder() + 3;
3397 real_t diff = fabs(vals(j) - exsol.
Eval(*T, ip));
3400 diff = pow(diff,
p);
3403 diff *= weight->
Eval(*T, ip);
3411 diff *= weight->
Eval(*T, ip);
3413 error[i] = std::max(error[i], diff);
3421 error[i] = -pow(-error[i], 1./
p);
3425 error[i] = pow(error[i], 1./
p);
3442 for (
int i = 0; i <
fes->
GetNE(); i++)
3452 int intorder = 2*fe->
GetOrder() + 3;
3457 exsol.
Eval(exact_vals, *T, *ir);
3464 vals.
Norm2(loc_errs);
3468 v_weight->
Eval(exact_vals, *T, *ir);
3471 for (
int j = 0; j < vals.
Width(); j++)
3474 for (
int d = 0; d < vals.
Height(); d++)
3476 errj += vals(d,j)*exact_vals(d,j);
3478 loc_errs(j) = fabs(errj);
3485 real_t errj = loc_errs(j);
3488 errj = pow(errj,
p);
3491 errj *= weight->
Eval(*T, ip);
3499 errj *= weight->
Eval(*T, ip);
3501 error = std::max(error, errj);
3511 error = -pow(-error, 1./
p);
3515 error = pow(error, 1./
p);
3530 "Incorrect size for result vector");
3538 for (
int i = 0; i <
fes->
GetNE(); i++)
3548 int intorder = 2*fe->
GetOrder() + 3;
3553 exsol.
Eval(exact_vals, *T, *ir);
3560 vals.
Norm2(loc_errs);
3564 v_weight->
Eval(exact_vals, *T, *ir);
3567 for (
int j = 0; j < vals.
Width(); j++)
3570 for (
int d = 0; d < vals.
Height(); d++)
3572 errj += vals(d,j)*exact_vals(d,j);
3574 loc_errs(j) = fabs(errj);
3581 real_t errj = loc_errs(j);
3584 errj = pow(errj,
p);
3587 errj *= weight->
Eval(*T, ip);
3595 errj *= weight->
Eval(*T, ip);
3597 error[i] = std::max(error[i], errj);
3605 error[i] = -pow(-error[i], 1./
p);
3609 error[i] = pow(error[i], 1./
p);
3636 os <<
"NURBS_patches\n";
3655 ofstream ofs(fname);
3656 ofs.precision(precision);
3660#ifdef MFEM_USE_ADIOS2
3662 const std::string& variable_name,
3665 os.
Save(*
this, variable_name, type);
3681 os <<
"SCALARS " << field_name <<
" double 1\n"
3682 <<
"LOOKUP_TABLE default\n";
3683 for (
int i = 0; i < mesh->
GetNE(); i++)
3690 for (
int j = 0; j < val.
Size(); j++)
3692 os << val(j) <<
'\n';
3696 else if ( (vec_dim == 2 || vec_dim == 3) && mesh->
SpaceDimension() > 1)
3699 os <<
"VECTORS " << field_name <<
" double\n";
3700 for (
int i = 0; i < mesh->
GetNE(); i++)
3709 for (
int j = 0; j < vval.
Width(); j++)
3711 os << vval(0, j) <<
' ' << vval(1, j) <<
' ';
3727 for (
int vd = 0; vd < vec_dim; vd++)
3729 os <<
"SCALARS " << field_name << vd <<
" double 1\n"
3730 <<
"LOOKUP_TABLE default\n";
3731 for (
int i = 0; i < mesh->
GetNE(); i++)
3738 for (
int j = 0; j < val.
Size(); j++)
3740 os << val(j) <<
'\n';
3751 real_t v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
3752 real_t v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
3753 real_t n[] = { v1[1] * v2[2] - v1[2] * v2[1],
3754 v1[2] * v2[0] - v1[0] * v2[2],
3755 v1[0] * v2[1] - v1[1] * v2[0]
3757 real_t rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
3758 n[0] *= rl; n[1] *= rl; n[2] *= rl;
3760 os <<
" facet normal " << n[0] <<
' ' << n[1] <<
' ' << n[2]
3762 <<
"\n vertex " << p1[0] <<
' ' << p1[1] <<
' ' << p1[2]
3763 <<
"\n vertex " << p2[0] <<
' ' << p2[1] <<
' ' << p2[2]
3764 <<
"\n vertex " << p3[0] <<
' ' << p3[1] <<
' ' << p3[2]
3765 <<
"\n endloop\n endfacet\n";
3783 os <<
"solid GridFunction\n";
3785 bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
3786 bbox[2][0] = bbox[2][1] = 0.0;
3787 for (i = 0; i < mesh->
GetNE(); i++)
3794 for (k = 0; k < RG.
Size()/n; k++)
3796 for (j = 0; j < n; j++)
3799 pts[j][0] = pointmat(0,l);
3800 pts[j][1] = pointmat(1,l);
3801 pts[j][2] = values(l);
3817 bbox[0][0] = pointmat(0,0);
3818 bbox[0][1] = pointmat(0,0);
3819 bbox[1][0] = pointmat(1,0);
3820 bbox[1][1] = pointmat(1,0);
3821 bbox[2][0] = values(0);
3822 bbox[2][1] = values(0);
3825 for (j = 0; j < values.
Size(); j++)
3827 if (bbox[0][0] > pointmat(0,j))
3829 bbox[0][0] = pointmat(0,j);
3831 if (bbox[0][1] < pointmat(0,j))
3833 bbox[0][1] = pointmat(0,j);
3835 if (bbox[1][0] > pointmat(1,j))
3837 bbox[1][0] = pointmat(1,j);
3839 if (bbox[1][1] < pointmat(1,j))
3841 bbox[1][1] = pointmat(1,j);
3843 if (bbox[2][0] > values(j))
3845 bbox[2][0] = values(j);
3847 if (bbox[2][1] < values(j))
3849 bbox[2][1] = values(j);
3854 mfem::out <<
"[xmin,xmax] = [" << bbox[0][0] <<
',' << bbox[0][1] <<
"]\n"
3855 <<
"[ymin,ymax] = [" << bbox[1][0] <<
',' << bbox[1][1] <<
"]\n"
3856 <<
"[zmin,zmax] = [" << bbox[2][0] <<
',' << bbox[2][1] <<
']'
3859 os <<
"endsolid GridFunction" << endl;
3876 MFEM_ASSERT(new_vertex.
Size() == mesh->
GetNV(),
"");
3880 for (
int i = 0; i < new_vertex.
Size(); i++)
3882 old_vertex[new_vertex[i]] = i;
3889 for (
int i = 0; i < mesh->
GetNV(); i++)
3894 for (
int j = 0; j < new_vdofs.
Size(); j++)
3896 tmp(new_vdofs[j]) = (*this)(old_vdofs[j]);
3902 for (
int i = 0; i < mesh->
GetNEdges(); i++)
3905 if (old_vertex[ev[0]] > old_vertex[ev[1]])
3910 for (
int k = 0; k < dofs.
Size(); k++)
3912 int new_dof = dofs[k];
3913 int old_dof = dofs[(ind[k] < 0) ? -1-ind[k] : ind[k]];
3920 real_t sign = (ind[k] < 0) ? -1.0 : 1.0;
3921 tmp(new_vdof) = sign * (*this)(old_vdof);
3934 int with_subdomains,
3942 int nfe = ufes->
GetNE();
3946 Vector ul, fl, fla, d_xyz;
3956 if (with_subdomains)
3961 real_t total_error = 0.0;
3962 for (
int s = 1;
s <= nsd;
s++)
3965 u.ComputeFlux(blfi, flux, with_coeff, (with_subdomains ?
s : -1));
3967 for (
int i = 0; i < nfe; i++)
3969 if (with_subdomains && ufes->
GetAttribute(i) !=
s) {
continue; }
3974 u.GetSubVector(udofs, ul);
3987 *ffes->
GetFE(i), fl, with_coeff);
3992 (aniso_flags ? &d_xyz : NULL));
3994 error_estimates(i) = std::sqrt(eng);
4000 for (
int k = 0; k <
dim; k++)
4007 for (
int k = 0; k <
dim; k++)
4009 if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
4012 (*aniso_flags)[i] = flag;
4020 auto process_local_error = total_error;
4021 MPI_Allreduce(&process_local_error, &total_error, 1,
4023 MPI_SUM, pfes->GetComm());
4026 return std::sqrt(total_error);
4038 MFEM_VERIFY(
dim >= 1,
"dim must be positive");
4039 MFEM_VERIFY(
dim <= 3,
"dim cannot be greater than 3");
4040 MFEM_VERIFY(order >= 0,
"order cannot be negative");
4042 bool rotate = (angle != 0.0) || (midpoint->
Norml2() != 0.0);
4051 x[0] = tmp[0]*cos(-angle) - tmp[1]*sin(-angle);
4052 x[1] = tmp[0]*sin(-angle) + tmp[1]*cos(-angle);
4061 real_t x1 = (x(0) - xmin(0))/(xmax(0)-xmin(0)), x2, x3;
4062 Vector poly_x(order+1), poly_y(order+1), poly_z(order+1);
4066 x2 = (x(1)-xmin(1))/(xmax(1)-xmin(1));
4071 x3 = (x(2)-xmin(2))/(xmax(2)-xmin(2));
4075 int basis_dimension =
static_cast<int>(pow(order+1,
dim));
4076 poly.
SetSize(basis_dimension);
4081 for (
int i = 0; i <= order; i++)
4083 poly(i) = poly_x(i);
4089 for (
int j = 0; j <= order; j++)
4091 for (
int i = 0; i <= order; i++)
4093 int cnt = i + (order+1) * j;
4094 poly(cnt) = poly_x(i) * poly_y(j);
4101 for (
int k = 0; k <= order; k++)
4103 for (
int j = 0; j <= order; j++)
4105 for (
int i = 0; i <= order; i++)
4107 int cnt = i + (order+1) * j + (order+1) * (order+1) * k;
4108 poly(cnt) = poly_x(i) * poly_y(j) * poly_z(k);
4116 MFEM_ABORT(
"TensorProductLegendre: invalid value of dim");
4132 int num_elems = patch.
Size();
4142 if (
rotate && iface >= 0)
4148 physical_diff = 0.0;
4151 for (
int i = 0; i < 2; i++)
4154 Tr.
Transform(reference_pt, physical_pt);
4155 midpoint += physical_pt;
4156 physical_pt *= pow(-1.0,i);
4157 physical_diff += physical_pt;
4160 angle = atan2(physical_diff(1),physical_diff(0));
4163 for (
int i = 0; i < num_elems; i++)
4165 int ielem = patch[i];
4176 transip -= midpoint;
4179 transip[0] = tmp[0]*cos(-angle) - tmp[1]*sin(-angle);
4180 transip[1] = tmp[0]*sin(-angle) + tmp[1]*cos(-angle);
4182 for (
int d = 0; d <
dim; d++) { xmax(d) = max(xmax(d), transip(d)); }
4183 for (
int d = 0; d <
dim; d++) { xmin(d) = min(xmin(d), transip(d)); }
4191 bool subdomain_reconstruction,
4195 MFEM_VERIFY(tichonov_coeff >= 0.0,
"tichonov_coeff cannot be negative");
4202 int nfe = ufes->
GetNE();
4203 int nfaces = ufes->
GetNF();
4210 error_estimates = 0.0;
4221 if (subdomain_reconstruction)
4226 real_t total_error = 0.0;
4227 for (
int iface = 0; iface < nfaces; iface++)
4234 patch[0] = el1; patch[1] = el2;
4237 if (el1 == -1 || el2 == -1)
4248 if (el1_attr != el2_attr) {
continue; }
4257 int num_basis_functions =
static_cast<int>(pow(patch_order+1,
dim));
4258 int flux_order = 2*patch_order + 1;
4267 xmin, xmax, angle, midpoint, iface);
4272 for (
int i = 0; i < patch.
Size(); i++)
4274 int ielem = patch[i];
4280 u.GetSubVector(udofs, ul);
4288 *dummy, fl, with_coeff, ir);
4292 for (
int k = 0; k < num_integration_pts; k++)
4304 for (
int l = 0; l < num_basis_functions; l++)
4307 for (
int n = 0; n < sdim; n++)
4309 b[l + n * num_basis_functions] +=
p(l) * fl(k + n * num_integration_pts);
4321 for (
int i = 0; i < num_basis_functions; i++)
4323 A(i,i) += tichonov_coeff;
4330 if (!lu.
Factor(num_basis_functions,TOL))
4333 mfem::out <<
"LSZZErrorEstimator: Matrix A is singular.\t"
4334 <<
"Consider increasing tichonov_coeff." << endl;
4335 for (
int i = 0; i < num_basis_functions; i++)
4339 lu.
Factor(num_basis_functions,TOL);
4341 lu.
Solve(num_basis_functions, sdim,
b);
4349 for (
int i = 0; i < num_basis_functions; i++)
4351 for (
int j = 0; j < sdim; j++)
4353 f(j) +=
b[i + j * num_basis_functions] *
p(i);
4360 real_t element_error = 0.0;
4361 real_t patch_error = 0.0;
4362 for (
int i = 0; i < patch.
Size(); i++)
4364 int ielem = patch[i];
4365 element_error =
u.ComputeElementGradError(ielem, &global_poly);
4366 element_error *= element_error;
4367 patch_error += element_error;
4368 error_estimates(ielem) += element_error;
4372 total_error += patch_error;
4380 for (
int ielem = 0; ielem < nfe; ielem++)
4382 if (counters[ielem] == 0)
4384 error_estimates(ielem) =
infinity();
4388 error_estimates(ielem) /= counters[ielem]/2.0;
4389 error_estimates(ielem) = sqrt(error_estimates(ielem));
4392 return std::sqrt(total_error/
dim);
4414 for (
int j = 0; j < nip; j++)
4426 errj = pow(errj,
p);
4431 norm = std::max(norm, errj);
4440 norm = -pow(-norm, 1./
p);
4444 norm = pow(norm, 1./
p);
4458 return sol_in.
Eval(*T_in, ip);
4469 string cname = name;
4470 if (cname ==
"Linear")
4474 else if (cname ==
"Quadratic")
4478 else if (cname ==
"Cubic")
4482 else if (!strncmp(name,
"H1_", 3))
4486 else if (!strncmp(name,
"H1Pos_", 6))
4491 else if (!strncmp(name,
"L2_T", 4))
4495 else if (!strncmp(name,
"L2_", 3))
4501 mfem::err <<
"Extrude1DGridFunction : unknown FE collection : "
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Piecewise-(bi)cubic continuous finite elements.
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
real_t Scale()
Return the scale factor times the optional time dependent function. Returns with when not set by th...
real_t Tol()
Return the tolerance used to identify the mesh vertices.
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void GetColumnReference(int c, Vector &col)
real_t * Data() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void Norm2(real_t *v) const
Take the 2-norm of the columns of A and store in v.
Class used for extruding scalar GridFunctions.
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i].
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name.
@ CONTINUOUS
Field is continuous across element interfaces.
virtual int GetContType() const =0
virtual const char * Name() const
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
void Save(std::ostream &out) const
Save finite element space to output stream out.
int GetNVDofs() const
Number of all scalar vertex dofs.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
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...
static void AdjustVDofs(Array< int > &vdofs)
Remove the orientation information encoded into an array of dofs Some basis function types have a rel...
bool Nonconforming() const
void GetVertexVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified vertices.
int GetNEDofs() const
Number of all scalar edge-interior dofs.
int GetAttribute(int i) const
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
int GetNBE() const
Returns number of boundary elements in the mesh.
const NURBSExtension * GetNURBSext() const
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
int GetBdrAttribute(int i) const
int GetLocalDofForDof(int i) const
Return the local dof index in the first element that contains dof i.
int GetElementForDof(int i) const
Return the index of the first element that contains dof i.
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller.
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Ordering::Type GetOrdering() const
Return the ordering method.
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
int GetNE() const
Returns number of elements in the mesh.
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
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...
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
int GetNFDofs() const
Number of all scalar face-interior dofs.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
const FiniteElementCollection * FEColl() const
Mesh * GetMesh() const
Returns the mesh.
int GetNV() const
Returns number of vertices in the mesh.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
int GetVDim() const
Returns vector dimension.
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
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...
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
Abstract class for all finite elements.
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
int GetDim() const
Returns the reference space dimension for the finite element.
virtual void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in reference space at the give...
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
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 CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in reference space at the gi...
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
const IntegrationRule * GetVertices(int GeomType) const
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType.
int NumBdr(int GeomType) const
Return the number of boundary "faces" of a given Geometry::Type.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps, int vdim=1) const
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, const Array< int > &bdr_attr, Array< int > &values_counter)
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
virtual void CountElementsPerVDof(Array< int > &elem_per_vdof) const
For each vdof, counts how many elements contain the vdof, as containment is determined by FiniteEleme...
virtual real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
void SaveVTK(std::ostream &out, const std::string &field_name, int ref)
Write the GridFunction in VTK format. Note that Mesh::PrintVTK must be called first....
virtual real_t ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, class JumpScaling jump_scaling, const IntegrationRule *irs[]=NULL) const
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
virtual real_t ComputeHCurlError(VectorCoefficient *exsol, VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(curl)-norm for ND elements.
virtual real_t ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, real_t Nu, int norm_type) const
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Extension of GetGradient(...) for a collection of IntegrationPoints.
FiniteElementCollection * fec
Used when the grid function is read from a file. It can also be set explicitly, see MakeOwner().
void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, const Array< int > &attr, Array< int > &values_counter)
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Compute the vector gradient with respect to the physical element variable.
virtual real_t ComputeMaxError(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
void ImposeBounds(int i, const Vector &weights, const Vector &lo_, const Vector &hi_)
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void MakeTRef(FiniteElementSpace *f, real_t *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
void AccumulateAndCountDerivativeValues(int comp, int der_comp, GridFunction &der, Array< int > &zones_per_dof)
Used for the serial and parallel implementations of the GetDerivative() method; see its documentation...
void GetVectorFieldValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr, int comp=0) const
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
void GetElementAverages(GridFunction &avgs) const
virtual real_t ComputeElementGradError(int ielem, VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 in element ielem for H1 or L2 elements.
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
void SaveSTL(std::ostream &out, int TimesToRefine=1)
Write the GridFunction in STL format. Note that the mesh dimension must be 2 and that quad elements w...
virtual void ComputeElementLpErrors(const real_t p, Coefficient &exsol, Vector &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
virtual void GetElementDofValues(int el, Vector &dof_vals) const
virtual real_t ComputeLpError(const real_t p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
FiniteElementSpace * FESpace()
void SaveSTLTri(std::ostream &out, real_t p1[], real_t p2[], real_t p3[])
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int GetFaceVectorValues(int i, int side, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
void GetValuesFrom(const GridFunction &orig_func)
void LegacyNCReorder()
Loading helper.
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, real_t &integral)
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
virtual real_t ComputeCurlError(VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns ||curl u_ex - curl u_h||_L2 for ND elements.
void GetBdrValuesFrom(const GridFunction &orig_func)
virtual real_t ComputeHDivError(VectorCoefficient *exsol, Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(div)-norm for RT elements.
virtual real_t ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, const Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
void GetDerivative(int comp, int der_comp, GridFunction &der)
Compute a certain derivative of a function's component. Derivatives of the function are computed at t...
int GetFaceValues(int i, int side, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, int vdim=1) const
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
GridFunction & operator=(const GridFunction &rhs)
Copy assignment. Only the data of the base class Vector is copied.
virtual real_t ComputeDivError(Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns ||div u_ex - div u_h||_L2 for RT elements.
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
virtual real_t ComputeGradError(VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 for H1 or L2 elements.
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
void RestrictConforming()
real_t GetDivergence(ElementTransformation &tr) const
void AccumulateAndCountZones(Coefficient &coeff, AvgType type, Array< int > &zones_per_vdof)
Accumulates (depending on type) the values of coeff at all shared vdofs and counts in how many zones ...
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
void GetCurl(ElementTransformation &tr, Vector &curl) const
void GetGradient(ElementTransformation &tr, Vector &grad) const
Gradient of a scalar function at a quadrature point.
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
Compute the vector gradient with respect to the reference element variable.
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
void GetHessians(int i, const IntegrationRule &ir, DenseMatrix &hess, int vdim=1) const
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh.
void GetVectorFieldNodalValues(Vector &val, int comp) const
void ProjectBdrCoefficient(Coefficient &coeff, const Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
void ReorderByNodes()
For a vector grid function, makes sure that the ordering is byNODES.
Arbitrary order H1-conforming (continuous) finite elements.
Class for integration point with weight.
void Set1w(const real_t x1, const real_t w)
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void SetRelTol(real_t rtol)
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
void SetAbsTol(real_t atol)
real_t Eval(real_t h, int p) const
Arbitrary order "L2-conforming" discontinuous finite elements.
virtual bool Factor(int m, real_t TOL=0.0)
Compute the LU factorization of the current matrix.
virtual void Solve(int m, int n, real_t *X) const
Piecewise-(bi/tri)linear continuous finite elements.
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
int GetNEdges() const
Return the number of edges.
void GetBdrElementFace(int i, int *f, int *o) const
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Geometry::Type GetElementGeometry(int i) const
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
bool Nonconforming() const
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.