15 #include "../mesh/nurbs.hpp" 16 #include "../general/text.hpp" 37 istream::int_type next_char = input.peek();
43 if (buff ==
"NURBS_patches")
46 "NURBS_patches requires NURBS FE space");
51 MFEM_ABORT(
"unknown section: " << buff);
83 int g_nddofs = g_ndofs - (g_nvdofs + g_nedofs + g_nfdofs);
85 vi = ei = fi = di = 0;
86 for (
int i = 0; i < num_pieces; i++)
93 int l_nddofs = l_ndofs - (l_nvdofs + l_nedofs + l_nfdofs);
94 const double *l_data = gf_array[i]->
GetData();
95 double *g_data =
data;
98 for (
int d = 0; d < vdim; d++)
100 memcpy(g_data+vi, l_data, l_nvdofs*
sizeof(
double));
103 memcpy(g_data+ei, l_data, l_nedofs*
sizeof(
double));
106 memcpy(g_data+fi, l_data, l_nfdofs*
sizeof(
double));
109 memcpy(g_data+di, l_data, l_nddofs*
sizeof(
double));
116 memcpy(g_data+vdim*vi, l_data, vdim*l_nvdofs*
sizeof(
double));
117 l_data += vdim*l_nvdofs;
118 g_data += vdim*g_nvdofs;
119 memcpy(g_data+vdim*ei, l_data, vdim*l_nedofs*
sizeof(
double));
120 l_data += vdim*l_nedofs;
121 g_data += vdim*g_nedofs;
122 memcpy(g_data+vdim*fi, l_data, vdim*l_nfdofs*
sizeof(
double));
123 l_data += vdim*l_nfdofs;
124 g_data += vdim*g_nfdofs;
125 memcpy(g_data+vdim*di, l_data, vdim*l_nddofs*
sizeof(
double));
126 l_data += vdim*l_nddofs;
127 g_data += vdim*g_nddofs;
155 MFEM_ABORT(
"Error in update sequence. GridFunction needs to be updated " 156 "right after the space is updated.");
164 old_data.
Swap(*
this);
166 T->
Mult(old_data, *
this);
242 int nfe = ufes->
GetNE();
250 for (
int i = 0; i < nfe; i++)
252 if (subdomain >= 0 && ufes->
GetAttribute(i) != subdomain)
264 *ffes->
GetFE(i), fl, wcoef);
269 for (
int j = 0; j < fdofs.
Size(); j++)
285 for (
int i = 0; i < count.Size(); i++)
287 if (count[i] != 0) { flux(i) /= count[i]; }
297 static const int geoms[3] =
354 int dof = FElem->
GetDof();
364 "invalid FE map type");
366 for (k = 0; k < n; k++)
369 nval[k] = shape * ((
const double *)loc_data + dof * vdim);
376 for (k = 0; k < n; k++)
380 nval[k] = loc_data * (&vshape(0,vdim));
397 return (DofVal * LocVec);
404 int dof = FElem->
GetDof();
412 "invalid FE map type");
417 for (
int k = 0; k < vdim; k++)
419 val(k) = shape * ((
const double *)loc_data + dof * k);
445 "invalid FE map type");
446 int dof = FElem->
GetDof();
447 Vector DofVal(dof), loc_data(dof);
449 for (
int k = 0; k < n; k++)
452 vals(k) = DofVal * loc_data;
522 int dof = FElem->
GetDof();
531 "invalid FE map type");
535 for (
int j = 0; j < nip; j++)
539 for (
int k = 0; k < vdim; k++)
541 vals(k,j) = shape * ((
const double *)loc_data + dof * k);
551 for (
int j = 0; j < nip; j++)
620 Vector shape, loc_values, orig_loc_values;
621 int i, j, d, ne, dof, odof, vdim;
625 for (i = 0; i < ne; i++)
634 loc_values.
SetSize(dof * vdim);
637 for (j = 0; j < dof; j++)
641 for (d = 0; d < vdim; d++)
643 loc_values(d*dof+j) =
644 shape * ((
const double *)orig_loc_values + d * odof) ;
657 Vector shape, loc_values, orig_loc_values;
658 int i, j, d, nbe, dof, odof, vdim;
662 for (i = 0; i < nbe; i++)
671 loc_values.
SetSize(dof * vdim);
674 for (j = 0; j < dof; j++)
678 for (d = 0; d < vdim; d++)
680 loc_values(d*dof+j) =
681 shape * ((
const double *)orig_loc_values + d * odof);
695 int d, j, k, n, sdim, dof, ind;
702 int *dofs = &vdofs[comp*dof];
708 for (k = 0; k < n; k++)
713 for (d = 0; d < sdim; d++)
716 for (j = 0; j < dof; j++)
717 if ( (ind=dofs[j]) >= 0 )
719 a += vshape(j, d) *
data[ind];
723 a -= vshape(j, d) *
data[-1-ind];
740 double *temp =
new double[
size];
743 for (j = 0; j < ndofs; j++)
744 for (i = 0; i < vdim; i++)
746 temp[j+i*ndofs] =
data[k++];
749 for (i = 0; i <
size; i++)
777 val(vertices[k]) += vals(k, comp);
778 overlap[vertices[k]]++;
782 for (i = 0; i < overlap.Size(); i++)
784 val(i) /= overlap[i];
792 int d, i, k, ind, dof, sdim;
801 for (i = 0; i < new_fes->
GetNE(); i++)
808 for (d = 0; d < sdim; d++)
810 for (k = 0; k < dof; k++)
812 if ( (ind=new_vdofs[dof*d+k]) < 0 )
814 ind = -1-ind, vals(k, d) = - vals(k, d);
816 vec_field(ind) += vals(k, d);
822 for (i = 0; i < overlap.Size(); i++)
824 vec_field(i) /= overlap[i];
836 int i, j, k,
dim, dof, der_dof, ind;
839 for (i = 0; i < overlap.Size(); i++)
846 for (i = 0; i < der_fes->
GetNE(); i++)
855 der_dof = der_fe->
GetDof();
861 for (j = 0; j < dof; j++)
862 loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
864 for (k = 0; k < der_dof; k++)
872 for (j = 0; j <
dim; j++)
874 a += inv_jac(j, der_comp) * pt_grad(j);
876 der(der_dofs[k]) += a;
877 overlap[der_dofs[k]]++;
881 for (i = 0; i < overlap.Size(); i++)
883 der(i) /= overlap[i];
904 MultAtB(loc_data_mat, dshape, gh);
915 "invalid FE map type");
922 for (
int i = 0; i < Jinv.Width(); i++)
924 for (
int j = 0; j < Jinv.Height(); j++)
926 div_v += grad_hat(i, j) * Jinv(j, i);
938 div_v = (loc_data * divshape) / tr.
Weight();
950 "invalid FE map type");
957 Mult(grad_hat, Jinv, grad);
958 MFEM_ASSERT(grad.Height() == grad.Width(),
"");
959 if (grad.Height() == 3)
962 curl(0) = grad(2,1) - grad(1,2);
963 curl(1) = grad(0,2) - grad(2,0);
964 curl(2) = grad(1,0) - grad(0,1);
966 else if (grad.Height() == 2)
969 curl(0) = grad(1,0) - grad(0,1);
981 curl.
SetSize(curl_shape.Width());
982 if (curl_shape.Width() == 3)
985 curl_shape.MultTranspose(loc_data, curl_hat);
990 curl_shape.MultTranspose(loc_data, curl);
1010 dshape.MultTranspose(lval, gh);
1032 dshape.MultTranspose(lval, gh);
1033 Tr->SetIntPoint(&ip);
1036 Jinv.MultTranspose(gh, gcol);
1044 "invalid FE map type");
1050 grad.
SetSize(grad_hat.Height(), Jinv.Width());
1051 Mult(grad_hat, Jinv, grad);
1059 Vector loc_avgs, loc_this;
1064 for (
int i = 0; i <
fes->
GetNE(); i++)
1069 avgs.FESpace()->GetElementDofs(i, te_dofs);
1072 loc_mass.
Mult(loc_this, loc_avgs);
1073 avgs.AddElementVector(te_dofs, loc_avgs);
1075 loc_mass.
Mult(loc_this, loc_avgs);
1076 int_psi.AddElementVector(te_dofs, loc_avgs);
1078 for (
int i = 0; i < avgs.Size(); i++)
1080 avgs(i) /= int_psi(i);
1099 mfem_error(
"GridFunction::ProjectGridFunction() :" 1100 " incompatible vector dimensions!");
1104 for (
int i = 0; i < mesh->
GetNE(); i++)
1108 for (
int vd = 0; vd < vdim; vd++)
1126 MFEM_ASSERT(weights.
Size() ==
size,
"Different # of weights and dofs.");
1127 MFEM_ASSERT(_lo.
Size() ==
size,
"Different # of lower bounds and dofs.");
1128 MFEM_ASSERT(_hi.
Size() ==
size,
"Different # of upper bounds and dofs.");
1131 double tol = 1.e-12;
1139 slbqp.
Mult(vals, new_vals);
1145 double _min,
double _max)
1153 double max_val = vals.
Max();
1154 double min_val = vals.
Min();
1156 if (max_val <= _min)
1163 if (_min <= min_val && max_val <= _max)
1169 minv = (_min > min_val) ? _min : min_val;
1170 maxv = (_max < max_val) ? _max : max_val;
1189 for (j = 0; j < vertices.
Size(); j++)
1191 nval(vertices[j]) += values[j];
1192 overlap[vertices[j]]++;
1195 for (i = 0; i < overlap.Size(); i++)
1197 nval(i) /= overlap[i];
1212 for (
int i = 0; i <
fes->
GetNE(); i++)
1220 for (
int j = 0; j < vdofs.
Size(); j++)
1224 MFEM_VERIFY(vals[j] != 0.0,
1225 "Coefficient has zeros, harmonic avg is undefined!");
1226 (*this)(vdofs[j]) += 1.0 / vals[j];
1230 (*this)(vdofs[j]) += vals[j];
1232 else { MFEM_ABORT(
"Not implemented"); }
1234 zones_per_vdof[vdofs[j]]++;
1250 for (
int i = 0; i <
fes->
GetNE(); i++)
1258 for (
int j = 0; j < vdofs.
Size(); j++)
1275 MFEM_VERIFY(vals[j] != 0.0,
1276 "Coefficient has zeros, harmonic avg is undefined!");
1277 (*this)(ldof) += isign / vals[j];
1281 (*this)(ldof) += isign*vals[j];
1284 else { MFEM_ABORT(
"Not implemented"); }
1286 zones_per_vdof[ldof]++;
1296 for (
int i = 0; i <
size; i++)
1298 (*this)(i) /= zones_per_vdof[i];
1303 for (
int i = 0; i <
size; i++)
1305 (*this)(i) = zones_per_vdof[i]/(*
this)(i);
1310 MFEM_ABORT(
"invalud AvgType");
1325 const double *center = delta_coeff.
Center();
1326 const double *vert = mesh->
GetVertex(0);
1327 double min_dist, dist;
1332 for (
int i = 0; i < mesh->
GetNV(); i++)
1336 if (dist < min_dist)
1346 if (min_dist >= delta_coeff.
Tol())
1355 Vector vals, loc_mass_vals;
1356 for (
int i = 0; i < mesh->
GetNE(); i++)
1359 for (
int j = 0; j < vertices.
Size(); j++)
1360 if (vertices[j] == v_idx)
1370 loc_mass.Mult(vals, loc_mass_vals);
1371 integral += loc_mass_vals.
Sum();
1381 if (delta_c == NULL)
1386 for (
int i = 0; i <
fes->
GetNE(); i++)
1400 (*this) *= (delta_c->
Scale() / integral);
1411 for (
int i = 0; i < dofs.
Size(); i++)
1424 (*this)(vdof) = coeff.
Eval(*T, ip);
1452 for (
int i = 0; i < dofs.
Size(); i++)
1464 vcoeff.
Eval(val, *T, ip);
1468 (*this)(vdof) = val(vd);
1475 int i, j, fdof, d, ind, vdim;
1489 for (j = 0; j < fdof; j++)
1493 for (d = 0; d < vdim; d++)
1495 if (!coeff[d]) {
continue; }
1497 val = coeff[d]->
Eval(*transf, ip);
1498 if ( (ind = vdofs[fdof*d+j]) < 0 )
1500 val = -val, ind = -1-ind;
1519 for (
int i = 0; i <
fes->
GetNE(); i++)
1528 for (
int j = 0; j < vdofs.
Size(); j++)
1530 if (attr > dof_attr[vdofs[j]])
1532 (*this)(vdofs[j]) = vals[j];
1533 dof_attr[vdofs[j]] = attr;
1568 int i, j, fdof, d, ind, vdim;
1585 for (j = 0; j < fdof; j++)
1589 for (d = 0; d < vdim; d++)
1591 if (!coeff[d]) {
continue; }
1593 val = coeff[d]->
Eval(*transf, ip);
1594 if ( (ind = vdofs[fdof*d+j]) < 0 )
1596 val = -val, ind = -1-ind;
1622 for (i = 0; i < bdr_edges.
Size(); i++)
1624 int edge = bdr_edges[i];
1626 if (vdofs.
Size() == 0) {
continue; }
1632 for (d = 0; d < vdim; d++)
1634 if (!coeff[d]) {
continue; }
1636 fe->
Project(*coeff[d], *transf, vals);
1637 for (
int k = 0; k < vals.
Size(); k++)
1639 (*this)(vdofs[d*vals.
Size()+k]) = vals(k);
1676 vcoeff.
Eval(vc, *T, ip);
1679 lvec.
Add(ip.
weight * (vc * nor), shape);
1707 vcoeff.
Eval(vc, *T, ip);
1709 lvec(j) = (vc * nor);
1735 fe->
Project(vcoeff, *T, lvec);
1746 for (
int i = 0; i < bdr_edges.
Size(); i++)
1748 int edge = bdr_edges[i];
1750 if (dofs.
Size() == 0) {
continue; }
1756 fe->
Project(vcoeff, *T, lvec);
1765 double error = 0.0, a;
1770 int fdof, d, i, intorder, j, k;
1796 for (k = 0; k < fdof; k++)
1797 if (vdofs[fdof*d+k] >= 0)
1799 a += (*this)(vdofs[fdof*d+k]) * shape(k);
1803 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
1806 a -= exsol[d]->
Eval(*transf, ip);
1814 return -sqrt(-error);
1829 for (
int i = 0; i <
fes->
GetNE(); i++)
1831 if (elems != NULL && (*elems)[i] == 0) {
continue; }
1833 int intorder = 2*fe->
GetOrder() + 1;
1845 exsol.
Eval(exact_vals, *T, *ir);
1848 vals.
Norm2(loc_errs);
1853 error += ip.
weight * T->
Weight() * (loc_errs(j) * loc_errs(j));
1859 return -sqrt(-error);
1866 Coefficient *ell_coeff,
double Nu,
int norm_type)
const 1869 int i, fdof,
dim, intorder, j, k;
1874 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
1887 for (i = 0; i < mesh->
GetNE(); i++)
1898 for (k = 0; k < fdof; k++)
1901 el_dofs(k) = (*this)(vdofs[k]);
1905 el_dofs(k) = - (*this)(-1-vdofs[k]);
1912 exgrad->
Eval(e_grad, *transf, ip);
1914 Mult(dshape, Jinv, dshapet);
1918 ell_coeff->
Eval(*transf, ip) *
1927 int i1 = face_elem_transf->
Elem1No;
1928 int i2 = face_elem_transf->
Elem2No;
1935 intorder = 2 * intorder;
1941 transf = face_elem_transf->
Elem1;
1947 for (k = 0; k < fdof; k++)
1950 el_dofs(k) = (*this)(vdofs[k]);
1954 el_dofs(k) = - (*this)(-1-vdofs[k]);
1961 ell_coeff_val(j) = ell_coeff->
Eval(*transf, eip);
1962 err_val(j) = exsol->
Eval(*transf, eip) - (shape * el_dofs);
1968 transf = face_elem_transf->
Elem2;
1974 for (k = 0; k < fdof; k++)
1977 el_dofs(k) = (*this)(vdofs[k]);
1981 el_dofs(k) = - (*this)(-1-vdofs[k]);
1988 ell_coeff_val(j) += ell_coeff->
Eval(*transf, eip);
1989 ell_coeff_val(j) *= 0.5;
1990 err_val(j) -= (exsol->
Eval(*transf, eip) - (shape * el_dofs));
1994 transf = face_elem_transf->
Face;
1999 error += (ip.
weight * Nu * ell_coeff_val(j) *
2000 pow(transf->
Weight(), 1.0-1.0/(
dim-1)) *
2001 err_val(j) * err_val(j));
2007 return -sqrt(-error);
2015 double error = 0.0, a;
2020 int fdof, d, i, intorder, j, k;
2047 for (k = 0; k < fdof; k++)
2048 if (vdofs[fdof*d+k] >= 0)
2050 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2054 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2056 a -= exsol[d]->
Eval(*transf, ip);
2074 int i, fdof,
dim, intorder, j, k;
2078 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
2081 double a, error = 0.0;
2090 for (i = 0; i < mesh->
GetNE(); i++)
2092 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2109 for (k = 0; k < fdof; k++)
2112 el_dofs(k) = (*this)(vdofs[k]);
2116 el_dofs(k) = -(*this)(-1-vdofs[k]);
2123 a = (el_dofs * shape) - (exsol->
Eval(*transf, ip));
2129 for (i = 0; i < mesh->
GetNE(); i++)
2131 if (elems != NULL && (*elems)[i] == 0) {
continue; }
2149 for (k = 0; k < fdof; k++)
2152 el_dofs(k) = (*this)(vdofs[k]);
2156 el_dofs(k) = -(*this)(-1-vdofs[k]);
2163 exgrad->
Eval(e_grad, *transf, ip);
2165 Mult(dshape, Jinv, dshapet);
2184 for (
int i = 0; i <
fes->
GetNE(); i++)
2194 int intorder = 2*fe->
GetOrder() + 1;
2203 double err = fabs(vals(j) - exsol.
Eval(*T, ip));
2219 error = std::max(error,
err);
2229 error = -pow(-error, 1./p);
2233 error = pow(error, 1./p);
2251 for (
int i = 0; i <
fes->
GetNE(); i++)
2261 int intorder = 2*fe->
GetOrder() + 1;
2266 exsol.
Eval(exact_vals, *T, *ir);
2273 vals.
Norm2(loc_errs);
2277 v_weight->
Eval(exact_vals, *T, *ir);
2280 for (
int j = 0; j < vals.
Width(); j++)
2283 for (
int d = 0; d < vals.
Height(); d++)
2285 err += vals(d,j)*exact_vals(d,j);
2287 loc_errs(j) = fabs(
err);
2294 double err = loc_errs(j);
2310 error = std::max(error,
err);
2320 error = -pow(-error, 1./p);
2324 error = pow(error, 1./p);
2357 out <<
"NURBS_patches\n";
2386 out <<
"SCALARS " << field_name <<
" double 1\n" 2387 <<
"LOOKUP_TABLE default\n";
2388 for (
int i = 0; i < mesh->
GetNE(); i++)
2395 for (
int j = 0; j < val.
Size(); j++)
2397 out << val(j) <<
'\n';
2401 else if ( (vec_dim == 2 || vec_dim == 3) && mesh->
SpaceDimension() > 1)
2404 out <<
"VECTORS " << field_name <<
" double\n";
2405 for (
int i = 0; i < mesh->
GetNE(); i++)
2412 for (
int j = 0; j < vval.
Width(); j++)
2414 out << vval(0, j) <<
' ' << vval(1, j) <<
' ';
2430 for (
int vd = 0; vd < vec_dim; vd++)
2432 out <<
"SCALARS " << field_name << vd <<
" double 1\n" 2433 <<
"LOOKUP_TABLE default\n";
2434 for (
int i = 0; i < mesh->
GetNE(); i++)
2441 for (
int j = 0; j < val.
Size(); j++)
2443 out << val(j) <<
'\n';
2454 double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
2455 double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
2456 double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
2457 v1[2] * v2[0] - v1[0] * v2[2],
2458 v1[0] * v2[1] - v1[1] * v2[0]
2460 double rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2461 n[0] *= rl; n[1] *= rl; n[2] *= rl;
2463 out <<
" facet normal " << n[0] <<
' ' << n[1] <<
' ' << n[2]
2465 <<
"\n vertex " << p1[0] <<
' ' << p1[1] <<
' ' << p1[2]
2466 <<
"\n vertex " << p2[0] <<
' ' << p2[1] <<
' ' << p2[2]
2467 <<
"\n vertex " << p3[0] <<
' ' << p3[1] <<
' ' << p3[2]
2468 <<
"\n endloop\n endfacet\n";
2484 double pts[4][3], bbox[3][2];
2486 out <<
"solid GridFunction\n";
2488 bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
2489 bbox[2][0] = bbox[2][1] = 0.0;
2490 for (i = 0; i < mesh->
GetNE(); i++)
2497 for (k = 0; k < RG.
Size()/n; k++)
2499 for (j = 0; j < n; j++)
2502 pts[j][0] = pointmat(0,l);
2503 pts[j][1] = pointmat(1,l);
2504 pts[j][2] = values(l);
2520 bbox[0][0] = pointmat(0,0);
2521 bbox[0][1] = pointmat(0,0);
2522 bbox[1][0] = pointmat(1,0);
2523 bbox[1][1] = pointmat(1,0);
2524 bbox[2][0] = values(0);
2525 bbox[2][1] = values(0);
2528 for (j = 0; j < values.
Size(); j++)
2530 if (bbox[0][0] > pointmat(0,j))
2532 bbox[0][0] = pointmat(0,j);
2534 if (bbox[0][1] < pointmat(0,j))
2536 bbox[0][1] = pointmat(0,j);
2538 if (bbox[1][0] > pointmat(1,j))
2540 bbox[1][0] = pointmat(1,j);
2542 if (bbox[1][1] < pointmat(1,j))
2544 bbox[1][1] = pointmat(1,j);
2546 if (bbox[2][0] > values(j))
2548 bbox[2][0] = values(j);
2550 if (bbox[2][1] < values(j))
2552 bbox[2][1] = values(j);
2557 mfem::out <<
"[xmin,xmax] = [" << bbox[0][0] <<
',' << bbox[0][1] <<
"]\n" 2558 <<
"[ymin,ymax] = [" << bbox[1][0] <<
',' << bbox[1][1] <<
"]\n" 2559 <<
"[zmin,zmax] = [" << bbox[2][0] <<
',' << bbox[2][1] <<
']' 2562 out <<
"endsolid GridFunction" << endl;
2574 const char *msg =
"invalid input stream";
2580 in >> ident; MFEM_VERIFY(ident ==
"VDim:", msg);
2607 out <<
"VDim: " <<
vdim <<
'\n' 2624 int with_subdomains)
2626 const int with_coeff = 0;
2632 int nfe = ufes->
GetNE();
2636 Vector ul, fl, fla, d_xyz;
2646 if (with_subdomains)
2648 for (
int i = 0; i < nfe; i++)
2651 if (attr > nsd) { nsd = attr; }
2655 double total_error = 0.0;
2656 for (
int s = 1; s <= nsd; s++)
2659 u.
ComputeFlux(blfi, flux, with_coeff, (with_subdomains ? s : -1));
2661 for (
int i = 0; i < nfe; i++)
2663 if (with_subdomains && ufes->
GetAttribute(i) != s) {
continue; }
2673 *ffes->
GetFE(i), fl, with_coeff);
2678 (aniso_flags ? &d_xyz : NULL));
2680 error_estimates(i) = std::sqrt(
err);
2686 for (
int k = 0; k <
dim; k++)
2691 double thresh = 0.15 * 3.0/
dim;
2693 for (
int k = 0; k <
dim; k++)
2695 if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
2698 (*aniso_flags)[i] = flag;
2703 return std::sqrt(total_error);
2726 for (
int j = 0; j < nip; j++)
2743 norm = std::max(norm,
err);
2752 norm = -pow(-norm, 1./p);
2756 norm = pow(norm, 1./p);
2770 return sol_in.
Eval(*T_in, ip);
2781 string cname = name;
2782 if (cname ==
"Linear")
2786 else if (cname ==
"Quadratic")
2790 else if (cname ==
"Cubic")
2794 else if (!strncmp(name,
"H1_", 3))
2798 else if (!strncmp(name,
"H1Pos_", 6))
2803 else if (!strncmp(name,
"L2_T", 4))
2807 else if (!strncmp(name,
"L2_", 3))
2813 mfem::err <<
"Extrude1DGridFunction : unknown FE collection : " void SaveSTLTri(std::ostream &out, double p1[], double p2[], double p3[])
Abstract class for Finite Elements.
std::ostream & operator<<(std::ostream &out, const Mesh &mesh)
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
For scalar fields; preserves point values.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
void GetVectorValues(ElementTransformation &T, const IntegrationRule &ir, DenseMatrix &vals) const
int GetNPoints() const
Returns the number of the points in the integration rule.
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.
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
Class for grid function - Vector with associated FE space.
Class used for extruding scalar GridFunctions.
void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
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)
virtual double ComputeLpError(const double p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
void PrintSolution(const GridFunction &sol, std::ostream &out) const
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 void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
bool Nonconforming() const
double ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains)
void SetSize(int s)
Resize the vector to size s.
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'.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
double Scale()
Return the scale set by SetScale() multiplied by the time-dependent function specified by SetFunction...
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, int wcoef=1, int subdomain=-1)
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)
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
void GetGradient(ElementTransformation &tr, Vector &grad) const
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
RefinedGeometry * Refine(int Geom, int Times, int ETimes=1)
int Size() const
Returns the size of the vector.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
void GetElementAverages(GridFunction &avgs) const
Piecewise-(bi)linear continuous finite elements.
Data type dense matrix using column-major storage.
Delta function coefficient.
const NURBSExtension * GetNURBSext() const
void ImposeBounds(int i, const Vector &weights, const Vector &_lo, const Vector &_hi)
int GetNDofs() const
Returns number of degrees of freedom.
double Tol()
See SetTol() for description of the tolerance parameter.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int vdim
Vector dimension.
int GetSize() const
Return the total number of quadrature points.
void LoadSolution(std::istream &input, GridFunction &sol) const
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
void CalcOrtho(const DenseMatrix &J, Vector &n)
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)
int GetGeomType() const
Returns the Geometry::Type of the reference element.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i'th element.
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr)
Vector & operator=(const double *v)
Copy Size() entries from v.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Piecewise-(bi)cubic continuous finite elements.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &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)
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...
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, int wcoef, int subdomain)
void SetLinearConstraint(const Vector &_w, double _a)
void GetCurl(ElementTransformation &tr, Vector &curl) const
void SetPrintLevel(int print_lvl)
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
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.
virtual double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
void GetGradients(const int elem, const IntegrationRule &ir, DenseMatrix &grad) const
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.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
virtual const SparseMatrix * GetRestrictionMatrix() const
virtual void ProjectDelta(int vertex, Vector &dofs) 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...
double Sum() const
Return the sum of the vector entries.
QuadratureFunction()
Create an empty QuadratureFunction.
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
virtual const char * Name() const
long GetSequence() const
Return update counter (see Mesh::sequence)
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
const SparseMatrix * GetConformingProlongation() const
GeometryRefiner GlobGeometryRefiner
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.
double GetDivergence(ElementTransformation &tr) const
int GetVDim()
Returns dimension of the vector.
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
FiniteElementSpace * FESpace()
void Mult(const double *x, double *y) const
Matrix vector multiplication.
int GetNE() const
Returns number of elements in the mesh.
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 (element) subvector to the vector.
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
const FiniteElement * GetEdgeElement(int i) const
QuadratureFunction & operator=(double value)
Redefine '=' for QuadratureFunction = constant.
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.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Mesh * GetMesh() const
Returns the mesh.
void SetAbsTol(double atol)
void SetRelTol(double rtol)
int GetDim() const
Returns the reference space dimension for the finite element.
double Distance(const double *x, const double *y, const int n)
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Base class Coefficient that may optionally depend on time.
FiniteElementSpace * fes
FE space on which grid function lives.
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
void Save(std::ostream &out) const
Write the QuadratureFunction to the stream out.
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
double Min() const
Returns the minimal element of the vector.
QuadratureSpace * qspace
Associated QuadratureSpace.
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
int GetElementForDof(int i) const
void GetVectorFieldNodalValues(Vector &val, int comp) const
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.
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().
int GetElementBaseGeometry(int i=0) const
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
double Max() const
Returns the maximal element of the vector.
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
int GetAttribute(int i) const
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
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)
int GetNE() const
Returns number of elements.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
virtual const Operator * GetProlongationMatrix() const
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Class for integration point with weight.
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i'th boundary element.
virtual void Mult(const Vector &xt, Vector &x) const
Operator application: y=A(x).
void SaveSTL(std::ostream &out, int TimesToRefine=1)
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
int GetBdrAttribute(int i) const
Ordering::Type GetOrdering() const
Return the ordering method.
void Save(std::ostream &out) const
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 GetEdgeDofs(int i, Array< int > &dofs) const
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th edge.
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.
void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
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().
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)
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.
NCMesh * ncmesh
Optional non-conforming mesh extension.
int Size() const
Logical size of the array.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
void filter_dos(std::string &line)
GridFunction & operator=(double value)
Redefine '=' for GridFunction = constant.
int DofToVDof(int dof, int vd, int ndofs=-1) const
int GetNV() const
Returns number of vertices in the mesh.
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
int GetNFaces() const
Return the number of faces in a 3D mesh.
Arbitrary order H1-conforming (continuous) finite elements.
Class representing the storage layout of a QuadratureFunction.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
const IntegrationRule & GetNodes() const
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Class representing a function through its values (scalar or vector) at quadrature points...
virtual double ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i'th boundary element.
void GetValuesFrom(const GridFunction &orig_func)
int GetLocalDofForDof(int i) const
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
void GetBdrValuesFrom(const GridFunction &orig_func)
static void AdjustVDofs(Array< int > &vdofs)
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 SaveVTK(std::ostream &out, const std::string &field_name, int ref)
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const =0
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.