16 #include "../mesh/nurbs.hpp" 17 #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, l_nvdofs*
sizeof(
double)*vdim);
134 l_data += vdim*l_nvdofs;
135 g_data += vdim*g_nvdofs;
136 memcpy(g_data+vdim*ei, l_data, l_nedofs*
sizeof(
double)*vdim);
137 l_data += vdim*l_nedofs;
138 g_data += vdim*g_nedofs;
139 memcpy(g_data+vdim*fi, l_data, l_nfdofs*
sizeof(
double)*vdim);
140 l_data += vdim*l_nfdofs;
141 g_data += vdim*g_nfdofs;
142 memcpy(g_data+vdim*di, l_data, l_nddofs*
sizeof(
double)*vdim);
143 l_data += vdim*l_nddofs;
144 g_data += vdim*g_nddofs;
182 old_data.
Swap(*
this);
185 T->Mult(old_data, *
this);
213 MFEM_ASSERT(v.
Size() >= v_offset +
f->GetVSize(),
"");
245 MFEM_ASSERT(tv.
Size() >= tv_offset +
f->GetTrueVSize(),
"");
266 int nfe = ufes->
GetNE();
274 for (
int i = 0; i < nfe; i++)
276 if (subdomain >= 0 && ufes->
GetAttribute(i) != subdomain)
284 u.GetSubVector(udofs, ul);
292 *ffes->
GetFE(i), fl, wcoef);
301 for (
int j = 0; j < fdofs.
Size(); j++)
317 for (
int i = 0; i < count.Size(); i++)
319 if (count[i] != 0) { flux(i) /= count[i]; }
403 int dof = FElem->
GetDof();
419 for (
int k = 0; k < n; k++)
422 nval[k] = shape * (&loc_data[dof * vdim]);
428 for (
int k = 0; k < n; k++)
432 nval[k] = shape * (&loc_data[dof * vdim]);
440 for (
int k = 0; k < n; k++)
444 nval[k] = loc_data * (&vshape(0,vdim));
465 fe->CalcPhysShape(*Tr, DofVal);
473 return (DofVal * LocVec);
480 int dof = FElem->
GetDof();
504 for (
int k = 0; k < vdim; k++)
506 val(k) = shape * (&loc_data[dof * k]);
531 int dof = FElem->
GetDof();
532 Vector DofVal(dof), loc_data(dof);
540 for (
int k = 0; k < n; k++)
543 vals(k) = DofVal * loc_data;
549 for (
int k = 0; k < n; k++)
553 vals(k) = DofVal * loc_data;
582 "invalid FE map type");
584 int dof = FElem->
GetDof();
585 Vector DofLap(dof), loc_data(dof);
587 for (
int k = 0; k < n; k++)
592 laps(k) = DofLap * loc_data;
625 "invalid FE map type");
627 int dof = FElem->
GetDof();
635 for (
int k = 0; k < n; k++)
641 for (
int j = 0; j <
size; j++)
643 for (
int d = 0; d < dof; d++)
645 hess(k,j) += DofHes(d,j) * loc_data[d];
730 fip.
x = 1.0 - ip.
x - ip.
y;
736 fip.
y = 1.0 - ip.
x - ip.
y;
775 int comp,
Vector *tr)
const 780 T.Transform(ip, *tr);
786 switch (T.ElementType)
801 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \"" 803 <<
"on mesh edges.");
816 MFEM_ABORT(
"GridFunction::GetValue: Field continuity type \"" 818 <<
"on mesh faces.");
870 MFEM_ABORT(
"GridFunction::GetValue: Unsupported element type \"" 871 << T.ElementType <<
"\"");
888 return (DofVal * LocVec);
898 T.Transform(ir, *tr);
903 for (
int j = 0; j < nip; j++)
918 T.Transform(ip, *tr);
925 switch (T.ElementType)
940 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \"" 942 <<
"on mesh edges.");
955 MFEM_ABORT(
"GridFunction::GetVectorValue: Field continuity type \"" 957 <<
"on mesh faces.");
1009 MFEM_ABORT(
"GridFunction::GetVectorValue: Unsupported element type \"" 1010 << T.ElementType <<
"\"");
1011 if (val.
Size() > 0) { val = NAN; }
1036 for (
int k = 0; k < vdim; k++)
1038 val(k) = shape * (&loc_data[dof * k]);
1044 int vdim = std::max(spaceDim, fe->
GetVDim());
1048 vshape.MultTranspose(loc_data, val);
1059 T.Transform(ir, *tr);
1063 int dof = FElem->
GetDof();
1081 for (
int j = 0; j < nip; j++)
1087 for (
int k = 0; k < vdim; k++)
1089 vals(k,j) = shape * (&loc_data[dof * k]);
1096 int vdim = std::max(spaceDim, FElem->
GetVDim());
1102 for (
int j = 0; j < nip; j++)
1109 vshape.MultTranspose(loc_data, val_j);
1165 Vector shape, loc_values, orig_loc_values;
1166 int i, j, d, ne, dof, odof, vdim;
1170 for (i = 0; i < ne; i++)
1182 odof = orig_fe->
GetDof();
1183 loc_values.
SetSize(dof * vdim);
1186 for (j = 0; j < dof; j++)
1190 for (d = 0; d < vdim; d++)
1192 loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
1211 Vector shape, loc_values, loc_values_t, orig_loc_values, orig_loc_values_t;
1212 int i, j, d, nbe, dof, odof, vdim;
1216 for (i = 0; i < nbe; i++)
1224 odof = orig_fe->
GetDof();
1225 loc_values.
SetSize(dof * vdim);
1228 for (j = 0; j < dof; j++)
1232 for (d = 0; d < vdim; d++)
1234 loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
1248 int d, k, n, sdim, dof;
1260 Vector loc_data, val(sdim);
1266 for (k = 0; k < n; k++)
1272 for (d = 0; d < sdim; d++)
1289 double *temp =
new double[
size];
1292 for (j = 0; j < ndofs; j++)
1293 for (i = 0; i < vdim; i++)
1295 temp[j+i*ndofs] =
data[k++];
1298 for (i = 0; i <
size; i++)
1326 val(vertices[k]) += vals(k, comp);
1327 overlap[vertices[k]]++;
1331 for (i = 0; i < overlap.Size(); i++)
1333 val(i) /= overlap[i];
1341 int d, i, k, ind, dof, sdim;
1350 for (i = 0; i < new_fes->
GetNE(); i++)
1357 for (d = 0; d < sdim; d++)
1359 for (k = 0; k < dof; k++)
1361 if ( (ind=new_vdofs[dof*d+k]) < 0 )
1363 ind = -1-ind, vals(k, d) = - vals(k, d);
1365 vec_field(ind) += vals(k, d);
1371 for (i = 0; i < overlap.Size(); i++)
1373 vec_field(i) /= overlap[i];
1386 Vector pt_grad, loc_func;
1387 int i, j, k,
dim, dof, der_dof, ind;
1394 for (i = 0; i < der_fes->
GetNE(); i++)
1403 der_dof = der_fe->
GetDof();
1409 for (j = 0; j < dof; j++)
1410 loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
1412 for (k = 0; k < der_dof; k++)
1420 for (j = 0; j <
dim; j++)
1422 a += inv_jac(j, der_comp) * pt_grad(j);
1424 der(der_dofs[k]) +=
a;
1425 zones_per_dof[der_dofs[k]]++;
1435 for (
int i = 0; i < overlap.
Size(); i++)
1437 der(i) /= overlap[i];
1454 MultAtB(loc_data_mat, dshape, gh);
1459 switch (T.ElementType)
1463 int elNo = T.ElementNo;
1468 "invalid FE map type");
1473 for (
int i = 0; i < Jinv.
Width(); i++)
1475 for (
int j = 0; j < Jinv.
Height(); j++)
1477 div_v += grad_hat(i, j) * Jinv(j, i);
1494 return (loc_data * divshape) / T.Weight();
1540 MFEM_ABORT(
"GridFunction::GetDivergence: Unsupported element type \"" 1541 << T.ElementType <<
"\"");
1549 switch (T.ElementType)
1553 int elNo = T.ElementNo;
1558 "invalid FE map type");
1564 Mult(grad_hat, Jinv, grad);
1565 MFEM_ASSERT(grad.Height() == grad.Width(),
"");
1566 if (grad.Height() == 3)
1569 curl(0) = grad(2,1) - grad(1,2);
1570 curl(1) = grad(0,2) - grad(2,0);
1571 curl(2) = grad(1,0) - grad(0,1);
1573 else if (grad.Height() == 2)
1576 curl(0) = grad(1,0) - grad(0,1);
1591 curl.
SetSize(curl_shape.Width());
1593 curl_shape.MultTranspose(loc_data, curl);
1637 MFEM_ABORT(
"GridFunction::GetCurl: Unsupported element type \"" 1638 << T.ElementType <<
"\"");
1645 switch (T.ElementType)
1651 "invalid FE map type");
1652 MFEM_ASSERT(
fes->
GetVDim() == 1,
"Defined for scalar functions.");
1662 T.InverseJacobian().MultTranspose(gh, grad);
1706 MFEM_ABORT(
"GridFunction::GetGradient: Unsupported element type \"" 1707 << T.ElementType <<
"\"");
1728 dshape.MultTranspose(lval, gh);
1739 switch (T.ElementType)
1749 Mult(grad_hat, Jinv, grad);
1792 MFEM_ABORT(
"GridFunction::GetVectorGradient: " 1793 "Unsupported element type \"" << T.ElementType <<
"\"");
1805 Vector loc_avgs, loc_this;
1810 for (
int i = 0; i <
fes->
GetNE(); i++)
1815 te_doftrans = avgs.FESpace()->GetElementDofs(i, te_dofs);
1822 loc_mass.
Mult(loc_this, loc_avgs);
1827 avgs.AddElementVector(te_dofs, loc_avgs);
1829 loc_mass.
Mult(loc_this, loc_avgs);
1830 int_psi.AddElementVector(te_dofs, loc_avgs);
1832 for (
int i = 0; i < avgs.Size(); i++)
1834 avgs(i) /= int_psi(i);
1855 if (!mesh->
GetNE()) {
return; }
1866 MFEM_VERIFY(vdim == src.
fes->
GetVDim(),
"incompatible vector dimensions!");
1871 for (
int i = 0; i < mesh->
GetNE(); i++)
1878 dest_lvec.SetSize(vdim*P.
Height());
1888 for (
int vd = 0; vd < vdim; vd++)
1915 MFEM_ASSERT(weights.
Size() ==
size,
"Different # of weights and dofs.");
1916 MFEM_ASSERT(lo_.
Size() ==
size,
"Different # of lower bounds and dofs.");
1917 MFEM_ASSERT(hi_.
Size() ==
size,
"Different # of upper bounds and dofs.");
1920 double tol = 1.e-12;
1928 slbqp.
Mult(vals, new_vals);
1938 double min_,
double max_)
1950 double max_val = vals.
Max();
1951 double min_val = vals.
Min();
1953 if (max_val <= min_)
1964 if (min_ <= min_val && max_val <= max_)
1970 minv = (min_ > min_val) ? min_ : min_val;
1971 maxv = (max_ < max_val) ? max_ : max_val;
1984 R->
Mult(*
this, tmp);
1985 P->
Mult(tmp, *
this);
2003 for (j = 0; j < vertices.
Size(); j++)
2005 nval(vertices[j]) += values[j];
2006 overlap[vertices[j]]++;
2009 for (i = 0; i < overlap.Size(); i++)
2011 nval(i) /= overlap[i];
2022 for (
int i = 0; i <
fes->
GetNE(); i++)
2026 for (
int j = 0; j < vdofs.
Size(); j++)
2028 elem_per_vdof[vdofs[j]]++;
2047 for (
int i = 0; i <
fes->
GetNE(); i++)
2055 for (
int j = 0; j < vdofs.
Size(); j++)
2059 MFEM_VERIFY(vals[j] != 0.0,
2060 "Coefficient has zeros, harmonic avg is undefined!");
2061 (*this)(vdofs[j]) += 1.0 / vals[j];
2065 (*this)(vdofs[j]) += vals[j];
2067 else { MFEM_ABORT(
"Not implemented"); }
2069 zones_per_vdof[vdofs[j]]++;
2088 for (
int i = 0; i <
fes->
GetNE(); i++)
2096 for (
int j = 0; j < vdofs.
Size(); j++)
2113 MFEM_VERIFY(vals[j] != 0.0,
2114 "Coefficient has zeros, harmonic avg is undefined!");
2115 (*this)(ldof) += isign / vals[j];
2119 (*this)(ldof) += isign*vals[j];
2122 else { MFEM_ABORT(
"Not implemented"); }
2124 zones_per_vdof[ldof]++;
2133 int i, j, fdof, d, ind, vdim;
2157 for (j = 0; j < fdof; j++)
2161 if (vcoeff) { vcoeff->
Eval(vc, *transf, ip); }
2162 for (d = 0; d < vdim; d++)
2164 if (!vcoeff && !coeff[d]) {
continue; }
2166 val = vcoeff ? vc(d) : coeff[d]->
Eval(*transf, ip);
2167 if ( (ind = vdofs[fdof*d+j]) < 0 )
2169 val = -val, ind = -1-ind;
2171 if (++values_counter[ind] == 1)
2177 (*this)(ind) += val;
2200 for (i = 0; i < bdr_edges.
Size(); i++)
2202 int edge = bdr_edges[i];
2204 if (vdofs.
Size() == 0) {
continue; }
2212 for (d = 0; d < vdim; d++)
2214 if (!coeff[d]) {
continue; }
2216 fe->
Project(*coeff[d], *transf, vals);
2217 for (
int k = 0; k < vals.
Size(); k++)
2219 ind = vdofs[d*vals.
Size()+k];
2220 if (++values_counter[ind] == 1)
2222 (*this)(ind) = vals(k);
2226 (*this)(ind) += vals(k);
2234 fe->
Project(*vcoeff, *transf, vals);
2235 for (
int k = 0; k < vals.
Size(); k++)
2238 if (++values_counter[ind] == 1)
2240 (*this)(ind) = vals(k);
2244 (*this)(ind) += vals(k);
2255 for (
int i = 0; i < dofs.
Size(); i++)
2258 double val = vals(i);
2259 if (k < 0) { k = -1 - k; val = -val; }
2260 if (++values_counter[k] == 1)
2295 fe->
Project(vcoeff, *T, lvec);
2297 accumulate_dofs(dofs, lvec, *
this, values_counter);
2307 for (
int i = 0; i < bdr_edges.
Size(); i++)
2309 int edge = bdr_edges[i];
2311 if (dofs.
Size() == 0) {
continue; }
2317 fe->
Project(vcoeff, *T, lvec);
2318 accumulate_dofs(dofs, lvec, *
this, values_counter);
2328 for (
int i = 0; i <
size; i++)
2330 const int nz = zones_per_vdof[i];
2331 if (nz) { (*this)(i) /= nz; }
2336 for (
int i = 0; i <
size; i++)
2338 const int nz = zones_per_vdof[i];
2339 if (nz) { (*this)(i) = nz/(*
this)(i); }
2344 MFEM_ABORT(
"invalid AvgType");
2359 const double *center = delta_coeff.
Center();
2360 const double *vert = mesh->
GetVertex(0);
2361 double min_dist, dist;
2366 for (
int i = 0; i < mesh->
GetNV(); i++)
2370 if (dist < min_dist)
2380 if (min_dist >= delta_coeff.
Tol())
2389 Vector vals, loc_mass_vals;
2390 for (
int i = 0; i < mesh->
GetNE(); i++)
2393 for (
int j = 0; j < vertices.
Size(); j++)
2394 if (vertices[j] == v_idx)
2408 loc_mass.Mult(vals, loc_mass_vals);
2409 integral += loc_mass_vals.
Sum();
2420 if (delta_c == NULL)
2425 for (
int i = 0; i <
fes->
GetNE(); i++)
2443 (*this) *= (delta_c->
Scale() / integral);
2456 for (
int i = 0; i < dofs.
Size(); i++)
2468 T->SetIntPoint(&ip);
2469 (*this)(vdof) = coeff.
Eval(*T, ip);
2505 for (
int i = 0; i < dofs.
Size(); i++)
2516 T->SetIntPoint(&ip);
2517 vcoeff.
Eval(val, *T, ip);
2521 (*this)(vdof) = val(vd);
2554 int i, j, fdof, d, ind, vdim;
2570 for (j = 0; j < fdof; j++)
2574 for (d = 0; d < vdim; d++)
2576 if (!coeff[d]) {
continue; }
2578 val = coeff[d]->
Eval(*transf, ip);
2579 if ( (ind = vdofs[fdof*d+j]) < 0 )
2581 val = -val, ind = -1-ind;
2601 for (
int i = 0; i <
fes->
GetNE(); i++)
2610 for (
int j = 0; j < vdofs.
Size(); j++)
2612 if (attr > dof_attr[vdofs[j]])
2614 (*this)(vdofs[j]) = vals[j];
2615 dof_attr[vdofs[j]] = attr;
2657 for (
int i = 0; i < values_counter.
Size(); i++)
2659 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2674 ess_vdofs_marker = 0;
2678 if (!coeff[i]) {
continue; }
2680 for (
int j = 0; j<
Size(); j++)
2682 ess_vdofs_marker[j] = bool(ess_vdofs_marker[j]) ||
2683 bool(component_dof_marker[j]);
2686 for (
int i = 0; i < values_counter.
Size(); i++)
2688 MFEM_ASSERT(
bool(values_counter[i]) == ess_vdofs_marker[i],
2723 T->SetIntPoint(&ip);
2724 vcoeff.
Eval(vc, *T, ip);
2727 lvec.
Add(ip.
weight * (vc * nor), shape);
2754 T->SetIntPoint(&ip);
2755 vcoeff.
Eval(vc, *T, ip);
2757 lvec(j) = (vc * nor);
2778 for (
int i = 0; i < values_counter.
Size(); i++)
2780 MFEM_ASSERT(
bool(values_counter[i]) ==
bool(ess_vdofs_marker[i]),
2790 double error = 0.0,
a;
2795 int fdof, d, i, intorder, j, k;
2799 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2822 for (k = 0; k < fdof; k++)
2823 if (vdofs[fdof*d+k] >= 0)
2825 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2829 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2832 a -= exsol[d]->
Eval(*transf, ip);
2838 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2851 for (
int i = 0; i <
fes->
GetNE(); i++)
2853 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2855 int intorder = 2*fe->
GetOrder() + 3;
2867 exsol.
Eval(exact_vals, *T, *ir);
2870 vals.
Norm2(loc_errs);
2874 T->SetIntPoint(&ip);
2875 error += ip.
weight * T->Weight() * (loc_errs(j) * loc_errs(j));
2879 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2913 exgrad->
Eval(vec,*Tr,ip);
2917 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2932 for (
int i = 0; i <
fes->
GetNE(); i++)
2952 exgrad->
Eval(vec,*Tr,ip);
2957 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2972 for (
int i = 0; i <
fes->
GetNE(); i++)
2992 excurl->
Eval(vec,*Tr,ip);
2998 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
3004 double error = 0.0,
a;
3010 for (
int i = 0; i <
fes->
GetNE(); i++)
3034 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
3042 int fdof, intorder, k;
3047 Vector shape, el_dofs, err_val, ell_coeff_val;
3069 intorder = 2 * intorder;
3083 transf = face_elem_transf->
Elem1;
3089 for (k = 0; k < fdof; k++)
3092 el_dofs(k) = (*this)(vdofs[k]);
3096 el_dofs(k) = - (*this)(-1-vdofs[k]);
3103 ell_coeff_val(j) = ell_coeff->
Eval(*transf, eip);
3104 err_val(j) = exsol->
Eval(*transf, eip) - (shape * el_dofs);
3110 transf = face_elem_transf->
Elem2;
3116 for (k = 0; k < fdof; k++)
3119 el_dofs(k) = (*this)(vdofs[k]);
3123 el_dofs(k) = - (*this)(-1-vdofs[k]);
3130 ell_coeff_val(j) += ell_coeff->
Eval(*transf, eip);
3131 ell_coeff_val(j) *= 0.5;
3132 err_val(j) -= (exsol->
Eval(*transf, eip) - (shape * el_dofs));
3136 transf = face_elem_transf;
3141 double nu = jump_scaling.
Eval(h,
p);
3142 error += (ip.
weight * nu * ell_coeff_val(j) *
3144 err_val(j) * err_val(j));
3148 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
3163 int norm_type)
const 3165 double error1 = 0.0;
3166 double error2 = 0.0;
3174 return sqrt(error1 * error1 + error2 * error2);
3183 return sqrt(L2error*L2error + GradError*GradError);
3192 return sqrt(L2error*L2error + DivError*DivError);
3201 return sqrt(L2error*L2error + CurlError*CurlError);
3207 double error = 0.0,
a;
3212 int fdof, d, i, intorder, j, k;
3239 for (k = 0; k < fdof; k++)
3240 if (vdofs[fdof*d+k] >= 0)
3242 a += (*this)(vdofs[fdof*d+k]) * shape(k);
3246 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
3248 a -= exsol[d]->
Eval(*transf, ip);
3265 int i, fdof,
dim, intorder, j, k;
3269 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
3272 double a, error = 0.0;
3281 for (i = 0; i < mesh->
GetNE(); i++)
3283 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3300 for (k = 0; k < fdof; k++)
3303 el_dofs(k) = (*this)(vdofs[k]);
3307 el_dofs(k) = -(*this)(-1-vdofs[k]);
3314 a = (el_dofs * shape) - (exsol->
Eval(*transf, ip));
3320 for (i = 0; i < mesh->
GetNE(); i++)
3322 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3340 for (k = 0; k < fdof; k++)
3343 el_dofs(k) = (*this)(vdofs[k]);
3347 el_dofs(k) = -(*this)(-1-vdofs[k]);
3354 exgrad->
Eval(e_grad, *transf, ip);
3356 Mult(dshape, Jinv, dshapet);
3376 for (
int i = 0; i <
fes->
GetNE(); i++)
3378 if (elems != NULL && (*elems)[i] == 0) {
continue; }
3387 int intorder = 2*fe->
GetOrder() + 3;
3395 T->SetIntPoint(&ip);
3396 double diff = fabs(vals(j) - exsol.
Eval(*T, ip));
3399 diff = pow(diff,
p);
3402 diff *= weight->
Eval(*T, ip);
3404 error += ip.
weight * T->Weight() * diff;
3410 diff *= weight->
Eval(*T, ip);
3412 error = std::max(error, diff);
3422 error = -pow(-error, 1./
p);
3426 error = pow(error, 1./
p);
3439 "Incorrect size for result vector");
3446 for (
int i = 0; i <
fes->
GetNE(); i++)
3456 int intorder = 2*fe->
GetOrder() + 3;
3464 T->SetIntPoint(&ip);
3465 double diff = fabs(vals(j) - exsol.
Eval(*T, ip));
3468 diff = pow(diff,
p);
3471 diff *= weight->
Eval(*T, ip);
3473 error[i] += ip.
weight * T->Weight() * diff;
3479 diff *= weight->
Eval(*T, ip);
3481 error[i] = std::max(error[i], diff);
3489 error[i] = -pow(-error[i], 1./
p);
3493 error[i] = pow(error[i], 1./
p);
3510 for (
int i = 0; i <
fes->
GetNE(); i++)
3520 int intorder = 2*fe->
GetOrder() + 3;
3525 exsol.
Eval(exact_vals, *T, *ir);
3532 vals.
Norm2(loc_errs);
3536 v_weight->
Eval(exact_vals, *T, *ir);
3539 for (
int j = 0; j < vals.
Width(); j++)
3542 for (
int d = 0; d < vals.
Height(); d++)
3544 errj += vals(d,j)*exact_vals(d,j);
3546 loc_errs(j) = fabs(errj);
3552 T->SetIntPoint(&ip);
3553 double errj = loc_errs(j);
3556 errj = pow(errj,
p);
3559 errj *= weight->
Eval(*T, ip);
3561 error += ip.
weight * T->Weight() * errj;
3567 errj *= weight->
Eval(*T, ip);
3569 error = std::max(error, errj);
3579 error = -pow(-error, 1./
p);
3583 error = pow(error, 1./
p);
3598 "Incorrect size for result vector");
3606 for (
int i = 0; i <
fes->
GetNE(); i++)
3616 int intorder = 2*fe->
GetOrder() + 3;
3621 exsol.
Eval(exact_vals, *T, *ir);
3628 vals.
Norm2(loc_errs);
3632 v_weight->
Eval(exact_vals, *T, *ir);
3635 for (
int j = 0; j < vals.
Width(); j++)
3638 for (
int d = 0; d < vals.
Height(); d++)
3640 errj += vals(d,j)*exact_vals(d,j);
3642 loc_errs(j) = fabs(errj);
3648 T->SetIntPoint(&ip);
3649 double errj = loc_errs(j);
3652 errj = pow(errj,
p);
3655 errj *= weight->
Eval(*T, ip);
3657 error[i] += ip.
weight * T->Weight() * errj;
3663 errj *= weight->
Eval(*T, ip);
3665 error[i] = std::max(error[i], errj);
3673 error[i] = -pow(-error[i], 1./
p);
3677 error[i] = pow(error[i], 1./
p);
3704 os <<
"NURBS_patches\n";
3723 ofstream ofs(fname);
3724 ofs.precision(precision);
3728 #ifdef MFEM_USE_ADIOS2 3730 const std::string& variable_name,
3733 os.
Save(*
this, variable_name, type);
3749 os <<
"SCALARS " << field_name <<
" double 1\n" 3750 <<
"LOOKUP_TABLE default\n";
3751 for (
int i = 0; i < mesh->
GetNE(); i++)
3758 for (
int j = 0; j < val.
Size(); j++)
3760 os << val(j) <<
'\n';
3764 else if ( (vec_dim == 2 || vec_dim == 3) && mesh->
SpaceDimension() > 1)
3767 os <<
"VECTORS " << field_name <<
" double\n";
3768 for (
int i = 0; i < mesh->
GetNE(); i++)
3777 for (
int j = 0; j < vval.
Width(); j++)
3779 os << vval(0, j) <<
' ' << vval(1, j) <<
' ';
3795 for (
int vd = 0; vd < vec_dim; vd++)
3797 os <<
"SCALARS " << field_name << vd <<
" double 1\n" 3798 <<
"LOOKUP_TABLE default\n";
3799 for (
int i = 0; i < mesh->
GetNE(); i++)
3806 for (
int j = 0; j < val.
Size(); j++)
3808 os << val(j) <<
'\n';
3819 double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
3820 double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
3821 double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
3822 v1[2] * v2[0] - v1[0] * v2[2],
3823 v1[0] * v2[1] - v1[1] * v2[0]
3825 double rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
3826 n[0] *= rl; n[1] *= rl; n[2] *= rl;
3828 os <<
" facet normal " << n[0] <<
' ' << n[1] <<
' ' << n[2]
3830 <<
"\n vertex " << p1[0] <<
' ' << p1[1] <<
' ' << p1[2]
3831 <<
"\n vertex " << p2[0] <<
' ' << p2[1] <<
' ' << p2[2]
3832 <<
"\n vertex " << p3[0] <<
' ' << p3[1] <<
' ' << p3[2]
3833 <<
"\n endloop\n endfacet\n";
3849 double pts[4][3], bbox[3][2];
3851 os <<
"solid GridFunction\n";
3853 bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
3854 bbox[2][0] = bbox[2][1] = 0.0;
3855 for (i = 0; i < mesh->
GetNE(); i++)
3862 for (k = 0; k < RG.
Size()/n; k++)
3864 for (j = 0; j < n; j++)
3867 pts[j][0] = pointmat(0,l);
3868 pts[j][1] = pointmat(1,l);
3869 pts[j][2] = values(l);
3885 bbox[0][0] = pointmat(0,0);
3886 bbox[0][1] = pointmat(0,0);
3887 bbox[1][0] = pointmat(1,0);
3888 bbox[1][1] = pointmat(1,0);
3889 bbox[2][0] = values(0);
3890 bbox[2][1] = values(0);
3893 for (j = 0; j < values.
Size(); j++)
3895 if (bbox[0][0] > pointmat(0,j))
3897 bbox[0][0] = pointmat(0,j);
3899 if (bbox[0][1] < pointmat(0,j))
3901 bbox[0][1] = pointmat(0,j);
3903 if (bbox[1][0] > pointmat(1,j))
3905 bbox[1][0] = pointmat(1,j);
3907 if (bbox[1][1] < pointmat(1,j))
3909 bbox[1][1] = pointmat(1,j);
3911 if (bbox[2][0] > values(j))
3913 bbox[2][0] = values(j);
3915 if (bbox[2][1] < values(j))
3917 bbox[2][1] = values(j);
3922 mfem::out <<
"[xmin,xmax] = [" << bbox[0][0] <<
',' << bbox[0][1] <<
"]\n" 3923 <<
"[ymin,ymax] = [" << bbox[1][0] <<
',' << bbox[1][1] <<
"]\n" 3924 <<
"[zmin,zmax] = [" << bbox[2][0] <<
',' << bbox[2][1] <<
']' 3927 os <<
"endsolid GridFunction" << endl;
3944 MFEM_ASSERT(new_vertex.Size() == mesh->
GetNV(),
"");
3947 old_vertex.SetSize(new_vertex.Size());
3948 for (
int i = 0; i < new_vertex.Size(); i++)
3950 old_vertex[new_vertex[i]] = i;
3957 for (
int i = 0; i < mesh->
GetNV(); i++)
3962 for (
int j = 0; j < new_vdofs.
Size(); j++)
3964 tmp(new_vdofs[j]) = (*this)(old_vdofs[j]);
3970 for (
int i = 0; i < mesh->
GetNEdges(); i++)
3973 if (old_vertex[ev[0]] > old_vertex[ev[1]])
3978 for (
int k = 0; k < dofs.
Size(); k++)
3980 int new_dof = dofs[k];
3981 int old_dof = dofs[(ind[k] < 0) ? -1-ind[k] : ind[k]];
3988 double sign = (ind[k] < 0) ? -1.0 : 1.0;
3989 tmp(new_vdof) = sign * (*this)(old_vdof);
4002 int with_subdomains,
4010 int nfe = ufes->
GetNE();
4014 Vector ul, fl, fla, d_xyz;
4024 if (with_subdomains)
4029 double total_error = 0.0;
4030 for (
int s = 1;
s <= nsd;
s++)
4033 u.ComputeFlux(blfi, flux, with_coeff, (with_subdomains ?
s : -1));
4035 for (
int i = 0; i < nfe; i++)
4037 if (with_subdomains && ufes->
GetAttribute(i) !=
s) {
continue; }
4042 u.GetSubVector(udofs, ul);
4055 *ffes->
GetFE(i), fl, with_coeff);
4060 (aniso_flags ? &d_xyz : NULL));
4062 error_estimates(i) = std::sqrt(eng);
4068 for (
int k = 0; k <
dim; k++)
4073 double thresh = 0.15 * 3.0/
dim;
4075 for (
int k = 0; k <
dim; k++)
4077 if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
4080 (*aniso_flags)[i] = flag;
4088 auto process_local_error = total_error;
4089 MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
4090 MPI_SUM, pfes->GetComm());
4092 #endif // MFEM_USE_MPI 4093 return std::sqrt(total_error);
4105 MFEM_VERIFY(
dim >= 1,
"dim must be positive");
4106 MFEM_VERIFY(
dim <= 3,
"dim cannot be greater than 3");
4107 MFEM_VERIFY(order >= 0,
"order cannot be negative");
4109 bool rotate = (angle != 0.0) || (midpoint->
Norml2() != 0.0);
4118 x[0] = tmp[0]*cos(-angle) - tmp[1]*sin(-angle);
4119 x[1] = tmp[0]*sin(-angle) + tmp[1]*cos(-angle);
4128 double x1 = (x(0) - xmin(0))/(xmax(0)-xmin(0)), x2, x3;
4129 Vector poly_x(order+1), poly_y(order+1), poly_z(order+1);
4133 x2 = (x(1)-xmin(1))/(xmax(1)-xmin(1));
4138 x3 = (x(2)-xmin(2))/(xmax(2)-xmin(2));
4142 int basis_dimension =
static_cast<int>(pow(order+1,
dim));
4143 poly.
SetSize(basis_dimension);
4148 for (
int i = 0; i <= order; i++)
4150 poly(i) = poly_x(i);
4156 for (
int j = 0; j <= order; j++)
4158 for (
int i = 0; i <= order; i++)
4160 int cnt = i + (order+1) * j;
4161 poly(cnt) = poly_x(i) * poly_y(j);
4168 for (
int k = 0; k <= order; k++)
4170 for (
int j = 0; j <= order; j++)
4172 for (
int i = 0; i <= order; i++)
4174 int cnt = i + (order+1) * j + (order+1) * (order+1) * k;
4175 poly(cnt) = poly_x(i) * poly_y(j) * poly_z(k);
4183 MFEM_ABORT(
"TensorProductLegendre: invalid value of dim");
4199 int num_elems = patch.
Size();
4209 if (
rotate && iface >= 0)
4215 physical_diff = 0.0;
4218 for (
int i = 0; i < 2; i++)
4220 reference_pt.
Set1w((
double)i, 0.0);
4221 Tr.
Transform(reference_pt, physical_pt);
4222 midpoint += physical_pt;
4223 physical_pt *= pow(-1.0,i);
4224 physical_diff += physical_pt;
4227 angle = atan2(physical_diff(1),physical_diff(0));
4230 for (
int i = 0; i < num_elems; i++)
4232 int ielem = patch[i];
4243 transip -= midpoint;
4246 transip[0] = tmp[0]*cos(-angle) - tmp[1]*sin(-angle);
4247 transip[1] = tmp[0]*sin(-angle) + tmp[1]*cos(-angle);
4249 for (
int d = 0; d <
dim; d++) { xmax(d) = max(xmax(d), transip(d)); }
4250 for (
int d = 0; d <
dim; d++) { xmin(d) = min(xmin(d), transip(d)); }
4258 bool subdomain_reconstruction,
4260 double tichonov_coeff)
4262 MFEM_VERIFY(tichonov_coeff >= 0.0,
"tichonov_coeff cannot be negative");
4269 int nfe = ufes->
GetNE();
4270 int nfaces = ufes->
GetNF();
4277 error_estimates = 0.0;
4288 if (subdomain_reconstruction)
4293 double total_error = 0.0;
4294 for (
int iface = 0; iface < nfaces; iface++)
4301 patch[0] = el1; patch[1] = el2;
4304 if (el1 == -1 || el2 == -1)
4315 if (el1_attr != el2_attr) {
continue; }
4324 int num_basis_functions =
static_cast<int>(pow(patch_order+1,
dim));
4325 int flux_order = 2*patch_order + 1;
4334 xmin, xmax, angle, midpoint, iface);
4339 for (
int i = 0; i < patch.
Size(); i++)
4341 int ielem = patch[i];
4347 u.GetSubVector(udofs, ul);
4355 *dummy, fl, with_coeff, ir);
4359 for (
int k = 0; k < num_integration_pts; k++)
4371 for (
int l = 0; l < num_basis_functions; l++)
4374 for (
int n = 0; n < sdim; n++)
4376 b[l + n * num_basis_functions] +=
p(l) * fl(k + n * num_integration_pts);
4388 for (
int i = 0; i < num_basis_functions; i++)
4390 A(i,i) += tichonov_coeff;
4397 if (!lu.Factor(num_basis_functions,TOL))
4400 mfem::out <<
"LSZZErrorEstimator: Matrix A is singular.\t" 4401 <<
"Consider increasing tichonov_coeff." << endl;
4402 for (
int i = 0; i < num_basis_functions; i++)
4406 lu.Factor(num_basis_functions,TOL);
4408 lu.Solve(num_basis_functions, sdim,
b);
4416 for (
int i = 0; i < num_basis_functions; i++)
4418 for (
int j = 0; j < sdim; j++)
4420 f(j) +=
b[i + j * num_basis_functions] *
p(i);
4427 double element_error = 0.0;
4428 double patch_error = 0.0;
4429 for (
int i = 0; i < patch.
Size(); i++)
4431 int ielem = patch[i];
4432 element_error =
u.ComputeElementGradError(ielem, &global_poly);
4433 element_error *= element_error;
4434 patch_error += element_error;
4435 error_estimates(ielem) += element_error;
4439 total_error += patch_error;
4447 for (
int ielem = 0; ielem < nfe; ielem++)
4449 if (counters[ielem] == 0)
4451 error_estimates(ielem) =
infinity();
4455 error_estimates(ielem) /= counters[ielem]/2.0;
4456 error_estimates(ielem) = sqrt(error_estimates(ielem));
4459 return std::sqrt(total_error/
dim);
4481 for (
int j = 0; j < nip; j++)
4484 T->SetIntPoint(&ip);
4490 double errj = val1.
Norml2();
4493 errj = pow(errj,
p);
4494 norm += ip.
weight * T->Weight() * errj;
4498 norm = std::max(norm, errj);
4507 norm = -pow(-norm, 1./
p);
4511 norm = pow(norm, 1./
p);
4525 return sol_in.
Eval(*T_in, ip);
4536 string cname = name;
4537 if (cname ==
"Linear")
4541 else if (cname ==
"Quadratic")
4545 else if (cname ==
"Cubic")
4549 else if (!strncmp(name,
"H1_", 3))
4553 else if (!strncmp(name,
"H1Pos_", 6))
4558 else if (!strncmp(name,
"L2_T", 4))
4562 else if (!strncmp(name,
"L2_", 3))
4568 mfem::err <<
"Extrude1DGridFunction : unknown FE collection : " void SaveSTLTri(std::ostream &out, double p1[], double p2[], double p3[])
Abstract class for all finite elements.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
virtual double ComputeGradError(VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 for H1 or L2 elements.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
virtual void CountElementsPerVDof(Array< int > &elem_per_vdof) const
For each vdof, counts how many elements contain the vdof, as containment is determined by FiniteEleme...
bool IsLegacyLoaded() const
I/O: Return true if the mesh was loaded from the legacy v1.1 format.
int GetNPoints() const
Returns the number of the points in the integration rule.
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 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.
Class used for extruding scalar GridFunctions.
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
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)
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Base class for vector Coefficients that optionally depend on time and space.
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Geometry::Type GetElementBaseGeometry(int i) 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 ...
bool Nonconforming() const
int Dimension() const
Dimension of the reference space used within the elements.
virtual int GetContType() const =0
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void SetSize(int s)
Resize the vector to size s.
void RestrictConforming()
GridFunction * Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d, GridFunction *sol, const int ny)
Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Compute the vector gradient with respect to the physical element variable.
void LegacyToNewVertexOrdering(Array< int > &order) const
I/O: Return a map from old (v1.1) vertex indices to new vertex indices.
double Scale()
Return the scale factor times the optional time dependent function. Returns with when not set by th...
void GetVertexVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified vertices.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
virtual void ComputeElementLpErrors(const double p, Coefficient &exsol, Vector &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
void GetGradient(ElementTransformation &tr, Vector &grad) const
Gradient of a scalar function at a quadrature point.
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
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...
bool Nonconforming() const
int Size() const
Returns the size of the vector.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void GetElementAverages(GridFunction &avgs) const
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 ...
const NURBSExtension * GetNURBSext() const
double * Data() const
Returns the matrix data array.
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
double Tol()
Return the tolerance used to identify the mesh vertices.
int GetNEdges() const
Return the number of edges.
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.
void LoadSolution(std::istream &input, GridFunction &sol) const
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...
double Eval(double h, int p) const
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
static void CalcLegendre(const int p, const double x, double *u)
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
void CalcOrtho(const DenseMatrix &J, Vector &n)
void GetHessians(int i, const IntegrationRule &ir, DenseMatrix &hess, int vdim=1) const
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr)
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
void ImposeBounds(int i, const Vector &weights, const Vector &lo_, const Vector &hi_)
std::function< double(const Vector &)> f(double mass_coeff)
Vector & operator=(const double *v)
Copy Size() entries from v.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
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 GetNEDofs() const
Number of all scalar edge-interior dofs.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
virtual double ComputeMaxError(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void GetDerivative(int comp, int der_comp, GridFunction &der)
Compute a certain derivative of a function's component. Derivatives of the function are computed at t...
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
void GetCurl(ElementTransformation &tr, Vector &curl) const
GeometryRefiner GlobGeometryRefiner
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
double GetElementSize(ElementTransformation *T, int type=0)
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
void Swap(Vector &other)
Swap the contents of two Vectors.
const FiniteElementCollection * FEColl() const
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
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 MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
void SetMaxIter(int max_it)
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
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 BoundingBox(const Array< int > &patch, FiniteElementSpace *ufes, int order, Vector &xmin, Vector &xmax, double &angle, Vector &midpoint, int iface)
Defines the bounding box for the face patches used by NewZZErorrEstimator.
virtual double ComputeHDivError(VectorCoefficient *exsol, Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(div)-norm for RT elements.
void SetLinearConstraint(const Vector &w_, double a_)
double Sum() const
Return the sum of the vector entries.
virtual double ComputeCurlError(VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns ||curl u_ex - curl u_h||_L2 for ND elements.
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
virtual const char * Name() const
double LSZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, Vector &error_estimates, bool subdomain_reconstruction, bool with_coeff, double tichonov_coeff)
A ‘‘true’’ ZZ error estimator that uses face-based patches for flux reconstruction.
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
void AccumulateAndCountDerivativeValues(int comp, int der_comp, GridFunction &der, Array< int > &zones_per_dof)
Used for the serial and parallel implementations of the GetDerivative() method; see its documentation...
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
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.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
double GetDivergence(ElementTransformation &tr) const
int GetVDim()
Returns dimension of the vector.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
FiniteElementSpace * FESpace()
GridFunction & operator=(const GridFunction &rhs)
Copy assignment. Only the data of the base class Vector is copied.
void Mult(const double *x, double *y) const
Matrix vector multiplication.
int GetNE() const
Returns number of elements in the mesh.
virtual double ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, const Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
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...
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
double * GetData() const
Return a pointer to the beginning of the Vector data.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
int GetVDim() const
Returns vector dimension.
A general vector function coefficient.
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Geometry::Type GetElementGeometry(int i) const
virtual void GetElementDofValues(int el, Vector &dof_vals) const
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
bool IsIdentityProlongation(const Operator *P)
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
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]...
Mesh * GetMesh() const
Returns the mesh.
void SetAbsTol(double atol)
double p(const Vector &x, double t)
void SetRelTol(double rtol)
int GetDim() const
Returns the reference space dimension for the finite element.
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
double Distance(const double *x, const double *y, const int n)
virtual double ComputeElementGradError(int ielem, VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 in element ielem for H1 or L2 elements.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
virtual double ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, class JumpScaling jump_scaling, const IntegrationRule *irs[]=NULL) const
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...
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
virtual double ComputeDivError(Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns ||div u_ex - div u_h||_L2 for RT elements.
double Min() const
Returns the minimal element of the vector.
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int GetElementForDof(int i) const
Return the index of the first element that contains dof i.
void GetVectorFieldNodalValues(Vector &val, int comp) const
int GetEdgeDofs(int edge, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
int NumBdr(int GeomType)
Return the number of boundary "faces" of a given Geometry::Type.
void GetColumnReference(int c, Vector &col)
FiniteElementCollection * fec
Used when the grid function is read from a file. It can also be set explicitly, see MakeOwner()...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
Compute the vector gradient with respect to the reference element variable.
void TensorProductLegendre(int dim, int order, const Vector &x_in, const Vector &xmax, const Vector &xmin, Vector &poly, double angle, const Vector *midpoint)
Defines the global tensor product polynomial space used by NewZZErorrEstimator.
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
double Max() const
Returns the maximal element of the vector.
int GetNVDofs() const
Number of all scalar vertex dofs.
double Norml1() const
Returns the l_1 norm of the vector.
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.
int SpaceDimension() const
Dimension of the physical space containing the mesh.
int GetAttribute(int i) 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.
int GetNE() const
Returns number of elements.
virtual double ComputeLpError(const double p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) 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.
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
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 ...
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Class for integration point with weight.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
virtual void Mult(const Vector &xt, Vector &x) const
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)
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
int GetBdrAttribute(int i) const
Ordering::Type GetOrdering() const
Return the ordering method.
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
void Save(std::ostream &out) const
Save finite element space to output stream out.
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
void GetVectorFieldValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr, int comp=0) const
void NewMemoryAndSize(const Memory< double > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Extension of GetGradient(...) for a collection of IntegrationPoints.
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
int GetNBE() const
Returns number of boundary elements in the mesh.
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...
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
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.
int GetFaceValues(int i, int side, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, int vdim=1) const
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
int GetFaceVectorValues(int i, int side, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
double Norml2() const
Returns the l2 norm of the vector.
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 nonconforming mesh extension.
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps, int vdim=1) const
int Size() const
Return the logical size of the array.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
int GetNV() const
Returns number of vertices in the mesh.
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)
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Arbitrary order H1-conforming (continuous) finite elements.
void pts(int iphi, int t, double x[])
Field is continuous across element interfaces.
double u(const Vector &xvec)
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
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 Save(std::ostream &out) const
Save the GridFunction to an output stream.
void SetBounds(const Vector &lo_, const Vector &hi_)
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
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...
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...
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
void Set1w(const double x1, const double w)
int GetNFDofs() const
Number of all scalar face-interior dofs.
void GetValuesFrom(const GridFunction &orig_func)
int GetLocalDofForDof(int i) const
Return the local dof index in the first element that contains dof i.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
void GetBdrValuesFrom(const GridFunction &orig_func)
static void AdjustVDofs(Array< int > &vdofs)
Remove the orientation information encoded into an array of dofs Some basis function types have a rel...
Array< int > attributes
A list of all unique element attributes used by the Mesh.
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 AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr, Array< int > &values_counter)
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...
Arbitrary order "L2-conforming" discontinuous finite elements.
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.