15 #include "../mesh/nurbs.hpp"
16 #include "../general/text.hpp"
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);
111 const double *l_data = gf_array[i]->
GetData();
112 double *g_data =
data;
115 for (
int d = 0; d < vdim; d++)
117 memcpy(g_data+vi, l_data, l_nvdofs*
sizeof(
double));
120 memcpy(g_data+ei, l_data, l_nedofs*
sizeof(
double));
123 memcpy(g_data+fi, l_data, l_nfdofs*
sizeof(
double));
126 memcpy(g_data+di, l_data, l_nddofs*
sizeof(
double));
133 memcpy(g_data+vdim*vi, l_data, vdim*l_nvdofs*
sizeof(
double));
134 l_data += vdim*l_nvdofs;
135 g_data += vdim*g_nvdofs;
136 memcpy(g_data+vdim*ei, l_data, vdim*l_nedofs*
sizeof(
double));
137 l_data += vdim*l_nedofs;
138 g_data += vdim*g_nedofs;
139 memcpy(g_data+vdim*fi, l_data, vdim*l_nfdofs*
sizeof(
double));
140 l_data += vdim*l_nfdofs;
141 g_data += vdim*g_nfdofs;
142 memcpy(g_data+vdim*di, l_data, vdim*l_nddofs*
sizeof(
double));
143 l_data += vdim*l_nddofs;
144 g_data += vdim*g_nddofs;
182 old_data.
Swap(*
this);
185 T->
Mult(old_data, *
this);
262 int nfe = ufes->
GetNE();
270 for (
int i = 0; i < nfe; i++)
272 if (subdomain >= 0 && ufes->
GetAttribute(i) != subdomain)
284 *ffes->
GetFE(i), fl, wcoef);
289 for (
int j = 0; j < fdofs.
Size(); j++)
305 for (
int i = 0; i < count.Size(); i++)
307 if (count[i] != 0) { flux(i) /= count[i]; }
371 int dof = FElem->
GetDof();
381 "invalid FE map type");
383 for (k = 0; k < n; k++)
386 nval[k] = shape * ((
const double *)loc_data + dof * vdim);
393 for (k = 0; k < n; k++)
397 nval[k] = loc_data * (&vshape(0,vdim));
418 fe->CalcPhysShape(*Tr, DofVal);
422 return (DofVal * LocVec);
429 int dof = FElem->
GetDof();
449 for (
int k = 0; k < vdim; k++)
451 val(k) = shape * ((
const double *)loc_data + dof * k);
476 int dof = FElem->
GetDof();
477 Vector DofVal(dof), loc_data(dof);
481 for (
int k = 0; k < n; k++)
484 vals(k) = DofVal * loc_data;
490 for (
int k = 0; k < n; k++)
494 vals(k) = DofVal * loc_data;
523 "invalid FE map type");
525 int dof = FElem->
GetDof();
526 Vector DofLap(dof), loc_data(dof);
528 for (
int k = 0; k < n; k++)
533 laps(k) = DofLap * loc_data;
563 int size = (dim*(dim+1))/2;
566 "invalid FE map type");
568 int dof = FElem->
GetDof();
576 for (
int k = 0; k < n; k++)
582 for (
int i = 0; i <
size; i++)
584 for (
int d = 0; d < dof; d++)
586 hess(k,i) += DofHes(d,i) * loc_data[d];
671 fip.
x = 1.0 - ip.
x - ip.
y;
677 fip.
y = 1.0 - ip.
x - ip.
y;
716 int comp,
Vector *tr)
const
742 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \""
744 <<
"on mesh edges.");
757 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \""
759 <<
"on mesh faces.");
811 MFEM_ABORT(
"GridFunction::GetValue: Unsupported element type \""
829 return (DofVal * LocVec);
844 for (
int j = 0; j < nip; j++)
880 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \""
882 <<
"on mesh edges.");
895 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \""
897 <<
"on mesh faces.");
949 MFEM_ABORT(
"GridFunction::GetVectorValue: Unsupported element type \""
951 if (val.
Size() > 0) { val = NAN; }
972 for (
int k = 0; k < vdim; k++)
974 val(k) = shape * ((
const double *)loc_data + dof * k);
998 int dof = FElem->
GetDof();
1012 for (
int j = 0; j < nip; j++)
1018 for (
int k = 0; k < vdim; k++)
1020 vals(k,j) = shape * ((
const double *)loc_data + dof * k);
1032 for (
int j = 0; j < nip; j++)
1093 Vector shape, loc_values, orig_loc_values;
1094 int i, j, d, ne, dof, odof, vdim;
1098 for (i = 0; i < ne; i++)
1106 odof = orig_fe->
GetDof();
1107 loc_values.
SetSize(dof * vdim);
1110 for (j = 0; j < dof; j++)
1114 for (d = 0; d < vdim; d++)
1116 loc_values(d*dof+j) =
1117 shape * ((
const double *)orig_loc_values + d * odof) ;
1130 Vector shape, loc_values, orig_loc_values;
1131 int i, j, d, nbe, dof, odof, vdim;
1135 for (i = 0; i < nbe; i++)
1143 odof = orig_fe->
GetDof();
1144 loc_values.
SetSize(dof * vdim);
1147 for (j = 0; j < dof; j++)
1151 for (d = 0; d < vdim; d++)
1153 loc_values(d*dof+j) =
1154 shape * ((
const double *)orig_loc_values + d * odof);
1168 int d, j, k, n, sdim, dof, ind;
1175 int *dofs = &vdofs[comp*dof];
1181 for (k = 0; k < n; k++)
1186 for (d = 0; d < sdim; d++)
1189 for (j = 0; j < dof; j++)
1190 if ( (ind=dofs[j]) >= 0 )
1192 a += vshape(j, d) *
data[ind];
1196 a -= vshape(j, d) *
data[-1-ind];
1213 double *temp =
new double[
size];
1216 for (j = 0; j < ndofs; j++)
1217 for (i = 0; i < vdim; i++)
1219 temp[j+i*ndofs] =
data[k++];
1222 for (i = 0; i <
size; i++)
1250 val(vertices[k]) += vals(k, comp);
1251 overlap[vertices[k]]++;
1255 for (i = 0; i < overlap.Size(); i++)
1257 val(i) /= overlap[i];
1265 int d, i, k, ind, dof, sdim;
1274 for (i = 0; i < new_fes->
GetNE(); i++)
1281 for (d = 0; d < sdim; d++)
1283 for (k = 0; k < dof; k++)
1285 if ( (ind=new_vdofs[dof*d+k]) < 0 )
1287 ind = -1-ind, vals(k, d) = - vals(k, d);
1289 vec_field(ind) += vals(k, d);
1295 for (i = 0; i < overlap.Size(); i++)
1297 vec_field(i) /= overlap[i];
1308 Vector pt_grad, loc_func;
1309 int i, j, k,
dim, dof, der_dof, ind;
1312 for (i = 0; i < overlap.Size(); i++)
1319 for (i = 0; i < der_fes->
GetNE(); i++)
1328 der_dof = der_fe->
GetDof();
1334 for (j = 0; j < dof; j++)
1335 loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
1337 for (k = 0; k < der_dof; k++)
1345 for (j = 0; j <
dim; j++)
1347 a += inv_jac(j, der_comp) * pt_grad(j);
1349 der(der_dofs[k]) +=
a;
1350 overlap[der_dofs[k]]++;
1354 for (i = 0; i < overlap.Size(); i++)
1356 der(i) /= overlap[i];
1377 MultAtB(loc_data_mat, dshape, gh);
1391 "invalid FE map type");
1396 for (
int i = 0; i < Jinv.
Width(); i++)
1398 for (
int j = 0; j < Jinv.
Height(); j++)
1400 div_v += grad_hat(i, j) * Jinv(j, i);
1413 return (loc_data * divshape) / T.
Weight();
1459 MFEM_ABORT(
"GridFunction::GetDivergence: Unsupported element type \""
1477 "invalid FE map type");
1483 Mult(grad_hat, Jinv, grad);
1484 MFEM_ASSERT(grad.Height() == grad.Width(),
"");
1485 if (grad.Height() == 3)
1488 curl(0) = grad(2,1) - grad(1,2);
1489 curl(1) = grad(0,2) - grad(2,0);
1490 curl(2) = grad(1,0) - grad(0,1);
1492 else if (grad.Height() == 2)
1495 curl(0) = grad(1,0) - grad(0,1);
1507 curl.
SetSize(curl_shape.Width());
1508 if (curl_shape.Width() == 3)
1511 curl_shape.MultTranspose(loc_data, curl_hat);
1516 curl_shape.MultTranspose(loc_data, curl);
1562 MFEM_ABORT(
"GridFunction::GetCurl: Unsupported element type \""
1576 "invalid FE map type");
1630 MFEM_ABORT(
"GridFunction::GetGradient: Unsupported element type \""
1653 dshape.MultTranspose(lval, gh);
1674 Mult(grad_hat, Jinv, grad);
1717 MFEM_ABORT(
"GridFunction::GetVectorGradient: "
1718 "Unsupported element type \"" << T.
ElementType <<
"\"");
1728 Vector loc_avgs, loc_this;
1733 for (
int i = 0; i <
fes->
GetNE(); i++)
1738 avgs.FESpace()->GetElementDofs(i, te_dofs);
1741 loc_mass.
Mult(loc_this, loc_avgs);
1742 avgs.AddElementVector(te_dofs, loc_avgs);
1744 loc_mass.
Mult(loc_this, loc_avgs);
1745 int_psi.AddElementVector(te_dofs, loc_avgs);
1747 for (
int i = 0; i < avgs.Size(); i++)
1749 avgs(i) /= int_psi(i);
1766 if (!mesh->
GetNE()) {
return; }
1777 MFEM_VERIFY(vdim == src.
fes->
GetVDim(),
"incompatible vector dimensions!");
1782 for (
int i = 0; i < mesh->
GetNE(); i++)
1789 dest_lvec.SetSize(vdim*P.
Height());
1795 for (
int vd = 0; vd < vdim; vd++)
1810 Vector vals, new_vals(size);
1813 MFEM_ASSERT(weights.
Size() ==
size,
"Different # of weights and dofs.");
1814 MFEM_ASSERT(lo_.
Size() ==
size,
"Different # of lower bounds and dofs.");
1815 MFEM_ASSERT(hi_.
Size() ==
size,
"Different # of upper bounds and dofs.");
1818 double tol = 1.e-12;
1826 slbqp.
Mult(vals, new_vals);
1832 double min_,
double max_)
1837 Vector vals, new_vals(size);
1840 double max_val = vals.
Max();
1841 double min_val = vals.
Min();
1843 if (max_val <= min_)
1850 if (min_ <= min_val && max_val <= max_)
1855 Vector minv(size), maxv(size);
1856 minv = (min_ > min_val) ? min_ : min_val;
1857 maxv = (max_ < max_val) ? max_ : max_val;
1870 R->
Mult(*
this, tmp);
1871 P->
Mult(tmp, *
this);
1889 for (j = 0; j < vertices.
Size(); j++)
1891 nval(vertices[j]) += values[j];
1892 overlap[vertices[j]]++;
1895 for (i = 0; i < overlap.Size(); i++)
1897 nval(i) /= overlap[i];
1915 for (
int i = 0; i <
fes->
GetNE(); i++)
1923 for (
int j = 0; j < vdofs.
Size(); j++)
1927 MFEM_VERIFY(vals[j] != 0.0,
1928 "Coefficient has zeros, harmonic avg is undefined!");
1929 (*this)(vdofs[j]) += 1.0 / vals[j];
1933 (*this)(vdofs[j]) += vals[j];
1935 else { MFEM_ABORT(
"Not implemented"); }
1937 zones_per_vdof[vdofs[j]]++;
1956 for (
int i = 0; i <
fes->
GetNE(); i++)
1964 for (
int j = 0; j < vdofs.
Size(); j++)
1981 MFEM_VERIFY(vals[j] != 0.0,
1982 "Coefficient has zeros, harmonic avg is undefined!");
1983 (*this)(ldof) += isign / vals[j];
1987 (*this)(ldof) += isign*vals[j];
1990 else { MFEM_ABORT(
"Not implemented"); }
1992 zones_per_vdof[ldof]++;
2001 int i, j, fdof, d, ind, vdim;
2025 for (j = 0; j < fdof; j++)
2029 if (vcoeff) { vcoeff->
Eval(vc, *transf, ip); }
2030 for (d = 0; d < vdim; d++)
2032 if (!vcoeff && !coeff[d]) {
continue; }
2034 val = vcoeff ? vc(d) : coeff[d]->
Eval(*transf, ip);
2035 if ( (ind = vdofs[fdof*d+j]) < 0 )
2037 val = -val, ind = -1-ind;
2039 if (++values_counter[ind] == 1)
2045 (*this)(ind) += val;
2068 for (i = 0; i < bdr_edges.
Size(); i++)
2070 int edge = bdr_edges[i];
2072 if (vdofs.
Size() == 0) {
continue; }
2080 for (d = 0; d < vdim; d++)
2082 if (!coeff[d]) {
continue; }
2084 fe->
Project(*coeff[d], *transf, vals);
2085 for (
int k = 0; k < vals.
Size(); k++)
2087 ind = vdofs[d*vals.
Size()+k];
2088 if (++values_counter[ind] == 1)
2090 (*this)(ind) = vals(k);
2094 (*this)(ind) += vals(k);
2102 fe->
Project(*vcoeff, *transf, vals);
2103 for (
int k = 0; k < vals.
Size(); k++)
2106 if (++values_counter[ind] == 1)
2108 (*this)(ind) = vals(k);
2112 (*this)(ind) += vals(k);
2123 for (
int i = 0; i < dofs.
Size(); i++)
2126 double val = vals(i);
2127 if (k < 0) { k = -1 - k; val = -val; }
2128 if (++values_counter[k] == 1)
2163 fe->
Project(vcoeff, *T, lvec);
2164 accumulate_dofs(dofs, lvec, *
this, values_counter);
2174 for (
int i = 0; i < bdr_edges.
Size(); i++)
2176 int edge = bdr_edges[i];
2178 if (dofs.
Size() == 0) {
continue; }
2184 fe->
Project(vcoeff, *T, lvec);
2185 accumulate_dofs(dofs, lvec, *
this, values_counter);
2195 for (
int i = 0; i <
size; i++)
2197 const int nz = zones_per_vdof[i];
2198 if (nz) { (*this)(i) /= nz; }
2203 for (
int i = 0; i <
size; i++)
2205 const int nz = zones_per_vdof[i];
2206 if (nz) { (*this)(i) = nz/(*
this)(i); }
2211 MFEM_ABORT(
"invalid AvgType");
2226 const double *center = delta_coeff.
Center();
2227 const double *vert = mesh->
GetVertex(0);
2228 double min_dist, dist;
2232 min_dist =
Distance(center, vert, dim);
2233 for (
int i = 0; i < mesh->
GetNV(); i++)
2236 dist =
Distance(center, vert, dim);
2237 if (dist < min_dist)
2247 if (min_dist >= delta_coeff.
Tol())
2256 Vector vals, loc_mass_vals;
2257 for (
int i = 0; i < mesh->
GetNE(); i++)
2260 for (
int j = 0; j < vertices.
Size(); j++)
2261 if (vertices[j] == v_idx)
2271 loc_mass.Mult(vals, loc_mass_vals);
2272 integral += loc_mass_vals.
Sum();
2282 if (delta_c == NULL)
2287 for (
int i = 0; i <
fes->
GetNE(); i++)
2301 (*this) *= (delta_c->
Scale() / integral);
2314 for (
int i = 0; i < dofs.
Size(); i++)
2327 (*this)(vdof) = coeff.
Eval(*T, ip);
2357 for (
int i = 0; i < dofs.
Size(); i++)
2369 vcoeff.
Eval(val, *T, ip);
2373 (*this)(vdof) = val(vd);
2400 int i, j, fdof, d, ind, vdim;
2414 for (j = 0; j < fdof; j++)
2418 for (d = 0; d < vdim; d++)
2420 if (!coeff[d]) {
continue; }
2422 val = coeff[d]->
Eval(*transf, ip);
2423 if ( (ind = vdofs[fdof*d+j]) < 0 )
2425 val = -val, ind = -1-ind;
2445 for (
int i = 0; i <
fes->
GetNE(); i++)
2454 for (
int j = 0; j < vdofs.
Size(); j++)
2456 if (attr > dof_attr[vdofs[j]])
2458 (*this)(vdofs[j]) = vals[j];
2459 dof_attr[vdofs[j]] = attr;
2501 for (
int i = 0; i < values_counter.
Size(); i++)
2503 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2519 for (
int i = 0; i < values_counter.
Size(); i++)
2521 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2537 Vector vc(dim), nor(dim), lvec, shape;
2557 vcoeff.
Eval(vc, *T, ip);
2560 lvec.
Add(ip.
weight * (vc * nor), shape);
2572 Vector vc(dim), nor(dim), lvec;
2588 vcoeff.
Eval(vc, *T, ip);
2590 lvec(j) = (vc * nor);
2607 for (
int i = 0; i < values_counter.
Size(); i++)
2609 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2618 double error = 0.0,
a;
2623 int fdof, d, i, intorder, j, k;
2649 for (k = 0; k < fdof; k++)
2650 if (vdofs[fdof*d+k] >= 0)
2652 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2656 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2659 a -= exsol[d]->
Eval(*transf, ip);
2665 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2678 for (
int i = 0; i <
fes->
GetNE(); i++)
2680 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2682 int intorder = 2*fe->
GetOrder() + 3;
2694 exsol.
Eval(exact_vals, *T, *ir);
2697 vals.
Norm2(loc_errs);
2702 error += ip.
weight * T->
Weight() * (loc_errs(j) * loc_errs(j));
2706 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2721 for (
int i = 0; i <
fes->
GetNE(); i++)
2741 exgrad->
Eval(vec,*Tr,ip);
2746 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2759 int n = (dim == 3) ? dim : 1;
2762 for (
int i = 0; i <
fes->
GetNE(); i++)
2782 excurl->
Eval(vec,*Tr,ip);
2788 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2794 double error = 0.0,
a;
2800 for (
int i = 0; i <
fes->
GetNE(); i++)
2824 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2832 int fdof, intorder, k;
2837 Vector shape, el_dofs, err_val, ell_coeff_val;
2859 intorder = 2 * intorder;
2873 transf = face_elem_transf->
Elem1;
2879 for (k = 0; k < fdof; k++)
2882 el_dofs(k) = (*this)(vdofs[k]);
2886 el_dofs(k) = - (*this)(-1-vdofs[k]);
2893 ell_coeff_val(j) = ell_coeff->
Eval(*transf, eip);
2894 err_val(j) = exsol->
Eval(*transf, eip) - (shape * el_dofs);
2900 transf = face_elem_transf->
Elem2;
2906 for (k = 0; k < fdof; k++)
2909 el_dofs(k) = (*this)(vdofs[k]);
2913 el_dofs(k) = - (*this)(-1-vdofs[k]);
2920 ell_coeff_val(j) += ell_coeff->
Eval(*transf, eip);
2921 ell_coeff_val(j) *= 0.5;
2922 err_val(j) -= (exsol->
Eval(*transf, eip) - (shape * el_dofs));
2926 transf = face_elem_transf;
2931 double nu = jump_scaling.
Eval(h, p);
2932 error += (ip.
weight * nu * ell_coeff_val(j) *
2934 err_val(j) * err_val(j));
2938 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2953 int norm_type)
const
2955 double error1 = 0.0;
2956 double error2 = 0.0;
2964 return sqrt(error1 * error1 + error2 * error2);
2973 return sqrt(L2error*L2error + GradError*GradError);
2982 return sqrt(L2error*L2error + DivError*DivError);
2991 return sqrt(L2error*L2error + CurlError*CurlError);
2997 double error = 0.0,
a;
3002 int fdof, d, i, intorder, j, k;
3029 for (k = 0; k < fdof; k++)
3030 if (vdofs[fdof*d+k] >= 0)
3032 a += (*this)(vdofs[fdof*d+k]) * shape(k);
3036 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
3038 a -= exsol[d]->
Eval(*transf, ip);
3055 int i, fdof,
dim, intorder, j, k;
3059 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
3062 double a, error = 0.0;
3071 for (i = 0; i < mesh->
GetNE(); i++)
3073 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3090 for (k = 0; k < fdof; k++)
3093 el_dofs(k) = (*this)(vdofs[k]);
3097 el_dofs(k) = -(*this)(-1-vdofs[k]);
3104 a = (el_dofs * shape) - (exsol->
Eval(*transf, ip));
3110 for (i = 0; i < mesh->
GetNE(); i++)
3112 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3130 for (k = 0; k < fdof; k++)
3133 el_dofs(k) = (*this)(vdofs[k]);
3137 el_dofs(k) = -(*this)(-1-vdofs[k]);
3144 exgrad->
Eval(e_grad, *transf, ip);
3146 Mult(dshape, Jinv, dshapet);
3165 for (
int i = 0; i <
fes->
GetNE(); i++)
3175 int intorder = 2*fe->
GetOrder() + 3;
3184 double err = fabs(vals(j) - exsol.
Eval(*T, ip));
3190 err *= weight->
Eval(*T, ip);
3198 err *= weight->
Eval(*T, ip);
3200 error = std::max(error, err);
3210 error = -pow(-error, 1./p);
3214 error = pow(error, 1./p);
3227 "Incorrect size for result vector");
3234 for (
int i = 0; i <
fes->
GetNE(); i++)
3244 int intorder = 2*fe->
GetOrder() + 3;
3253 double err = fabs(vals(j) - exsol.
Eval(*T, ip));
3259 err *= weight->
Eval(*T, ip);
3267 err *= weight->
Eval(*T, ip);
3269 error[i] = std::max(error[i], err);
3277 error[i] = -pow(-error[i], 1./p);
3281 error[i] = pow(error[i], 1./p);
3298 for (
int i = 0; i <
fes->
GetNE(); i++)
3308 int intorder = 2*fe->
GetOrder() + 3;
3313 exsol.
Eval(exact_vals, *T, *ir);
3320 vals.
Norm2(loc_errs);
3324 v_weight->
Eval(exact_vals, *T, *ir);
3327 for (
int j = 0; j < vals.
Width(); j++)
3330 for (
int d = 0; d < vals.
Height(); d++)
3332 err += vals(d,j)*exact_vals(d,j);
3334 loc_errs(j) = fabs(err);
3341 double err = loc_errs(j);
3347 err *= weight->
Eval(*T, ip);
3355 err *= weight->
Eval(*T, ip);
3357 error = std::max(error, err);
3367 error = -pow(-error, 1./p);
3371 error = pow(error, 1./p);
3386 "Incorrect size for result vector");
3394 for (
int i = 0; i <
fes->
GetNE(); i++)
3404 int intorder = 2*fe->
GetOrder() + 3;
3409 exsol.
Eval(exact_vals, *T, *ir);
3416 vals.
Norm2(loc_errs);
3420 v_weight->
Eval(exact_vals, *T, *ir);
3423 for (
int j = 0; j < vals.
Width(); j++)
3426 for (
int d = 0; d < vals.
Height(); d++)
3428 err += vals(d,j)*exact_vals(d,j);
3430 loc_errs(j) = fabs(err);
3437 double err = loc_errs(j);
3443 err *= weight->
Eval(*T, ip);
3451 err *= weight->
Eval(*T, ip);
3453 error[i] = std::max(error[i], err);
3461 error[i] = -pow(-error[i], 1./p);
3465 error[i] = pow(error[i], 1./p);
3492 out <<
"NURBS_patches\n";
3511 ofstream ofs(fname);
3512 ofs.precision(precision);
3516 #ifdef MFEM_USE_ADIOS2
3518 const std::string& variable_name,
3521 out.
Save(*
this, variable_name, type);
3537 out <<
"SCALARS " << field_name <<
" double 1\n"
3538 <<
"LOOKUP_TABLE default\n";
3539 for (
int i = 0; i < mesh->
GetNE(); i++)
3546 for (
int j = 0; j < val.
Size(); j++)
3548 out << val(j) <<
'\n';
3552 else if ( (vec_dim == 2 || vec_dim == 3) && mesh->
SpaceDimension() > 1)
3555 out <<
"VECTORS " << field_name <<
" double\n";
3556 for (
int i = 0; i < mesh->
GetNE(); i++)
3565 for (
int j = 0; j < vval.
Width(); j++)
3567 out << vval(0, j) <<
' ' << vval(1, j) <<
' ';
3583 for (
int vd = 0; vd < vec_dim; vd++)
3585 out <<
"SCALARS " << field_name << vd <<
" double 1\n"
3586 <<
"LOOKUP_TABLE default\n";
3587 for (
int i = 0; i < mesh->
GetNE(); i++)
3594 for (
int j = 0; j < val.
Size(); j++)
3596 out << val(j) <<
'\n';
3607 double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
3608 double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
3609 double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
3610 v1[2] * v2[0] - v1[0] * v2[2],
3611 v1[0] * v2[1] - v1[1] * v2[0]
3613 double rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
3614 n[0] *= rl; n[1] *= rl; n[2] *= rl;
3616 out <<
" facet normal " << n[0] <<
' ' << n[1] <<
' ' << n[2]
3618 <<
"\n vertex " << p1[0] <<
' ' << p1[1] <<
' ' << p1[2]
3619 <<
"\n vertex " << p2[0] <<
' ' << p2[1] <<
' ' << p2[2]
3620 <<
"\n vertex " << p3[0] <<
' ' << p3[1] <<
' ' << p3[2]
3621 <<
"\n endloop\n endfacet\n";
3637 double pts[4][3], bbox[3][2];
3639 out <<
"solid GridFunction\n";
3641 bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
3642 bbox[2][0] = bbox[2][1] = 0.0;
3643 for (i = 0; i < mesh->
GetNE(); i++)
3650 for (k = 0; k < RG.
Size()/n; k++)
3652 for (j = 0; j < n; j++)
3655 pts[j][0] = pointmat(0,l);
3656 pts[j][1] = pointmat(1,l);
3657 pts[j][2] = values(l);
3673 bbox[0][0] = pointmat(0,0);
3674 bbox[0][1] = pointmat(0,0);
3675 bbox[1][0] = pointmat(1,0);
3676 bbox[1][1] = pointmat(1,0);
3677 bbox[2][0] = values(0);
3678 bbox[2][1] = values(0);
3681 for (j = 0; j < values.
Size(); j++)
3683 if (bbox[0][0] > pointmat(0,j))
3685 bbox[0][0] = pointmat(0,j);
3687 if (bbox[0][1] < pointmat(0,j))
3689 bbox[0][1] = pointmat(0,j);
3691 if (bbox[1][0] > pointmat(1,j))
3693 bbox[1][0] = pointmat(1,j);
3695 if (bbox[1][1] < pointmat(1,j))
3697 bbox[1][1] = pointmat(1,j);
3699 if (bbox[2][0] > values(j))
3701 bbox[2][0] = values(j);
3703 if (bbox[2][1] < values(j))
3705 bbox[2][1] = values(j);
3710 mfem::out <<
"[xmin,xmax] = [" << bbox[0][0] <<
',' << bbox[0][1] <<
"]\n"
3711 <<
"[ymin,ymax] = [" << bbox[1][0] <<
',' << bbox[1][1] <<
"]\n"
3712 <<
"[zmin,zmax] = [" << bbox[2][0] <<
',' << bbox[2][1] <<
']'
3715 out <<
"endsolid GridFunction" << endl;
3732 MFEM_ASSERT(new_vertex.Size() == mesh->
GetNV(),
"");
3735 old_vertex.SetSize(new_vertex.Size());
3736 for (
int i = 0; i < new_vertex.Size(); i++)
3738 old_vertex[new_vertex[i]] = i;
3745 for (
int i = 0; i < mesh->
GetNV(); i++)
3750 for (
int j = 0; j < new_vdofs.
Size(); j++)
3752 tmp(new_vdofs[j]) = (*this)(old_vdofs[j]);
3758 for (
int i = 0; i < mesh->
GetNEdges(); i++)
3761 if (old_vertex[ev[0]] > old_vertex[ev[1]])
3766 for (
int k = 0; k < dofs.
Size(); k++)
3768 int new_dof = dofs[k];
3769 int old_dof = dofs[(ind[k] < 0) ? -1-ind[k] : ind[k]];
3776 double sign = (ind[k] < 0) ? -1.0 : 1.0;
3777 tmp(new_vdof) = sign * (*this)(old_vdof);
3789 const char *msg =
"invalid input stream";
3795 in >> ident; MFEM_VERIFY(ident ==
"VDim:", msg);
3822 out <<
"VDim: " <<
vdim <<
'\n'
3839 int with_subdomains,
3847 int nfe = ufes->
GetNE();
3851 Vector ul, fl, fla, d_xyz;
3861 if (with_subdomains)
3866 double total_error = 0.0;
3867 for (
int s = 1;
s <= nsd;
s++)
3870 u.
ComputeFlux(blfi, flux, with_coeff, (with_subdomains ?
s : -1));
3872 for (
int i = 0; i < nfe; i++)
3874 if (with_subdomains && ufes->
GetAttribute(i) !=
s) {
continue; }
3884 *ffes->
GetFE(i), fl, with_coeff);
3889 (aniso_flags ? &d_xyz : NULL));
3891 error_estimates(i) = std::sqrt(err);
3897 for (
int k = 0; k <
dim; k++)
3902 double thresh = 0.15 * 3.0/
dim;
3904 for (
int k = 0; k <
dim; k++)
3906 if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
3909 (*aniso_flags)[i] = flag;
3917 auto process_local_error = total_error;
3918 MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
3919 MPI_SUM, pfes->GetComm());
3921 #endif // MFEM_USE_MPI
3922 return std::sqrt(total_error);
3945 for (
int j = 0; j < nip; j++)
3962 norm = std::max(norm, err);
3971 norm = -pow(-norm, 1./p);
3975 norm = pow(norm, 1./p);
3989 return sol_in.
Eval(*T_in, ip);
4000 string cname = name;
4001 if (cname ==
"Linear")
4005 else if (cname ==
"Quadratic")
4009 else if (cname ==
"Cubic")
4013 else if (!strncmp(name,
"H1_", 3))
4017 else if (!strncmp(name,
"H1Pos_", 6))
4022 else if (!strncmp(name,
"L2_T", 4))
4026 else if (!strncmp(name,
"L2_", 3))
4032 mfem::err <<
"Extrude1DGridFunction : unknown FE collection : "
int GetNPoints() const
Returns the number of the points in the integration rule.
void SaveSTLTri(std::ostream &out, double p1[], double p2[], double p3[])
Abstract class for all finite elements.
void LegacyToNewVertexOrdering(Array< int > &order) const
I/O: Return a map from old (v1.1) vertex indices to new vertex indices.
Ordering::Type GetOrdering() const
Return the ordering method.
int Size() const
Return the logical size of the array.
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
int GetDim() const
Returns the reference space dimension for the finite element.
int GetAttribute(int i) const
int GetNDofs() const
Returns number of degrees of freedom.
void Save(const GridFunction &grid_function, const std::string &variable_name, const data_type type)
Class for an integration rule - an Array of IntegrationPoint.
void GetElementAverages(GridFunction &avgs) const
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Class for grid function - Vector with associated FE space.
int DofToVDof(int dof, int vd, int ndofs=-1) const
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
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...
Class used for extruding scalar GridFunctions.
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
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.
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Base class for vector Coefficients that optionally depend on time and space.
virtual double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
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 Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
virtual void ComputeElementLpErrors(const double p, Coefficient &exsol, Vector &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
virtual double ComputeCurlError(VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns ||curl u_ex - curl u_h||_L2 for ND elements.
virtual int GetContType() const =0
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
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.
void RestrictConforming()
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
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...
GridFunction * Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d, GridFunction *sol, const int ny)
Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
double Scale()
Return the scale factor times the optional time dependent function. Returns with when not set by th...
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
bool own_qspace
QuadratureSpace ownership flag.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int GetSize() const
Return the total number of quadrature points.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
double Norml2() const
Returns the l2 norm of the vector.
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
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...
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
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...
Piecewise-(bi/tri)linear continuous finite elements.
Data type dense matrix using column-major storage.
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
int Size() const
Returns the size of the vector.
int GetBdrAttribute(int i) const
int GetNE() const
Returns number of elements.
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...
double Tol()
Return the tolerance used to identify the mesh vertices.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, Array< int > &bdr_attr, Array< int > &values_counter)
Abstract parallel finite element space.
virtual double ComputeHDivError(VectorCoefficient *exsol, Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(div)-norm for RT elements.
int vdim
Vector dimension.
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
int GetNV() const
Returns number of vertices in the mesh.
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
virtual void GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
double * GetData() const
Return a pointer to the beginning of the Vector data.
void CalcOrtho(const DenseMatrix &J, Vector &n)
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
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...
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Geometry::Type GetElementBaseGeometry(int i) const
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr)
void GetVectorFieldNodalValues(Vector &val, int comp) const
void ImposeBounds(int i, const Vector &weights, const Vector &lo_, const Vector &hi_)
Vector & operator=(const double *v)
Copy Size() entries from v.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
virtual void GetElementDofValues(int el, Vector &dof_vals) const
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Piecewise-(bi)cubic continuous finite elements.
int GetNE() const
Returns number of elements in the mesh.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void GetDerivative(int comp, int der_comp, GridFunction &der)
bool Nonconforming() const
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
void SetPrintLevel(int print_lvl)
double f(const Vector &xvec)
void GetHessians(int i, const IntegrationRule &ir, DenseMatrix &hess, int vdim=1) const
double GetElementSize(ElementTransformation *T, int type=0)
int GetNBE() const
Returns number of boundary elements in the mesh.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void Swap(Vector &other)
Swap the contents of two Vectors.
Mesh * GetMesh() const
Returns the mesh.
void GetVectorFieldValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr, int comp=0) const
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
void SetMaxIter(int max_it)
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...
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetLinearConstraint(const Vector &w_, double a_)
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
int GetElementForDof(int i) const
Return the index of the first element that contains dof i.
QuadratureFunction()
Create an empty QuadratureFunction.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
GeometryRefiner GlobGeometryRefiner
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
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 ...
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps, int vdim=1) const
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
int GetVDim()
Returns dimension of the vector.
virtual void GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'.
int GetVDim() const
Returns vector dimension.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
void LoadSolution(std::istream &input, GridFunction &sol) const
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
FiniteElementSpace * FESpace()
GridFunction & operator=(const GridFunction &rhs)
Copy assignment. Only the data of the base class Vector is copied.
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
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...
virtual double ComputeLpError(const double p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
QuadratureFunction & operator=(double value)
Redefine '=' for QuadratureFunction = constant.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
int SpaceDimension() const
double Min() const
Returns the minimal element of the vector.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
double Norml1() const
Returns the l_1 norm of the vector.
virtual double ComputeGradError(VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 for H1 or L2 elements.
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
bool IsIdentityProlongation(const Operator *P)
bool IsLegacyLoaded() const
I/O: Return true if the mesh was loaded from the legacy v1.1 format.
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
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]...
void SetAbsTol(double atol)
double p(const Vector &x, double t)
void SetRelTol(double rtol)
const NURBSExtension * GetNURBSext() const
double Distance(const double *x, const double *y, const int n)
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
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.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
int GetLocalDofForDof(int i) const
Return the local dof index in the first element that contains dof i.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
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 ...
QuadratureSpace * qspace
Associated QuadratureSpace.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
virtual double ComputeMaxError(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
int NumBdr(int GeomType)
Return the number of boundary "faces" of a given Geometry::Type.
void GetColumnReference(int c, Vector &col)
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...
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Field is continuous across element interfaces.
FiniteElementCollection * fec
Used when the grid function is read from a file. It can also be set explicitly, see MakeOwner()...
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
double GetDivergence(ElementTransformation &tr) const
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
double ComputeElementLpDistance(double p, int i, GridFunction &gf1, GridFunction &gf2)
Compute the Lp distance between two grid functions on the given element.
virtual double ComputeDivError(Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns ||div u_ex - div u_h||_L2 for RT elements.
int GetFaceVectorValues(int i, int side, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
virtual double ComputeHCurlError(VectorCoefficient *exsol, VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(curl)-norm for ND elements.
bool Nonconforming() const
double ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains, bool with_coeff)
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 GetNVDofs() const
Number of all scalar vertex dofs.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
virtual const char * Name() const
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Class for integration point with weight.
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...
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
double Max() const
Returns the maximal element of the vector.
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
void NewMemoryAndSize(const Memory< double > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
int GetFaceValues(int i, int side, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, int vdim=1) const
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
void LegacyNCReorder()
Loading helper.
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
virtual void Mult(const Vector &xt, Vector &x) const
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Piecewise-(bi)quadratic continuous finite elements.
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
void ReorderByNodes()
For a vector grid function, makes sure that the ordering is byNODES.
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
NCMesh * ncmesh
Optional non-conforming mesh extension.
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th face element (2D and 3D).
virtual double ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th edge.
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.
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
void GetCurl(ElementTransformation &tr, Vector &curl) const
const FiniteElementCollection * FEColl() const
void Mult(const double *x, double *y) const
Matrix vector multiplication.
void be_to_bfe(Geometry::Type geom, int o, const IntegrationPoint &ip, IntegrationPoint &fip)
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Arbitrary order H1-conforming (continuous) finite elements.
void pts(int iphi, int t, double x[])
virtual double ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, class JumpScaling jump_scaling, const IntegrationRule *irs[]=NULL) const
Class representing the storage layout of a QuadratureFunction.
double u(const Vector &xvec)
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
void SetBounds(const Vector &lo_, const Vector &hi_)
double Eval(double h, int p) const
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
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...
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
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.
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
double Sum() const
Return the sum of the vector entries.
Class representing a function through its values (scalar or vector) at quadrature points...
void GetValuesFrom(const GridFunction &orig_func)
int GetNFDofs() const
Number of all scalar face-interior dofs.
int GetNEDofs() const
Number of all scalar edge-interior dofs.
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
void GetBdrValuesFrom(const GridFunction &orig_func)
static void AdjustVDofs(Array< int > &vdofs)
Array< int > attributes
A list of all unique element attributes used by the Mesh.
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr, Array< int > &values_counter)
void Save(std::ostream &out) const
Write the QuadratureFunction to the stream out.
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. The parameter ref > 0 must match the one used in Mesh::PrintVTK.
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
void GetGradient(ElementTransformation &tr, Vector &grad) const
Arbitrary order "L2-conforming" discontinuous finite elements.
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const