104 bool by_vdim =
false;
105 bool snapshot =
false;
106 bool vis_mesh =
false;
109 real_t amr_threshold = 0.6;
110 const char *keys =
"gAm";
111 const char *device_config =
"cpu";
112 const char *mesh_file =
"../../data/mobius-strip.mesh";
116class Surface:
public Mesh
128 Surface(Opt &opt,
const char *file):
Mesh(file, true), opt(opt) { }
138 Surface(Opt &opt,
int nv,
int ne,
int nbe):
153 mesh =
new Mesh(*
this,
true);
155 BoundaryConditions();
161 if (opt.vis) { opt.vis = glvis.
open(
vishost, opt.visport) == 0; }
167 ByVDim(*
this, opt).Solve();
171 ByNodes(*
this, opt).Solve();
173 if (opt.vis && opt.snapshot)
176 Visualize(glvis, opt, mesh, mesh->
GetNodes());
183 if (opt.vis) { glvis.
close(); }
184 delete mesh;
delete fec;
delete fes;
187 virtual void Prefix()
192 virtual void Transform() {
if (opt.Tptr) {
Mesh::Transform(opt.Tptr); } }
194 virtual void Postfix()
199 virtual void Refine()
201 for (
int l = 0; l < opt.refine; l++)
210 for (
int i = 0; i <
nodes.Size(); i++)
219 void SnapNodesToUnitSphere()
223 for (
int i = 0; i <
nodes.FESpace()->GetNDofs(); i++)
225 for (
int d = 0; d <
SDIM; d++)
227 node(d) =
nodes(
nodes.FESpace()->DofToVDof(i, d));
229 node /= node.Norml2();
230 for (
int d = 0; d <
SDIM; d++)
232 nodes(
nodes.FESpace()->DofToVDof(i, d)) = node(d);
237 virtual void BoundaryConditions()
250 Opt &opt,
const Mesh *mesh,
251 const int w,
const int h,
255 if (opt.vis_mesh) { glvis <<
"mesh\n" << *mesh; }
256 else { glvis <<
"solution\n" << *mesh << solution; }
258 glvis <<
"window_size " << w <<
" " << h <<
"\n";
259 glvis <<
"keys " << opt.keys <<
"\n";
261 if (opt.wait) { glvis <<
"pause\n"; }
267 const Opt &opt,
const Mesh *mesh,
271 if (opt.vis_mesh) { glvis <<
"mesh\n" << *mesh; }
272 else { glvis <<
"solution\n" << *mesh << solution; }
273 if (opt.wait) { glvis <<
"pause\n"; }
274 if (opt.snapshot && opt.keys) { glvis <<
"keys " << opt.keys <<
"\n"; }
281 const char *mesh_file =
"surface.mesh";
282 const char *sol_file =
"sol.gf";
285 mfem::out <<
"Printing " << mesh_file <<
", " << sol_file << std::endl;
288 std::ofstream mesh_ofs(mesh_file);
289 mesh_ofs.precision(8);
290 mesh->
Print(mesh_ofs);
293 std::ofstream sol_ofs(sol_file);
294 sol_ofs.precision(8);
311 const int print_iter = -1, max_num_iter = 2000;
314 Solver(Surface &S, Opt &opt): opt(opt), S(S), cg(),
315 a(S.fes), x(S.fes), x0(S.fes),
b(S.fes), one(1.0)
327 constexpr bool converged =
true;
328 if (opt.pa) {
a.SetAssemblyLevel(AssemblyLevel::PARTIAL);}
329 for (
int i=0; i < opt.niters; ++i)
331 if (opt.amr) { Amr(); }
332 if (opt.vis) { Surface::Visualize(S.glvis, opt, S.mesh); }
333 if (!opt.id) {
mfem::out <<
"Iteration " << i <<
": "; }
337 if (Step() == converged) {
break; }
339 if (opt.print) { Surface::Print(opt, S.mesh, S.mesh->
GetNodes()); }
342 virtual bool Step() = 0;
345 bool Converged(
const real_t rnorm) {
return rnorm <
NRM; }
351 a.FormLinearSystem(S.bc, x, b, A, X, B);
356 a.RecoverFEMSolution(X, b, x);
357 const bool by_vdim = opt.by_vdim;
362 if (!opt.id) {
mfem::out <<
"rnorm = " << rnorm << std::endl; }
363 const real_t lambda = opt.lambda;
366 MFEM_VERIFY(!opt.radial,
"'VDim solver can't use radial option!");
367 return Converged(rnorm);
375 for (
int i = 0; i <
delta.Size()/
SDIM; i++)
379 for (
int d = 0; d <
SDIM; d++)
381 ni(d) = (*nodes)(d*ndof + i);
382 di(d) =
delta(d*ndof + i);
385 const real_t ndotd = (ni*di) / (ni*ni);
388 for (
int d = 0; d <
SDIM; d++) {
delta(d*ndof + i) = di(d); }
393 add(lambda, *
nodes, (1.0-lambda), x, x);
394 return Converged(rnorm);
399 MFEM_VERIFY(opt.amr_threshold >= 0.0 && opt.amr_threshold <= 1.0,
"");
400 Mesh *smesh = S.mesh;
402 const int NE = smesh->
GetNE();
404 for (
int e = 0; e < NE; e++)
413 for (
int q = 0; q < NQ; q++)
421 minW = std::min(minW, w);
422 maxW = std::max(maxW, w);
424 if (std::abs(maxW) != 0.0)
426 const real_t rho = minW / maxW;
427 MFEM_VERIFY(rho <= 1.0,
"");
441 S.BoundaryConditions();
447 class ByNodes:
public Solver
450 ByNodes(Surface &S, Opt &opt):
Solver(S, opt)
456 bool converge = ParAXeqB();
458 return converge ? true :
false;
463 class ByVDim:
public Solver
468 auto d_Xi = Xi.
Read();
473 d_nodes[c*ndof + i] = d_Xi[i];
479 auto d_Xi = Xi.
Write();
484 d_Xi[i] = d_nodes[c*ndof + i];
488 ByVDim(Surface &S, Opt &opt):
Solver(S, opt)
493 bool cvg[
SDIM] {
false};
494 for (
int c=0; c <
SDIM; ++c)
501 const bool converged = cvg[0] && cvg[1] && cvg[2];
502 return converged ? true :
false;
508struct MeshFromFile:
public Surface
510 MeshFromFile(Opt &opt): Surface(opt, opt.mesh_file) { }
514struct Catenoid:
public Surface
516 Catenoid(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
518 void Prefix()
override
522 for (
int i = 0; i < v2v.Size(); i++) { v2v[i] = i; }
524 for (
int j = 0; j <= opt.ny; j++)
526 const int v_old = opt.nx + j * (opt.nx + 1);
527 const int v_new = j * (opt.nx + 1);
531 for (
int i = 0; i <
GetNE(); i++)
536 for (
int j = 0; j < nv; j++)
542 for (
int i = 0; i <
GetNBE(); i++)
547 for (
int j = 0; j < nv; j++)
569struct Helicoid:
public Surface
571 Helicoid(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
578 const real_t v = 2.0*
PI*(2.0*x[1]-1.0)/3.0;
586struct Enneper:
public Surface
588 Enneper(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
594 const real_t u = 4.0*(x[0]-0.5);
595 const real_t v = 4.0*(x[1]-0.5);
596 p[0] = +
u -
u*
u*
u/3.0 +
u*v*v;
597 p[1] = -v -
u*
u*v + v*v*v/3.0;
603struct Hold:
public Surface
605 Hold(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
607 void Prefix()
override
611 for (
int i = 0; i < v2v.Size(); i++) { v2v[i] = i; }
613 for (
int j = 0; j <= opt.ny; j++)
615 const int v_old = opt.nx + j * (opt.nx + 1);
616 const int v_new = j * (opt.nx + 1);
620 for (
int i = 0; i <
GetNE(); i++)
625 for (
int j = 0; j < nv; j++)
631 for (
int i = 0; i <
GetNBE(); i++)
636 for (
int j = 0; j < nv; j++)
651 p[0] = cos(
u)*(1.0 + 0.3*sin(3.*
u +
PI*v));
652 p[1] = sin(
u)*(1.0 + 0.3*sin(3.*
u +
PI*v));
660#define I cdouble(0.0, 1.0)
747 if (abs(q2n) <
EPS) { q2n = 0.0; }
770static real_t ALPHA[4] {0.0};
771struct Costa:
public Surface
773 Costa(Opt &opt): Surface((opt.Tptr = Parametrization, opt), false) { }
775 void Prefix()
override
778 const int nx = opt.nx, ny = opt.ny;
779 MFEM_VERIFY(nx>2 && ny>2,
"");
780 const int nXhalf = (nx%2)==0 ? 4 : 2;
781 const int nYhalf = (ny%2)==0 ? 4 : 2;
782 const int nxh = nXhalf + nYhalf;
783 const int NVert = (nx+1) * (ny+1);
784 const int NElem = nx*ny - 4 - nxh;
785 const int NBdrElem = 0;
788 for (
int j = 0; j <= ny; j++)
791 for (
int i = 0; i <= nx; i++)
799 for (
int j = 0; j < ny; j++)
801 for (
int i = 0; i < nx; i++)
803 if (i == 0 && j == 0) {
continue; }
804 if (i+1 == nx && j == 0) {
continue; }
805 if (i == 0 && j+1 == ny) {
continue; }
806 if (i+1 == nx && j+1 == ny) {
continue; }
807 if ((j == 0 || j+1 == ny) && (abs(nx-(i<<1)-1)<=1)) {
continue; }
808 if ((i == 0 || i+1 == nx) && (abs(ny-(j<<1)-1)<=1)) {
continue; }
809 const int i0 = i + j*(nx+1);
810 const int i1 = i+1 + j*(nx+1);
811 const int i2 = i+1 + (j+1)*(nx+1);
812 const int i3 = i + (j+1)*(nx+1);
813 const int ind[4] = {i0, i1, i2, i3};
825 const real_t tau = ALPHA[3];
832 const bool y_top = x[1] > 0.5;
833 const bool x_top = x[0] > 0.5;
836 if (y_top) { v = 1.0 - x[1]; }
837 if (x_top) {
u = 1.0 - x[0]; }
847 p[2] = sqrt(
PI/2.)*log(abs((pw-e1)/(pw+e1)));
848 if (y_top) {
p[1] *= -1.0; }
849 if (x_top) {
p[0] *= -1.0; }
850 const bool nan = std::isnan(
p[0]) || std::isnan(
p[1]) || std::isnan(
p[2]);
851 MFEM_VERIFY(!nan,
"nan");
852 ALPHA[0] = std::max(
p[0], ALPHA[0]);
853 ALPHA[1] = std::max(
p[1], ALPHA[1]);
854 ALPHA[2] = std::max(
p[2], ALPHA[2]);
860 MFEM_VERIFY(ALPHA[0] > 0.0,
"");
861 MFEM_VERIFY(ALPHA[1] > 0.0,
"");
862 MFEM_VERIFY(ALPHA[2] > 0.0,
"");
864 const real_t phi = (1.0 + sqrt(5.0)) / 2.0;
865 for (
int i = 0; i <
nodes.FESpace()->GetNDofs(); i++)
867 for (
int d = 0; d <
SDIM; d++)
870 const int vdof =
nodes.FESpace()->DofToVDof(i, d);
878struct Shell:
public Surface
880 Shell(Opt &opt): Surface((opt.niters = 1, opt.Tptr = Parametrization, opt)) { }
887 const real_t v = 21.0*x[1]-15.0;
888 p[0] = +1.0*pow(1.16,v)*cos(v)*(1.0+cos(
u));
889 p[1] = -1.0*pow(1.16,v)*sin(v)*(1.0+cos(
u));
890 p[2] = -2.0*pow(1.16,v)*(1.0+sin(
u));
895struct Scherk:
public Surface
906 p[2] = log(cos(v)/cos(
u));
909 Scherk(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
913struct FullPeach:
public Surface
915 static constexpr int NV = 8;
916 static constexpr int NE = 6;
919 Surface((opt.niters = std::min(4, opt.niters), opt), NV, NE, 0) { }
921 void Prefix()
override
925 {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
926 {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
928 const int quad_e[NE][4] =
930 {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
931 {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
934 for (
int j = 0; j < NV; j++) {
AddVertex(quad_v[j]); }
935 for (
int j = 0; j < NE; j++) {
AddQuad(quad_e[j], j+1); }
943 void Snap()
override { SnapNodesToUnitSphere(); }
945 void BoundaryConditions()
override
955 for (
int e = 0; e < fes->
GetNE(); e++)
962 for (
int dof = 0; dof < dofs.
Size(); dof++)
966 const int k = dofs[dof];
967 MFEM_ASSERT(k >= 0,
"");
968 PointMat.
Mult(one, X);
969 const bool halfX = std::abs(X[0]) <
EPS && X[1] <= 0.0;
970 const bool halfY = std::abs(X[2]) <
EPS && X[1] >= 0.0;
971 const bool is_on_bc = halfX || halfY;
972 for (
int c = 0; c <
SDIM; c++)
973 { ess_vdofs[fes->
DofToVDof(k, c)] = is_on_bc; }
991struct QuarterPeach:
public Surface
993 QuarterPeach(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
998 const real_t x = 2.0*X[0]-1.0;
1000 const real_t r = sqrt(x*x + y*y);
1001 const real_t t = (x==0.0) ?
PI/2.0 :
1002 (y==0.0 && x>0.0) ? 0. :
1003 (y==0.0 && x<0.0) ?
PI : acos(x/r);
1004 const real_t sqrtx = sqrt(1.0 + x*x);
1005 const real_t sqrty = sqrt(1.0 + y*y);
1006 const bool yaxis =
PI/4.0<t && t < 3.0*
PI/4.0;
1007 const real_t R = yaxis?sqrtx:sqrty;
1008 const real_t gamma = r/R;
1009 p[0] = gamma * cos(t);
1010 p[1] = gamma * sin(t);
1014 void Postfix()
override
1016 for (
int i = 0; i <
GetNBE(); i++)
1026 for (
int v = 0; v < 2; v++)
1030 for (
int d = 0; d <
SDIM; d++)
1032 nodes->GetNodalValues(nval, d+1);
1033 const real_t x = X[v][d] = nval[iv];
1034 if (d < 2) { R[v] += x*x; }
1037 if (std::abs(X[0][1])<=
EPS && std::abs(X[1][1])<=
EPS &&
1038 (R[0]>0.1 || R[1]>0.1))
1046struct SlottedSphere:
public Surface
1048 SlottedSphere(Opt &opt): Surface((opt.niters = 4, opt), 64, 40, 0) { }
1050 void Prefix()
override
1053 constexpr int NV1D = 4;
1054 constexpr int NV = NV1D*NV1D*NV1D;
1055 constexpr int NEPF = (NV1D-1)*(NV1D-1);
1056 constexpr int NE = NEPF*6;
1059 for (
int iv=0; iv<NV; ++iv)
1062 int iy = (iv / NV1D) % NV1D;
1063 int iz = (iv / NV1D) / NV1D;
1065 QV[iv][0] = V1D[ix];
1066 QV[iv][1] = V1D[iy];
1067 QV[iv][2] = V1D[iz];
1070 for (
int ix=0; ix<NV1D-1; ++ix)
1072 for (
int iy=0; iy<NV1D-1; ++iy)
1074 int el_offset = ix + iy*(NV1D-1);
1076 QE[0*NEPF + el_offset][0] = NV1D*ix + NV1D*NV1D*iy;
1077 QE[0*NEPF + el_offset][1] = NV1D*(ix+1) + NV1D*NV1D*iy;
1078 QE[0*NEPF + el_offset][2] = NV1D*(ix+1) + NV1D*NV1D*(iy+1);
1079 QE[0*NEPF + el_offset][3] = NV1D*ix + NV1D*NV1D*(iy+1);
1082 QE[1*NEPF + el_offset][3] = x_off + NV1D*ix + NV1D*NV1D*iy;
1083 QE[1*NEPF + el_offset][2] = x_off + NV1D*(ix+1) + NV1D*NV1D*iy;
1084 QE[1*NEPF + el_offset][1] = x_off + NV1D*(ix+1) + NV1D*NV1D*(iy+1);
1085 QE[1*NEPF + el_offset][0] = x_off + NV1D*ix + NV1D*NV1D*(iy+1);
1087 QE[2*NEPF + el_offset][0] = NV1D*NV1D*iy + ix;
1088 QE[2*NEPF + el_offset][1] = NV1D*NV1D*iy + ix + 1;
1089 QE[2*NEPF + el_offset][2] = NV1D*NV1D*(iy+1) + ix + 1;
1090 QE[2*NEPF + el_offset][3] = NV1D*NV1D*(iy+1) + ix;
1092 int y_off = NV1D*(NV1D-1);
1093 QE[3*NEPF + el_offset][0] = y_off + NV1D*NV1D*iy + ix;
1094 QE[3*NEPF + el_offset][1] = y_off + NV1D*NV1D*iy + ix + 1;
1095 QE[3*NEPF + el_offset][2] = y_off + NV1D*NV1D*(iy+1) + ix + 1;
1096 QE[3*NEPF + el_offset][3] = y_off + NV1D*NV1D*(iy+1) + ix;
1098 QE[4*NEPF + el_offset][0] = NV1D*iy + ix;
1099 QE[4*NEPF + el_offset][1] = NV1D*iy + ix + 1;
1100 QE[4*NEPF + el_offset][2] = NV1D*(iy+1) + ix + 1;
1101 QE[4*NEPF + el_offset][3] = NV1D*(iy+1) + ix;
1103 int z_off = NV1D*NV1D*(NV1D-1);
1104 QE[5*NEPF + el_offset][0] = z_off + NV1D*iy + ix;
1105 QE[5*NEPF + el_offset][1] = z_off + NV1D*iy + ix + 1;
1106 QE[5*NEPF + el_offset][2] = z_off + NV1D*(iy+1) + ix + 1;
1107 QE[5*NEPF + el_offset][3] = z_off + NV1D*(iy+1) + ix;
1111 QE[0*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1112 QE[0*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1114 QE[1*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1115 QE[1*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1117 QE[3*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1118 QE[3*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1120 QE[5*NEPF + 0 + 1*(NV1D-1)][0] = -1;
1121 QE[5*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1122 QE[5*NEPF + 2 + 1*(NV1D-1)][0] = -1;
1124 QE[4*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1125 QE[4*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1126 QE[4*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1128 QE[2*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1129 QE[2*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1130 for (
int j = 0; j < NV; j++) {
AddVertex(QV[j]); }
1131 for (
int j = 0; j < NE; j++)
1133 if (QE[j][0] < 0) {
continue; }
1142 void Snap()
override { SnapNodesToUnitSphere(); }
1145static int Problem0(Opt &opt)
1148 Surface *S =
nullptr;
1149 switch (opt.surface)
1151 case 0: S =
new MeshFromFile(opt);
break;
1152 case 1: S =
new Catenoid(opt);
break;
1153 case 2: S =
new Helicoid(opt);
break;
1154 case 3: S =
new Enneper(opt);
break;
1155 case 4: S =
new Hold(opt);
break;
1156 case 5: S =
new Costa(opt);
break;
1157 case 6: S =
new Shell(opt);
break;
1158 case 7: S =
new Scherk(opt);
break;
1159 case 8: S =
new FullPeach(opt);
break;
1160 case 9: S =
new QuarterPeach(opt);
break;
1161 case 10: S =
new SlottedSphere(opt);
break;
1162 default: MFEM_ABORT(
"Unknown surface (surface <= 10)!");
1172static real_t u0(
const Vector &x) {
return sin(3.0 *
PI * (x[1] + x[0])); }
1176static real_t qf(
const int order,
const int ker,
Mesh &m,
1183 const int NE(m.
GetNE());
1191 MFEM_VERIFY(ND == D1D*D1D,
"");
1192 MFEM_VERIFY(NQ == Q1D*Q1D,
"");
1194 Vector Eu(ND*NE), grad_u(
DIM*NQ*NE), sum(NE*NQ), one(NE*NQ);
1195 qi->SetOutputLayout(QVectorLayout::byVDIM);
1199 qi->Derivatives(Eu, grad_u);
1204 auto grdU =
Reshape(grad_u.Read(),
DIM, Q1D, Q1D, NE);
1205 auto S =
Reshape(sum.Write(), Q1D, Q1D, NE);
1209 MFEM_FOREACH_THREAD(qy,y,Q1D)
1211 MFEM_FOREACH_THREAD(qx,x,Q1D)
1213 const real_t w = W(qx, qy);
1214 const real_t J11 = J(qx, qy, 0, 0, e);
1215 const real_t J12 = J(qx, qy, 1, 0, e);
1216 const real_t J21 = J(qx, qy, 0, 1, e);
1217 const real_t J22 = J(qx, qy, 1, 1, e);
1218 const real_t det = detJ(qx, qy, e);
1219 const real_t area = w * det;
1220 const real_t gu0 = grdU(0, qx, qy, e);
1221 const real_t gu1 = grdU(1, qx, qy, e);
1222 const real_t tgu0 = (J22*gu0 - J12*gu1)/det;
1223 const real_t tgu1 = (J11*gu1 - J21*gu0)/det;
1224 const real_t ngu = tgu0*tgu0 + tgu1*tgu1;
1225 const real_t s = (ker ==
AREA) ? sqrt(1.0 + ngu) :
1226 (ker ==
NORM) ? ngu : 0.0;
1227 S(qx, qy, e) = area * s;
1235static int Problem1(Opt &opt)
1237 const int order = opt.order;
1252 u.ProjectCoefficient(u0_fc);
1254 if (opt.vis) { opt.vis = glvis.
open(
vishost, opt.visport) == 0; }
1255 if (opt.vis) { Surface::Visualize(glvis, opt, &mesh,
GLVIZ_W,
GLVIZ_H, &
u); }
1264 for (
int i = 0; i < opt.niters; i++)
1269 if (opt.pa) {
a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
1270 const real_t q_uold = qf(order,
AREA, mesh, fes, uold);
1271 MFEM_VERIFY(std::abs(q_uold) >
EPS,
"");
1275 a.FormLinearSystem(ess_tdof_list,
u,
b, A, X, B);
1278 a.RecoverFEMSolution(X,
b,
u);
1280 const real_t norm = sqrt(std::abs(qf(order,
NORM, mesh, fes, eps)));
1281 const real_t area = qf(order,
AREA, mesh, fes,
u);
1284 mfem::out <<
"Iteration " << i <<
", norm: " << norm
1285 <<
", area: " << area << std::endl;
1287 if (opt.vis) { Surface::Visualize(glvis, opt, &mesh, &
u); }
1288 if (opt.print) { Surface::Print(opt, &mesh, &
u); }
1289 if (norm <
NRM) {
break; }
1297 opt.sz = opt.id = 0;
1300 args.
AddOption(&opt.visport,
"-p",
"--send-port",
"Socket for GLVis.");
1301 args.
AddOption(&opt.pb,
"-p",
"--problem",
"Problem to solve.");
1302 args.
AddOption(&opt.mesh_file,
"-m",
"--mesh",
"Mesh file to use.");
1303 args.
AddOption(&opt.wait,
"-w",
"--wait",
"-no-w",
"--no-wait",
1304 "Enable or disable a GLVis pause.");
1305 args.
AddOption(&opt.radial,
"-rad",
"--radial",
"-no-rad",
"--no-radial",
1306 "Enable or disable radial constraints in solver.");
1307 args.
AddOption(&opt.nx,
"-x",
"--num-elements-x",
1308 "Number of elements in x-direction.");
1309 args.
AddOption(&opt.ny,
"-y",
"--num-elements-y",
1310 "Number of elements in y-direction.");
1311 args.
AddOption(&opt.order,
"-o",
"--order",
"Finite element order.");
1312 args.
AddOption(&opt.refine,
"-r",
"--ref-levels",
"Refinement");
1313 args.
AddOption(&opt.niters,
"-n",
"--niter-max",
"Max number of iterations");
1314 args.
AddOption(&opt.surface,
"-s",
"--surface",
"Choice of the surface.");
1315 args.
AddOption(&opt.pa,
"-pa",
"--partial-assembly",
"-no-pa",
1316 "--no-partial-assembly",
"Enable Partial Assembly.");
1317 args.
AddOption(&opt.tau,
"-t",
"--tau",
"Costa scale factor.");
1318 args.
AddOption(&opt.lambda,
"-l",
"--lambda",
"Lambda step toward solution.");
1319 args.
AddOption(&opt.amr,
"-a",
"--amr",
"-no-a",
"--no-amr",
"Enable AMR.");
1320 args.
AddOption(&opt.amr_threshold,
"-at",
"--amr-threshold",
"AMR threshold.");
1321 args.
AddOption(&opt.device_config,
"-d",
"--device",
1322 "Device configuration string, see Device::Configure().");
1323 args.
AddOption(&opt.keys,
"-k",
"--keys",
"GLVis configuration keys.");
1324 args.
AddOption(&opt.vis,
"-vis",
"--visualization",
"-no-vis",
1325 "--no-visualization",
"Enable or disable visualization.");
1326 args.
AddOption(&opt.vis_mesh,
"-vm",
"--vis-mesh",
"-no-vm",
1327 "--no-vis-mesh",
"Enable or disable mesh visualization.");
1328 args.
AddOption(&opt.by_vdim,
"-c",
"--solve-byvdim",
1329 "-no-c",
"--solve-bynodes",
1330 "Enable or disable the 'ByVdim' solver");
1331 args.
AddOption(&opt.print,
"-print",
"--print",
"-no-print",
"--no-print",
1332 "Enable or disable result output (files in mfem format).");
1333 args.
AddOption(&opt.snapshot,
"-ss",
"--snapshot",
"-no-ss",
"--no-snapshot",
1334 "Enable or disable GLVis snapshot.");
1337 MFEM_VERIFY(opt.lambda >= 0.0 && opt.lambda <= 1.0,
"");
1341 Device device(opt.device_config);
1342 if (!opt.id) { device.
Print(); }
1344 if (opt.pb == 0) { Problem0(opt); }
1346 if (opt.pb == 1) { Problem1(opt); }
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Conjugate gradient method.
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the Conjugate Gradient method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
A coefficient that is constant across space and time.
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void Transpose()
(*this) = (*this)^t
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Abstract data type element.
Geometry::Type GetGeometryType() const
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Type
Constants for the classes derived from Element.
virtual int GetNVertices() const =0
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int GetNE() const
Returns number of elements in the mesh.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Mesh * GetMesh() const
Returns the mesh.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
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 GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Vector detJ
Determinants of the Jacobians at all quadrature points.
Vector J
Jacobians of the element transformations at all quadrature points.
Class for grid function - Vector with associated FE space.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Arbitrary order H1-conforming (continuous) finite elements.
Class for an integration rule - an Array of IntegrationPoint.
int GetOrder() const
Returns the order of the integration rule.
int GetNPoints() const
Returns the number of the points in the integration rule.
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
void SetAbsTol(real_t atol)
void InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
Begin construction of a mesh.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally,...
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Geometry::Type GetTypicalElementGeometry() const
If the local mesh is not empty, return GetElementGeometry(0); otherwise, return a typical Geometry pr...
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
const Element * GetElement(int i) const
Return pointer to the i'th element object.
bool FaceIsTrueInterior(int FaceNo) const
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
int GetNE() const
Returns number of elements.
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
void GetNodes(Vector &node_coord) const
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
int GetNBE() const
Returns number of boundary elements.
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
void EnsureNCMesh(bool simplices_nonconforming=false)
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
void Transform(void(*f)(const Vector &, Vector &))
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void RemoveInternalBoundaries()
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Pointer to an Operator of a specified type.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
bool Good() const
Return true if the command line options were parsed successfully.
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
void SetSize(int s)
Resize the vector to size s.
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int close()
Close the socketstream.
real_t Surface()
Analytic surface integral.
cdouble LogEllipticTheta1Prime(const cdouble u, const cdouble q)
constexpr Element::Type QUAD
std::complex< real_t > cdouble
cdouble WeierstrassZeta(const cdouble z, const cdouble w1=0.5, const cdouble w3=real_t(0.5) *I)
cdouble EllipticTheta1Prime(const int k, const cdouble u, const cdouble q)
cdouble WeierstrassP(const cdouble z, const cdouble w1=0.5, const cdouble w3=real_t(0.5) *I)
cdouble EllipticTheta(const int a, const cdouble u, const cdouble q)
real_t u(const Vector &xvec)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void add(const Vector &v1, const Vector &v2, Vector &v)
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
void forall_2D(int N, int X, int Y, lambda &&body)
void subtract(const Vector &x, const Vector &y, Vector &z)
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
void forall(int N, lambda &&body)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
real_t p(const Vector &x, real_t t)
constexpr Element::Type QUAD
cdouble WeierstrassZeta(const cdouble z, const cdouble w1=0.5, const cdouble w3=real_t(0.5) *I)
cdouble WeierstrassP(const cdouble z, const cdouble w1=0.5, const cdouble w3=real_t(0.5) *I)
std::array< int, NCMesh::MaxFaceNodes > nodes