15 #include "../mesh/nurbs.hpp"
16 #include "../general/text.hpp"
40 istream::int_type next_char = input.peek();
46 if (buff ==
"NURBS_patches")
49 "NURBS_patches requires NURBS FE space");
54 MFEM_ABORT(
"unknown section: " << buff);
88 int g_nddofs = g_ndofs - (g_nvdofs + g_nedofs + g_nfdofs);
90 vi = ei = fi = di = 0;
91 for (
int i = 0; i < num_pieces; i++)
98 int l_nddofs = l_ndofs - (l_nvdofs + l_nedofs + l_nfdofs);
99 const double *l_data = gf_array[i]->
GetData();
100 double *g_data =
data;
103 for (
int d = 0; d < vdim; d++)
105 memcpy(g_data+vi, l_data, l_nvdofs*
sizeof(
double));
108 memcpy(g_data+ei, l_data, l_nedofs*
sizeof(
double));
111 memcpy(g_data+fi, l_data, l_nfdofs*
sizeof(
double));
114 memcpy(g_data+di, l_data, l_nddofs*
sizeof(
double));
121 memcpy(g_data+vdim*vi, l_data, vdim*l_nvdofs*
sizeof(
double));
122 l_data += vdim*l_nvdofs;
123 g_data += vdim*g_nvdofs;
124 memcpy(g_data+vdim*ei, l_data, vdim*l_nedofs*
sizeof(
double));
125 l_data += vdim*l_nedofs;
126 g_data += vdim*g_nedofs;
127 memcpy(g_data+vdim*fi, l_data, vdim*l_nfdofs*
sizeof(
double));
128 l_data += vdim*l_nfdofs;
129 g_data += vdim*g_nfdofs;
130 memcpy(g_data+vdim*di, l_data, vdim*l_nddofs*
sizeof(
double));
131 l_data += vdim*l_nddofs;
132 g_data += vdim*g_nddofs;
160 MFEM_ABORT(
"Error in update sequence. GridFunction needs to be updated "
161 "right after the space is updated.");
169 old_data.
Swap(*
this);
172 T->
Mult(old_data, *
this);
251 int nfe = ufes->
GetNE();
259 for (
int i = 0; i < nfe; i++)
261 if (subdomain >= 0 && ufes->
GetAttribute(i) != subdomain)
273 *ffes->
GetFE(i), fl, wcoef);
278 for (
int j = 0; j < fdofs.
Size(); j++)
294 for (
int i = 0; i < count.Size(); i++)
296 if (count[i] != 0) { flux(i) /= count[i]; }
327 tv.
MakeRef(const_cast<GridFunction &>(*
this), 0,
size);
360 int dof = FElem->
GetDof();
370 "invalid FE map type");
372 for (k = 0; k < n; k++)
375 nval[k] = shape * ((
const double *)loc_data + dof * vdim);
382 for (k = 0; k < n; k++)
386 nval[k] = loc_data * (&vshape(0,vdim));
407 fe->CalcPhysShape(*Tr, DofVal);
411 return (DofVal * LocVec);
418 int dof = FElem->
GetDof();
438 for (
int k = 0; k < vdim; k++)
440 val(k) = shape * ((
const double *)loc_data + dof * k);
466 "invalid FE map type");
467 int dof = FElem->
GetDof();
468 Vector DofVal(dof), loc_data(dof);
470 for (
int k = 0; k < n; k++)
473 vals(k) = DofVal * loc_data;
501 "invalid FE map type");
503 int dof = FElem->
GetDof();
504 Vector DofLap(dof), loc_data(dof);
506 for (
int k = 0; k < n; k++)
511 laps(k) = DofLap * loc_data;
541 int size = (dim*(dim+1))/2;
544 "invalid FE map type");
546 int dof = FElem->
GetDof();
554 for (
int k = 0; k < n; k++)
560 for (
int i = 0; i <
size; i++)
562 for (
int d = 0; d < dof; d++)
564 hess(k,i) += DofHes(d,i) * loc_data[d];
649 fip.
x = 1.0 - ip.
x - ip.
y;
655 fip.
y = 1.0 - ip.
x - ip.
y;
694 int comp,
Vector *tr)
const
720 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \""
722 <<
"on mesh edges.");
735 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \""
737 <<
"on mesh faces.");
789 MFEM_ABORT(
"GridFunction::GetValue: Unsupported element type \""
807 return (DofVal * LocVec);
822 for (
int j = 0; j < nip; j++)
858 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \""
860 <<
"on mesh edges.");
873 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \""
875 <<
"on mesh faces.");
927 MFEM_ABORT(
"GridFunction::GetVectorValue: Unsupported element type \""
929 if (val.
Size() > 0) { val = NAN; }
950 for (
int k = 0; k < vdim; k++)
952 val(k) = shape * ((
const double *)loc_data + dof * k);
976 int dof = FElem->
GetDof();
988 "invalid FE map type");
992 for (
int j = 0; j < nip; j++)
997 for (
int k = 0; k < vdim; k++)
999 vals(k,j) = shape * ((
const double *)loc_data + dof * k);
1011 for (
int j = 0; j < nip; j++)
1072 Vector shape, loc_values, orig_loc_values;
1073 int i, j, d, ne, dof, odof, vdim;
1077 for (i = 0; i < ne; i++)
1085 odof = orig_fe->
GetDof();
1086 loc_values.
SetSize(dof * vdim);
1089 for (j = 0; j < dof; j++)
1093 for (d = 0; d < vdim; d++)
1095 loc_values(d*dof+j) =
1096 shape * ((
const double *)orig_loc_values + d * odof) ;
1109 Vector shape, loc_values, orig_loc_values;
1110 int i, j, d, nbe, dof, odof, vdim;
1114 for (i = 0; i < nbe; i++)
1122 odof = orig_fe->
GetDof();
1123 loc_values.
SetSize(dof * vdim);
1126 for (j = 0; j < dof; j++)
1130 for (d = 0; d < vdim; d++)
1132 loc_values(d*dof+j) =
1133 shape * ((
const double *)orig_loc_values + d * odof);
1147 int d, j, k, n, sdim, dof, ind;
1154 int *dofs = &vdofs[comp*dof];
1160 for (k = 0; k < n; k++)
1165 for (d = 0; d < sdim; d++)
1168 for (j = 0; j < dof; j++)
1169 if ( (ind=dofs[j]) >= 0 )
1171 a += vshape(j, d) *
data[ind];
1175 a -= vshape(j, d) *
data[-1-ind];
1192 double *temp =
new double[
size];
1195 for (j = 0; j < ndofs; j++)
1196 for (i = 0; i < vdim; i++)
1198 temp[j+i*ndofs] =
data[k++];
1201 for (i = 0; i <
size; i++)
1229 val(vertices[k]) += vals(k, comp);
1230 overlap[vertices[k]]++;
1234 for (i = 0; i < overlap.Size(); i++)
1236 val(i) /= overlap[i];
1244 int d, i, k, ind, dof, sdim;
1253 for (i = 0; i < new_fes->
GetNE(); i++)
1260 for (d = 0; d < sdim; d++)
1262 for (k = 0; k < dof; k++)
1264 if ( (ind=new_vdofs[dof*d+k]) < 0 )
1266 ind = -1-ind, vals(k, d) = - vals(k, d);
1268 vec_field(ind) += vals(k, d);
1274 for (i = 0; i < overlap.Size(); i++)
1276 vec_field(i) /= overlap[i];
1287 Vector pt_grad, loc_func;
1288 int i, j, k,
dim, dof, der_dof, ind;
1291 for (i = 0; i < overlap.Size(); i++)
1298 for (i = 0; i < der_fes->
GetNE(); i++)
1307 der_dof = der_fe->
GetDof();
1313 for (j = 0; j < dof; j++)
1314 loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
1316 for (k = 0; k < der_dof; k++)
1324 for (j = 0; j <
dim; j++)
1326 a += inv_jac(j, der_comp) * pt_grad(j);
1328 der(der_dofs[k]) +=
a;
1329 overlap[der_dofs[k]]++;
1333 for (i = 0; i < overlap.Size(); i++)
1335 der(i) /= overlap[i];
1356 MultAtB(loc_data_mat, dshape, gh);
1370 "invalid FE map type");
1375 for (
int i = 0; i < Jinv.
Width(); i++)
1377 for (
int j = 0; j < Jinv.
Height(); j++)
1379 div_v += grad_hat(i, j) * Jinv(j, i);
1392 return (loc_data * divshape) / T.
Weight();
1438 MFEM_ABORT(
"GridFunction::GetDivergence: Unsupported element type \""
1456 "invalid FE map type");
1462 Mult(grad_hat, Jinv, grad);
1463 MFEM_ASSERT(grad.Height() == grad.Width(),
"");
1464 if (grad.Height() == 3)
1467 curl(0) = grad(2,1) - grad(1,2);
1468 curl(1) = grad(0,2) - grad(2,0);
1469 curl(2) = grad(1,0) - grad(0,1);
1471 else if (grad.Height() == 2)
1474 curl(0) = grad(1,0) - grad(0,1);
1486 curl.
SetSize(curl_shape.Width());
1487 if (curl_shape.Width() == 3)
1490 curl_shape.MultTranspose(loc_data, curl_hat);
1495 curl_shape.MultTranspose(loc_data, curl);
1541 MFEM_ABORT(
"GridFunction::GetCurl: Unsupported element type \""
1555 "invalid FE map type");
1611 MFEM_ABORT(
"GridFunction::GetGradient: Unsupported element type \""
1634 dshape.MultTranspose(lval, gh);
1655 Mult(grad_hat, Jinv, grad);
1698 MFEM_ABORT(
"GridFunction::GetVectorGradient: "
1699 "Unsupported element type \"" << T.
ElementType <<
"\"");
1709 Vector loc_avgs, loc_this;
1714 for (
int i = 0; i <
fes->
GetNE(); i++)
1719 avgs.FESpace()->GetElementDofs(i, te_dofs);
1722 loc_mass.
Mult(loc_this, loc_avgs);
1723 avgs.AddElementVector(te_dofs, loc_avgs);
1725 loc_mass.
Mult(loc_this, loc_avgs);
1726 int_psi.AddElementVector(te_dofs, loc_avgs);
1728 for (
int i = 0; i < avgs.Size(); i++)
1730 avgs(i) /= int_psi(i);
1740 if (!mesh->
GetNE()) {
return; }
1751 MFEM_VERIFY(vdim == src.
fes->
GetVDim(),
"incompatible vector dimensions!");
1756 for (
int i = 0; i < mesh->
GetNE(); i++)
1763 dest_lvec.SetSize(vdim*P.
Height());
1769 for (
int vd = 0; vd < vdim; vd++)
1784 Vector vals, new_vals(size);
1787 MFEM_ASSERT(weights.
Size() ==
size,
"Different # of weights and dofs.");
1788 MFEM_ASSERT(_lo.
Size() ==
size,
"Different # of lower bounds and dofs.");
1789 MFEM_ASSERT(_hi.
Size() ==
size,
"Different # of upper bounds and dofs.");
1792 double tol = 1.e-12;
1800 slbqp.
Mult(vals, new_vals);
1806 double _min,
double _max)
1811 Vector vals, new_vals(size);
1814 double max_val = vals.
Max();
1815 double min_val = vals.
Min();
1817 if (max_val <= _min)
1824 if (_min <= min_val && max_val <= _max)
1829 Vector minv(size), maxv(size);
1830 minv = (_min > min_val) ? _min : min_val;
1831 maxv = (_max < max_val) ? _max : max_val;
1844 R->
Mult(*
this, tmp);
1845 P->
Mult(tmp, *
this);
1863 for (j = 0; j < vertices.
Size(); j++)
1865 nval(vertices[j]) += values[j];
1866 overlap[vertices[j]]++;
1869 for (i = 0; i < overlap.Size(); i++)
1871 nval(i) /= overlap[i];
1889 for (
int i = 0; i <
fes->
GetNE(); i++)
1897 for (
int j = 0; j < vdofs.
Size(); j++)
1901 MFEM_VERIFY(vals[j] != 0.0,
1902 "Coefficient has zeros, harmonic avg is undefined!");
1903 (*this)(vdofs[j]) += 1.0 / vals[j];
1907 (*this)(vdofs[j]) += vals[j];
1909 else { MFEM_ABORT(
"Not implemented"); }
1911 zones_per_vdof[vdofs[j]]++;
1930 for (
int i = 0; i <
fes->
GetNE(); i++)
1938 for (
int j = 0; j < vdofs.
Size(); j++)
1955 MFEM_VERIFY(vals[j] != 0.0,
1956 "Coefficient has zeros, harmonic avg is undefined!");
1957 (*this)(ldof) += isign / vals[j];
1961 (*this)(ldof) += isign*vals[j];
1964 else { MFEM_ABORT(
"Not implemented"); }
1966 zones_per_vdof[ldof]++;
1975 int i, j, fdof, d, ind, vdim;
1999 for (j = 0; j < fdof; j++)
2003 if (vcoeff) { vcoeff->
Eval(vc, *transf, ip); }
2004 for (d = 0; d < vdim; d++)
2006 if (!vcoeff && !coeff[d]) {
continue; }
2008 val = vcoeff ? vc(d) : coeff[d]->
Eval(*transf, ip);
2009 if ( (ind = vdofs[fdof*d+j]) < 0 )
2011 val = -val, ind = -1-ind;
2013 if (++values_counter[ind] == 1)
2019 (*this)(ind) += val;
2042 for (i = 0; i < bdr_edges.
Size(); i++)
2044 int edge = bdr_edges[i];
2046 if (vdofs.
Size() == 0) {
continue; }
2054 for (d = 0; d < vdim; d++)
2056 if (!coeff[d]) {
continue; }
2058 fe->
Project(*coeff[d], *transf, vals);
2059 for (
int k = 0; k < vals.
Size(); k++)
2061 ind = vdofs[d*vals.
Size()+k];
2062 if (++values_counter[ind] == 1)
2064 (*this)(ind) = vals(k);
2068 (*this)(ind) += vals(k);
2076 fe->
Project(*vcoeff, *transf, vals);
2077 for (
int k = 0; k < vals.
Size(); k++)
2080 if (++values_counter[ind] == 1)
2082 (*this)(ind) = vals(k);
2086 (*this)(ind) += vals(k);
2097 for (
int i = 0; i < dofs.
Size(); i++)
2100 double val = vals(i);
2101 if (k < 0) { k = -1 - k; val = -val; }
2102 if (++values_counter[k] == 1)
2137 fe->
Project(vcoeff, *T, lvec);
2138 accumulate_dofs(dofs, lvec, *
this, values_counter);
2148 for (
int i = 0; i < bdr_edges.
Size(); i++)
2150 int edge = bdr_edges[i];
2152 if (dofs.
Size() == 0) {
continue; }
2158 fe->
Project(vcoeff, *T, lvec);
2159 accumulate_dofs(dofs, lvec, *
this, values_counter);
2169 for (
int i = 0; i <
size; i++)
2171 const int nz = zones_per_vdof[i];
2172 if (nz) { (*this)(i) /= nz; }
2177 for (
int i = 0; i <
size; i++)
2179 const int nz = zones_per_vdof[i];
2180 if (nz) { (*this)(i) = nz/(*
this)(i); }
2185 MFEM_ABORT(
"invalid AvgType");
2200 const double *center = delta_coeff.
Center();
2201 const double *vert = mesh->
GetVertex(0);
2202 double min_dist, dist;
2206 min_dist =
Distance(center, vert, dim);
2207 for (
int i = 0; i < mesh->
GetNV(); i++)
2210 dist =
Distance(center, vert, dim);
2211 if (dist < min_dist)
2221 if (min_dist >= delta_coeff.
Tol())
2230 Vector vals, loc_mass_vals;
2231 for (
int i = 0; i < mesh->
GetNE(); i++)
2234 for (
int j = 0; j < vertices.
Size(); j++)
2235 if (vertices[j] == v_idx)
2245 loc_mass.Mult(vals, loc_mass_vals);
2246 integral += loc_mass_vals.
Sum();
2256 if (delta_c == NULL)
2261 for (
int i = 0; i <
fes->
GetNE(); i++)
2275 (*this) *= (delta_c->
Scale() / integral);
2288 for (
int i = 0; i < dofs.
Size(); i++)
2301 (*this)(vdof) = coeff.
Eval(*T, ip);
2331 for (
int i = 0; i < dofs.
Size(); i++)
2343 vcoeff.
Eval(val, *T, ip);
2347 (*this)(vdof) = val(vd);
2354 int i, j, fdof, d, ind, vdim;
2368 for (j = 0; j < fdof; j++)
2372 for (d = 0; d < vdim; d++)
2374 if (!coeff[d]) {
continue; }
2376 val = coeff[d]->
Eval(*transf, ip);
2377 if ( (ind = vdofs[fdof*d+j]) < 0 )
2379 val = -val, ind = -1-ind;
2399 for (
int i = 0; i <
fes->
GetNE(); i++)
2408 for (
int j = 0; j < vdofs.
Size(); j++)
2410 if (attr > dof_attr[vdofs[j]])
2412 (*this)(vdofs[j]) = vals[j];
2413 dof_attr[vdofs[j]] = attr;
2454 for (
int i = 0; i < values_counter.
Size(); i++)
2456 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2471 for (
int i = 0; i < values_counter.
Size(); i++)
2473 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2489 Vector vc(dim), nor(dim), lvec, shape;
2509 vcoeff.
Eval(vc, *T, ip);
2512 lvec.
Add(ip.
weight * (vc * nor), shape);
2524 Vector vc(dim), nor(dim), lvec;
2540 vcoeff.
Eval(vc, *T, ip);
2542 lvec(j) = (vc * nor);
2559 for (
int i = 0; i < values_counter.
Size(); i++)
2561 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2570 double error = 0.0,
a;
2575 int fdof, d, i, intorder, j, k;
2601 for (k = 0; k < fdof; k++)
2602 if (vdofs[fdof*d+k] >= 0)
2604 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2608 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2611 a -= exsol[d]->
Eval(*transf, ip);
2617 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2630 for (
int i = 0; i <
fes->
GetNE(); i++)
2632 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2634 int intorder = 2*fe->
GetOrder() + 3;
2646 exsol.
Eval(exact_vals, *T, *ir);
2649 vals.
Norm2(loc_errs);
2654 error += ip.
weight * T->
Weight() * (loc_errs(j) * loc_errs(j));
2658 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2673 for (
int i = 0; i <
fes->
GetNE(); i++)
2693 exgrad->
Eval(vec,*Tr,ip);
2698 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2711 int n = (dim == 3) ? dim : 1;
2714 for (
int i = 0; i <
fes->
GetNE(); i++)
2734 excurl->
Eval(vec,*Tr,ip);
2740 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2746 double error = 0.0,
a;
2752 for (
int i = 0; i <
fes->
GetNE(); i++)
2776 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2783 int fdof,
dim, intorder, k;
2788 Vector shape, el_dofs, err_val, ell_coeff_val;
2799 int i1 = face_elem_transf->
Elem1No;
2800 int i2 = face_elem_transf->
Elem2No;
2807 intorder = 2 * intorder;
2820 transf = face_elem_transf->
Elem1;
2826 for (k = 0; k < fdof; k++)
2829 el_dofs(k) = (*this)(vdofs[k]);
2833 el_dofs(k) = - (*this)(-1-vdofs[k]);
2840 ell_coeff_val(j) = ell_coeff->
Eval(*transf, eip);
2841 err_val(j) = exsol->
Eval(*transf, eip) - (shape * el_dofs);
2847 transf = face_elem_transf->
Elem2;
2853 for (k = 0; k < fdof; k++)
2856 el_dofs(k) = (*this)(vdofs[k]);
2860 el_dofs(k) = - (*this)(-1-vdofs[k]);
2867 ell_coeff_val(j) += ell_coeff->
Eval(*transf, eip);
2868 ell_coeff_val(j) *= 0.5;
2869 err_val(j) -= (exsol->
Eval(*transf, eip) - (shape * el_dofs));
2873 transf = face_elem_transf;
2878 error += (ip.
weight * Nu * ell_coeff_val(j) *
2879 pow(transf->
Weight(), 1.0-1.0/(dim-1)) *
2880 err_val(j) * err_val(j));
2884 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2890 int norm_type)
const
2892 double error1 = 0.0;
2893 double error2 = 0.0;
2897 return sqrt(error1 * error1 + error2 * error2);
2906 return sqrt(L2error*L2error + GradError*GradError);
2915 return sqrt(L2error*L2error + DivError*DivError);
2924 return sqrt(L2error*L2error + CurlError*CurlError);
2930 double error = 0.0,
a;
2935 int fdof, d, i, intorder, j, k;
2962 for (k = 0; k < fdof; k++)
2963 if (vdofs[fdof*d+k] >= 0)
2965 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2969 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2971 a -= exsol[d]->
Eval(*transf, ip);
2988 int i, fdof,
dim, intorder, j, k;
2992 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
2995 double a, error = 0.0;
3004 for (i = 0; i < mesh->
GetNE(); i++)
3006 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3023 for (k = 0; k < fdof; k++)
3026 el_dofs(k) = (*this)(vdofs[k]);
3030 el_dofs(k) = -(*this)(-1-vdofs[k]);
3037 a = (el_dofs * shape) - (exsol->
Eval(*transf, ip));
3043 for (i = 0; i < mesh->
GetNE(); i++)
3045 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3063 for (k = 0; k < fdof; k++)
3066 el_dofs(k) = (*this)(vdofs[k]);
3070 el_dofs(k) = -(*this)(-1-vdofs[k]);
3077 exgrad->
Eval(e_grad, *transf, ip);
3079 Mult(dshape, Jinv, dshapet);
3098 for (
int i = 0; i <
fes->
GetNE(); i++)
3108 int intorder = 2*fe->
GetOrder() + 3;
3117 double err = fabs(vals(j) - exsol.
Eval(*T, ip));
3123 err *= weight->
Eval(*T, ip);
3131 err *= weight->
Eval(*T, ip);
3133 error = std::max(error, err);
3143 error = -pow(-error, 1./p);
3147 error = pow(error, 1./p);
3160 "Incorrect size for result vector");
3167 for (
int i = 0; i <
fes->
GetNE(); i++)
3177 int intorder = 2*fe->
GetOrder() + 3;
3186 double err = fabs(vals(j) - exsol.
Eval(*T, ip));
3192 err *= weight->
Eval(*T, ip);
3200 err *= weight->
Eval(*T, ip);
3202 error[i] = std::max(error[i], err);
3210 error[i] = -pow(-error[i], 1./p);
3214 error[i] = pow(error[i], 1./p);
3231 for (
int i = 0; i <
fes->
GetNE(); i++)
3241 int intorder = 2*fe->
GetOrder() + 3;
3246 exsol.
Eval(exact_vals, *T, *ir);
3253 vals.
Norm2(loc_errs);
3257 v_weight->
Eval(exact_vals, *T, *ir);
3260 for (
int j = 0; j < vals.
Width(); j++)
3263 for (
int d = 0; d < vals.
Height(); d++)
3265 err += vals(d,j)*exact_vals(d,j);
3267 loc_errs(j) = fabs(err);
3274 double err = loc_errs(j);
3280 err *= weight->
Eval(*T, ip);
3288 err *= weight->
Eval(*T, ip);
3290 error = std::max(error, err);
3300 error = -pow(-error, 1./p);
3304 error = pow(error, 1./p);
3319 "Incorrect size for result vector");
3327 for (
int i = 0; i <
fes->
GetNE(); i++)
3337 int intorder = 2*fe->
GetOrder() + 3;
3342 exsol.
Eval(exact_vals, *T, *ir);
3349 vals.
Norm2(loc_errs);
3353 v_weight->
Eval(exact_vals, *T, *ir);
3356 for (
int j = 0; j < vals.
Width(); j++)
3359 for (
int d = 0; d < vals.
Height(); d++)
3361 err += vals(d,j)*exact_vals(d,j);
3363 loc_errs(j) = fabs(err);
3370 double err = loc_errs(j);
3376 err *= weight->
Eval(*T, ip);
3384 err *= weight->
Eval(*T, ip);
3386 error[i] = std::max(error[i], err);
3394 error[i] = -pow(-error[i], 1./p);
3398 error[i] = pow(error[i], 1./p);
3425 out <<
"NURBS_patches\n";
3442 #ifdef MFEM_USE_ADIOS2
3444 const std::string& variable_name,
3447 out.
Save(*
this, variable_name, type);
3463 out <<
"SCALARS " << field_name <<
" double 1\n"
3464 <<
"LOOKUP_TABLE default\n";
3465 for (
int i = 0; i < mesh->
GetNE(); i++)
3472 for (
int j = 0; j < val.
Size(); j++)
3474 out << val(j) <<
'\n';
3478 else if ( (vec_dim == 2 || vec_dim == 3) && mesh->
SpaceDimension() > 1)
3481 out <<
"VECTORS " << field_name <<
" double\n";
3482 for (
int i = 0; i < mesh->
GetNE(); i++)
3491 for (
int j = 0; j < vval.
Width(); j++)
3493 out << vval(0, j) <<
' ' << vval(1, j) <<
' ';
3509 for (
int vd = 0; vd < vec_dim; vd++)
3511 out <<
"SCALARS " << field_name << vd <<
" double 1\n"
3512 <<
"LOOKUP_TABLE default\n";
3513 for (
int i = 0; i < mesh->
GetNE(); i++)
3520 for (
int j = 0; j < val.
Size(); j++)
3522 out << val(j) <<
'\n';
3533 double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
3534 double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
3535 double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
3536 v1[2] * v2[0] - v1[0] * v2[2],
3537 v1[0] * v2[1] - v1[1] * v2[0]
3539 double rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
3540 n[0] *= rl; n[1] *= rl; n[2] *= rl;
3542 out <<
" facet normal " << n[0] <<
' ' << n[1] <<
' ' << n[2]
3544 <<
"\n vertex " << p1[0] <<
' ' << p1[1] <<
' ' << p1[2]
3545 <<
"\n vertex " << p2[0] <<
' ' << p2[1] <<
' ' << p2[2]
3546 <<
"\n vertex " << p3[0] <<
' ' << p3[1] <<
' ' << p3[2]
3547 <<
"\n endloop\n endfacet\n";
3563 double pts[4][3], bbox[3][2];
3565 out <<
"solid GridFunction\n";
3567 bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
3568 bbox[2][0] = bbox[2][1] = 0.0;
3569 for (i = 0; i < mesh->
GetNE(); i++)
3576 for (k = 0; k < RG.
Size()/n; k++)
3578 for (j = 0; j < n; j++)
3581 pts[j][0] = pointmat(0,l);
3582 pts[j][1] = pointmat(1,l);
3583 pts[j][2] = values(l);
3599 bbox[0][0] = pointmat(0,0);
3600 bbox[0][1] = pointmat(0,0);
3601 bbox[1][0] = pointmat(1,0);
3602 bbox[1][1] = pointmat(1,0);
3603 bbox[2][0] = values(0);
3604 bbox[2][1] = values(0);
3607 for (j = 0; j < values.
Size(); j++)
3609 if (bbox[0][0] > pointmat(0,j))
3611 bbox[0][0] = pointmat(0,j);
3613 if (bbox[0][1] < pointmat(0,j))
3615 bbox[0][1] = pointmat(0,j);
3617 if (bbox[1][0] > pointmat(1,j))
3619 bbox[1][0] = pointmat(1,j);
3621 if (bbox[1][1] < pointmat(1,j))
3623 bbox[1][1] = pointmat(1,j);
3625 if (bbox[2][0] > values(j))
3627 bbox[2][0] = values(j);
3629 if (bbox[2][1] < values(j))
3631 bbox[2][1] = values(j);
3636 mfem::out <<
"[xmin,xmax] = [" << bbox[0][0] <<
',' << bbox[0][1] <<
"]\n"
3637 <<
"[ymin,ymax] = [" << bbox[1][0] <<
',' << bbox[1][1] <<
"]\n"
3638 <<
"[zmin,zmax] = [" << bbox[2][0] <<
',' << bbox[2][1] <<
']'
3641 out <<
"endsolid GridFunction" << endl;
3653 const char *msg =
"invalid input stream";
3659 in >> ident; MFEM_VERIFY(ident ==
"VDim:", msg);
3686 out <<
"VDim: " <<
vdim <<
'\n'
3703 int with_subdomains,
3711 int nfe = ufes->
GetNE();
3715 Vector ul, fl, fla, d_xyz;
3725 if (with_subdomains)
3730 double total_error = 0.0;
3731 for (
int s = 1; s <= nsd; s++)
3734 u.
ComputeFlux(blfi, flux, with_coeff, (with_subdomains ? s : -1));
3736 for (
int i = 0; i < nfe; i++)
3738 if (with_subdomains && ufes->
GetAttribute(i) != s) {
continue; }
3748 *ffes->
GetFE(i), fl, with_coeff);
3753 (aniso_flags ? &d_xyz : NULL));
3755 error_estimates(i) = std::sqrt(err);
3761 for (
int k = 0; k <
dim; k++)
3766 double thresh = 0.15 * 3.0/
dim;
3768 for (
int k = 0; k <
dim; k++)
3770 if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
3773 (*aniso_flags)[i] = flag;
3778 return std::sqrt(total_error);
3801 for (
int j = 0; j < nip; j++)
3818 norm = std::max(norm, err);
3827 norm = -pow(-norm, 1./p);
3831 norm = pow(norm, 1./p);
3845 return sol_in.
Eval(*T_in, ip);
3856 string cname = name;
3857 if (cname ==
"Linear")
3861 else if (cname ==
"Quadratic")
3865 else if (cname ==
"Cubic")
3869 else if (!strncmp(name,
"H1_", 3))
3873 else if (!strncmp(name,
"H1Pos_", 6))
3878 else if (!strncmp(name,
"L2_T", 4))
3882 else if (!strncmp(name,
"L2_", 3))
3888 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.
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().
double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
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.
virtual double ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, double Nu, const IntegrationRule *irs[]=NULL) const
Returns the Face Jumps error for L2 elements.
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 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)
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
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.
void MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of fec and fes.
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.
void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
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.
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
void ImposeBounds(int i, const Vector &weights, const Vector &_lo, const Vector &_hi)
int Size() const
Returns the size of the vector.
int GetBdrAttribute(int i) const
int GetNE() const
Returns number of elements.
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)
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...
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
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.
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Geometry::Type GetElementBaseGeometry(int i) const
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
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
virtual void GetFaceDofs(int i, Array< int > &dofs) const
Returns the indexes of the degrees of freedom for i'th face including the dofs for the edges and the ...
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.
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)
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
void SetLinearConstraint(const Vector &_w, double _a)
const FiniteElement * GetEdgeElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
void SetPrintLevel(int print_lvl)
void GetHessians(int i, const IntegrationRule &ir, DenseMatrix &hess, int vdim=1) const
bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
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.
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
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. If all dofs are true, then tv will be set to point to th...
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.
int GetVDim() const
Returns vector dimension.
A class for non-conforming AMR on higher-order hexahedral, prismatic, quadrilateral or triangular mes...
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
void GetEdgeDofs(int i, Array< int > &dofs) const
Returns the indexes of the degrees of freedom for i'th edge including the dofs for the vertices of th...
QuadratureFunction & operator=(double value)
Redefine '=' for QuadratureFunction = constant.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
int SpaceDimension() const
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
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.
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...
void SetBounds(const Vector &_lo, const Vector &_hi)
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)
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i'th boundary element.
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.
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
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...
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)
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).
double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
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.
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 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[])
Class representing the storage layout of a QuadratureFunction.
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.
long GetSequence() const
Return update counter (see Mesh::sequence)
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.
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.
double f(const Vector &p)
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const