66 #include "../../general/forall.hpp"
71 constexpr
int DIM = 2;
73 constexpr
double PI = M_PI;
74 constexpr
double NRM = 1.e-4;
75 constexpr
double EPS = 1.e-14;
77 constexpr
double NL_DMAX = std::numeric_limits<double>::max();
103 bool by_vdim =
false;
104 bool snapshot =
false;
105 bool vis_mesh =
false;
108 double amr_threshold = 0.6;
109 const char *keys =
"gAm";
110 const char *device_config =
"cpu";
111 const char *mesh_file =
"../../data/mobius-strip.mesh";
115 class Surface:
public Mesh
125 Surface(Opt &opt,
const char *file):
Mesh(file, true), opt(opt) { }
128 Surface(Opt &opt,
bool):
Mesh(), opt(opt) { }
132 :
Mesh(
Mesh::MakeCartesian2D(opt.nx, opt.ny,
QUAD, true)), opt(opt) { }
135 Surface(Opt &opt,
int nv,
int ne,
int nbe):
149 if (opt.amr) { EnsureNCMesh(); }
150 mesh =
new Mesh(*
this,
true);
152 BoundaryConditions();
164 ByVDim(*
this, opt).Solve();
168 ByNodes(*
this, opt).Solve();
170 if (opt.vis && opt.snapshot)
173 Visualize(opt, mesh, mesh->GetNodes());
180 if (opt.vis) { glvis.
close(); }
181 delete mesh;
delete fec;
delete fes;
184 virtual void Prefix()
189 virtual void Transform() {
if (opt.Tptr) {
Mesh::Transform(opt.Tptr); } }
191 virtual void Postfix()
196 virtual void Refine()
198 for (
int l = 0; l < opt.refine; l++)
207 for (
int i = 0; i < nodes.
Size(); i++)
209 if (std::abs(nodes(i)) <
EPS)
216 void SnapNodesToUnitSphere()
222 for (
int d = 0; d <
SDIM; d++)
227 for (
int d = 0; d <
SDIM; d++)
234 virtual void BoundaryConditions()
236 if (bdr_attributes.Size())
241 fes->GetEssentialTrueDofs(ess_bdr, bc);
246 static void Visualize(Opt &opt,
const Mesh *mesh,
247 const int w,
const int h,
251 if (opt.vis_mesh) { glvis <<
"mesh\n" << *mesh; }
252 else { glvis <<
"solution\n" << *mesh << solution; }
254 glvis <<
"window_size " << w <<
" " << h <<
"\n";
255 glvis <<
"keys " << opt.keys <<
"\n";
257 if (opt.wait) { glvis <<
"pause\n"; }
262 static void Visualize(
const Opt &opt,
const Mesh *mesh,
266 if (opt.vis_mesh) { glvis <<
"mesh\n" << *mesh; }
267 else { glvis <<
"solution\n" << *mesh << solution; }
268 if (opt.wait) { glvis <<
"pause\n"; }
269 if (opt.snapshot && opt.keys) { glvis <<
"keys " << opt.keys <<
"\n"; }
276 const char *mesh_file =
"surface.mesh";
277 const char *sol_file =
"sol.gf";
280 mfem::out <<
"Printing " << mesh_file <<
", " << sol_file << std::endl;
283 std::ofstream mesh_ofs(mesh_file);
284 mesh_ofs.precision(8);
285 mesh->
Print(mesh_ofs);
288 std::ofstream sol_ofs(sol_file);
289 sol_ofs.precision(8);
306 const int print_iter = -1, max_num_iter = 2000;
307 const double RTOLERANCE =
EPS, ATOLERANCE =
EPS*
EPS;
309 Solver(Surface &S, Opt &opt): opt(opt), S(S), cg(),
310 a(S.fes), x(S.fes), x0(S.fes),
b(S.fes), one(1.0)
312 cg.SetRelTol(RTOLERANCE);
313 cg.SetAbsTol(ATOLERANCE);
314 cg.SetMaxIter(max_num_iter);
315 cg.SetPrintLevel(print_iter);
322 constexpr
bool converged =
true;
324 for (
int i=0; i < opt.niters; ++i)
326 if (opt.amr) { Amr(); }
327 if (opt.vis) { Surface::Visualize(opt, S.mesh); }
328 if (!opt.id) {
mfem::out <<
"Iteration " << i <<
": "; }
329 S.mesh->DeleteGeometricFactors();
332 if (Step() == converged) {
break; }
334 if (opt.print) {
Surface::Print(opt, S.mesh, S.mesh->GetNodes()); }
337 virtual bool Step() = 0;
340 bool Converged(
const double rnorm) {
return rnorm <
NRM; }
346 a.FormLinearSystem(S.bc, x,
b, A, X, B);
348 if (M) { cg.SetPreconditioner(*M); }
351 a.RecoverFEMSolution(X,
b, x);
352 const bool by_vdim = opt.by_vdim;
353 GridFunction *nodes = by_vdim ? &x0 : S.fes->GetMesh()->GetNodes();
357 if (!opt.id) {
mfem::out <<
"rnorm = " << rnorm << std::endl; }
358 const double lambda = opt.lambda;
361 MFEM_VERIFY(!opt.radial,
"'VDim solver can't use radial option!");
362 return Converged(rnorm);
369 Vector ni(SDIM), di(SDIM);
370 for (
int i = 0; i <
delta.Size()/
SDIM; i++)
373 const int ndof = S.fes->GetNDofs();
374 for (
int d = 0; d <
SDIM; d++)
376 ni(d) = (*nodes)(d*ndof + i);
377 di(d) =
delta(d*ndof + i);
380 const double ndotd = (ni*di) / (ni*ni);
383 for (
int d = 0; d <
SDIM; d++) {
delta(d*ndof + i) = di(d); }
388 add(lambda, *nodes, (1.0-lambda), x, x);
389 return Converged(rnorm);
394 MFEM_VERIFY(opt.amr_threshold >= 0.0 && opt.amr_threshold <= 1.0,
"");
397 const int NE = mesh->
GetNE();
399 for (
int e = 0; e < NE; e++)
407 for (
int q = 0; q < NQ; q++)
414 const double w = Jadjt.
Weight();
415 minW = std::fmin(minW, w);
416 maxW = std::fmax(maxW, w);
418 if (std::fabs(maxW) != 0.0)
420 const double rho = minW / maxW;
421 MFEM_VERIFY(rho <= 1.0,
"");
435 S.BoundaryConditions();
441 class ByNodes:
public Solver
444 ByNodes(Surface &S, Opt &opt):
Solver(S, opt)
449 x = *S.fes->GetMesh()->GetNodes();
450 bool converge = ParAXeqB();
452 return converge ?
true :
false;
457 class ByVDim:
public Solver
462 auto d_Xi = Xi.
Read();
463 auto d_nodes = S.fes->GetMesh()->GetNodes()->Write();
464 const int ndof = S.fes->GetNDofs();
465 MFEM_FORALL(i, ndof, d_nodes[c*ndof + i] = d_Xi[i]; );
470 auto d_Xi = Xi.
Write();
471 const int ndof = S.fes->GetNDofs();
472 auto d_nodes = S.fes->GetMesh()->GetNodes()->Read();
473 MFEM_FORALL(i, ndof, d_Xi[i] = d_nodes[c*ndof + i]; );
476 ByVDim(Surface &S, Opt &opt):
Solver(S, opt)
481 bool cvg[
SDIM] {
false};
482 for (
int c=0; c <
SDIM; ++c)
489 const bool converged = cvg[0] && cvg[1] && cvg[2];
490 return converged ?
true :
false;
496 struct MeshFromFile:
public Surface
498 MeshFromFile(Opt &opt): Surface(opt, opt.mesh_file) { }
502 struct Catenoid:
public Surface
504 Catenoid(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
510 for (
int i = 0; i < v2v.Size(); i++) { v2v[i] = i; }
512 for (
int j = 0; j <= opt.ny; j++)
514 const int v_old = opt.nx + j * (opt.nx + 1);
515 const int v_new = j * (opt.nx + 1);
519 for (
int i = 0; i < GetNE(); i++)
524 for (
int j = 0; j < nv; j++)
530 for (
int i = 0; i < GetNBE(); i++)
532 Element *el = GetBdrElement(i);
535 for (
int j = 0; j < nv; j++)
540 RemoveUnusedVertices();
541 RemoveInternalBoundaries();
548 const double u = 2.0*
PI*x[0];
549 const double v =
PI*(x[1]-0.5)/3.;
557 struct Helicoid:
public Surface
559 Helicoid(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
565 const double u = 2.0*
PI*x[0];
566 const double v = 2.0*
PI*(2.0*x[1]-1.0)/3.0;
574 struct Enneper:
public Surface
576 Enneper(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
582 const double u = 4.0*(x[0]-0.5);
583 const double v = 4.0*(x[1]-0.5);
584 p[0] = +u - u*u*u/3.0 + u*v*v;
585 p[1] = -v - u*u*v + v*v*v/3.0;
591 struct Hold:
public Surface
593 Hold(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
599 for (
int i = 0; i < v2v.Size(); i++) { v2v[i] = i; }
601 for (
int j = 0; j <= opt.ny; j++)
603 const int v_old = opt.nx + j * (opt.nx + 1);
604 const int v_new = j * (opt.nx + 1);
608 for (
int i = 0; i < GetNE(); i++)
613 for (
int j = 0; j < nv; j++)
619 for (
int i = 0; i < GetNBE(); i++)
621 Element *el = GetBdrElement(i);
624 for (
int j = 0; j < nv; j++)
629 RemoveUnusedVertices();
630 RemoveInternalBoundaries();
637 const double u = 2.0*
PI*x[0];
638 const double v = x[1];
639 p[0] = cos(u)*(1.0 + 0.3*sin(3.*u +
PI*v));
640 p[1] = sin(u)*(1.0 + 0.3*sin(3.*u +
PI*v));
648 #define I cdouble(0.0, 1.0)
654 double delta = std::numeric_limits<double>::max();
658 for (
int n=0; delta >
EPS; n+=1)
660 const cdouble j = pow(-1,n)*pow(q,n*(n+1.0))*sin((2.0*n+1.0)*u);
664 return 2.0*pow(q,0.25)*J;
667 for (
int n=0; delta >
EPS; n+=1)
669 const cdouble j = pow(q,n*(n+1))*cos((2.0*n+1.0)*u);
673 return 2.0*pow(q,0.25)*J;
675 for (
int n=1; delta >
EPS; n+=1)
677 const cdouble j = pow(q,n*n)*cos(2.0*n*u);
683 for (
int n=1; delta >
EPS; n+=1)
685 const cdouble j =pow(-1,n)*pow(q,n*n)*cos(2.0*n*u);
700 const cdouble q = exp(I*M_PI*tau);
701 const cdouble e1 = M_PI*M_PI/(12.0*w1*w1)*
704 const cdouble u = M_PI*z / (2.0*w1);
713 double delta = std::numeric_limits<double>::max();
714 for (
int n=0; delta >
EPS; n+=1)
716 const double alpha = 2.0*n+1.0;
717 const cdouble Dcosine = pow(alpha,k)*sin(k*M_PI/2.0 + alpha*u);
718 const cdouble j = pow(-1,n)*pow(q,n*(n+1.0))*Dcosine;
722 return 2.0*pow(q,0.25)*J;
729 double delta = std::numeric_limits<double>::max();
730 for (
int n=1; delta >
EPS; n+=1)
733 if (abs(q2n) < EPS) { q2n = 0.0; }
734 const cdouble j = q2n/(1.0-q2n)*sin(2.0*n*u);
738 return 1.0/tan(u) + 4.0*J;
747 const cdouble q = exp(I*M_PI*tau);
748 const cdouble n1 = -M_PI*M_PI/(12.0*w1) *
751 const cdouble u = M_PI*z / (2.0*w1);
756 static double ALPHA[4] {0.0};
757 struct Costa:
public Surface
759 Costa(Opt &opt): Surface((opt.Tptr = Parametrization, opt), false) { }
764 const int nx = opt.nx, ny = opt.ny;
765 MFEM_VERIFY(nx>2 && ny>2,
"");
766 const int nXhalf = (nx%2)==0 ? 4 : 2;
767 const int nYhalf = (ny%2)==0 ? 4 : 2;
768 const int nxh = nXhalf + nYhalf;
769 const int NVert = (nx+1) * (ny+1);
770 const int NElem = nx*ny - 4 - nxh;
771 const int NBdrElem = 0;
772 InitMesh(
DIM, SDIM, NVert, NElem, NBdrElem);
774 for (
int j = 0; j <= ny; j++)
776 const double cy = ((double) j / ny) ;
777 for (
int i = 0; i <= nx; i++)
779 const double cx = ((double) i / nx);
780 const double coords[
SDIM] = {cx, cy, 0.0};
785 for (
int j = 0; j < ny; j++)
787 for (
int i = 0; i < nx; i++)
789 if (i == 0 && j == 0) {
continue; }
790 if (i+1 == nx && j == 0) {
continue; }
791 if (i == 0 && j+1 == ny) {
continue; }
792 if (i+1 == nx && j+1 == ny) {
continue; }
793 if ((j == 0 || j+1 == ny) && (abs(nx-(i<<1)-1)<=1)) {
continue; }
794 if ((i == 0 || i+1 == nx) && (abs(ny-(j<<1)-1)<=1)) {
continue; }
795 const int i0 = i + j*(nx+1);
796 const int i1 = i+1 + j*(nx+1);
797 const int i2 = i+1 + (j+1)*(nx+1);
798 const int i3 = i + (j+1)*(nx+1);
799 const int ind[4] = {i0, i1, i2, i3};
803 RemoveUnusedVertices();
804 FinalizeQuadMesh(
false, 0,
true);
811 const double tau = ALPHA[3];
818 const bool y_top = x[1] > 0.5;
819 const bool x_top = x[0] > 0.5;
822 if (y_top) { v = 1.0 - x[1]; }
823 if (x_top) { u = 1.0 - x[0]; }
831 p[0] = real(
PI*(u+
PI/(4.*e1))- zw +
PI/(2.*e1)*(dw));
832 p[1] = real(
PI*(v+
PI/(4.*e1))-I*zw-
PI*I/(2.*e1)*(dw));
833 p[2] = sqrt(
PI/2.)*log(abs((pw-e1)/(pw+e1)));
834 if (y_top) { p[1] *= -1.0; }
835 if (x_top) { p[0] *= -1.0; }
836 const bool nan = std::isnan(p[0]) || std::isnan(p[1]) || std::isnan(p[2]);
837 MFEM_VERIFY(!nan,
"nan");
838 ALPHA[0] = std::fmax(p[0], ALPHA[0]);
839 ALPHA[1] = std::fmax(p[1], ALPHA[1]);
840 ALPHA[2] = std::fmax(p[2], ALPHA[2]);
846 MFEM_VERIFY(ALPHA[0] > 0.0,
"");
847 MFEM_VERIFY(ALPHA[1] > 0.0,
"");
848 MFEM_VERIFY(ALPHA[2] > 0.0,
"");
850 const double phi = (1.0 + sqrt(5.0)) / 2.0;
853 for (
int d = 0; d <
SDIM; d++)
855 const double alpha = d==2 ? phi : 1.0;
857 nodes(vdof) /= alpha * ALPHA[d];
864 struct Shell:
public Surface
866 Shell(Opt &opt): Surface((opt.niters = 1, opt.Tptr = Parametrization, opt)) { }
872 const double u = 2.0*
PI*x[0];
873 const double v = 21.0*x[1]-15.0;
874 p[0] = +1.0*pow(1.16,v)*cos(v)*(1.0+cos(u));
875 p[1] = -1.0*pow(1.16,v)*sin(v)*(1.0+cos(u));
876 p[2] = -2.0*pow(1.16,v)*(1.0+sin(u));
881 struct Scherk:
public Surface
886 const double alpha = 0.49;
888 const double u = alpha*
PI*(2.0*x[0]-1.0);
889 const double v =
alpha*
PI*(2.0*x[1]-1.0);
892 p[2] = log(cos(v)/cos(u));
895 Scherk(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
899 struct FullPeach:
public Surface
901 static constexpr
int NV = 8;
902 static constexpr
int NE = 6;
905 Surface((opt.niters = std::min(4, opt.niters), opt), NV, NE, 0) { }
909 const double quad_v[NV][
SDIM] =
911 {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
912 {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
914 const int quad_e[NE][4] =
916 {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
917 {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
920 for (
int j = 0; j < NV; j++) { AddVertex(quad_v[j]); }
921 for (
int j = 0; j < NE; j++) { AddQuad(quad_e[j], j+1); }
923 FinalizeQuadMesh(
false, 0,
true);
924 FinalizeTopology(
false);
929 void Snap() { SnapNodesToUnitSphere(); }
931 void BoundaryConditions()
936 ess_vdofs.
SetSize(fes->GetVSize());
937 MFEM_VERIFY(fes->GetVSize() >= fes->GetTrueVSize(),
"");
941 for (
int e = 0; e < fes->GetNE(); e++)
943 fes->GetElementDofs(e, dofs);
948 for (
int dof = 0; dof < dofs.
Size(); dof++)
952 const int k = dofs[dof];
953 MFEM_ASSERT(k >= 0,
"");
954 PointMat.
Mult(one, X);
955 const bool halfX = std::fabs(X[0]) <
EPS && X[1] <= 0.0;
956 const bool halfY = std::fabs(X[2]) <
EPS && X[1] >= 0.0;
957 const bool is_on_bc = halfX || halfY;
958 for (
int c = 0; c <
SDIM; c++)
959 { ess_vdofs[fes->DofToVDof(k, c)] = is_on_bc; }
977 struct QuarterPeach:
public Surface
979 QuarterPeach(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
984 const double x = 2.0*X[0]-1.0;
985 const double y = X[1];
986 const double r = sqrt(x*x + y*y);
987 const double t = (x==0.0) ?
PI/2.0 :
988 (y==0.0 && x>0.0) ? 0. :
989 (y==0.0 && x<0.0) ?
PI : acos(x/r);
990 const double sqrtx = sqrt(1.0 + x*x);
991 const double sqrty = sqrt(1.0 + y*y);
992 const bool yaxis =
PI/4.0<t && t < 3.0*
PI/4.0;
993 const double R = yaxis?sqrtx:sqrty;
994 const double gamma = r/R;
995 p[0] = gamma * cos(t);
996 p[1] = gamma * sin(t);
1002 for (
int i = 0; i < GetNBE(); i++)
1004 Element *el = GetBdrElement(i);
1005 const int fn = GetBdrElementEdgeIndex(i);
1006 MFEM_VERIFY(!FaceIsTrueInterior(fn),
"");
1008 GetFaceVertices(fn, vertices);
1011 double R[2], X[2][
SDIM];
1012 for (
int v = 0; v < 2; v++)
1015 const int iv = vertices[v];
1016 for (
int d = 0; d <
SDIM; d++)
1019 const double x = X[v][d] = nval[iv];
1020 if (d < 2) { R[v] += x*x; }
1023 if (std::fabs(X[0][1])<=
EPS && std::fabs(X[1][1])<=
EPS &&
1024 (R[0]>0.1 || R[1]>0.1))
1032 struct SlottedSphere:
public Surface
1034 SlottedSphere(Opt &opt): Surface((opt.niters = 4, opt), 64, 40, 0) { }
1038 constexpr
double delta = 0.15;
1039 constexpr
int NV1D = 4;
1040 constexpr
int NV = NV1D*NV1D*NV1D;
1041 constexpr
int NEPF = (NV1D-1)*(NV1D-1);
1042 constexpr
int NE = NEPF*6;
1043 const double V1D[NV1D] = {-1.0, -
delta,
delta, 1.0};
1045 for (
int iv=0; iv<NV; ++iv)
1048 int iy = (iv / NV1D) % NV1D;
1049 int iz = (iv / NV1D) / NV1D;
1051 QV[iv][0] = V1D[ix];
1052 QV[iv][1] = V1D[iy];
1053 QV[iv][2] = V1D[iz];
1056 for (
int ix=0; ix<NV1D-1; ++ix)
1058 for (
int iy=0; iy<NV1D-1; ++iy)
1060 int el_offset = ix + iy*(NV1D-1);
1062 QE[0*NEPF + el_offset][0] = NV1D*ix + NV1D*NV1D*iy;
1063 QE[0*NEPF + el_offset][1] = NV1D*(ix+1) + NV1D*NV1D*iy;
1064 QE[0*NEPF + el_offset][2] = NV1D*(ix+1) + NV1D*NV1D*(iy+1);
1065 QE[0*NEPF + el_offset][3] = NV1D*ix + NV1D*NV1D*(iy+1);
1068 QE[1*NEPF + el_offset][3] = x_off + NV1D*ix + NV1D*NV1D*iy;
1069 QE[1*NEPF + el_offset][2] = x_off + NV1D*(ix+1) + NV1D*NV1D*iy;
1070 QE[1*NEPF + el_offset][1] = x_off + NV1D*(ix+1) + NV1D*NV1D*(iy+1);
1071 QE[1*NEPF + el_offset][0] = x_off + NV1D*ix + NV1D*NV1D*(iy+1);
1073 QE[2*NEPF + el_offset][0] = NV1D*NV1D*iy + ix;
1074 QE[2*NEPF + el_offset][1] = NV1D*NV1D*iy + ix + 1;
1075 QE[2*NEPF + el_offset][2] = NV1D*NV1D*(iy+1) + ix + 1;
1076 QE[2*NEPF + el_offset][3] = NV1D*NV1D*(iy+1) + ix;
1078 int y_off = NV1D*(NV1D-1);
1079 QE[3*NEPF + el_offset][0] = y_off + NV1D*NV1D*iy + ix;
1080 QE[3*NEPF + el_offset][1] = y_off + NV1D*NV1D*iy + ix + 1;
1081 QE[3*NEPF + el_offset][2] = y_off + NV1D*NV1D*(iy+1) + ix + 1;
1082 QE[3*NEPF + el_offset][3] = y_off + NV1D*NV1D*(iy+1) + ix;
1084 QE[4*NEPF + el_offset][0] = NV1D*iy + ix;
1085 QE[4*NEPF + el_offset][1] = NV1D*iy + ix + 1;
1086 QE[4*NEPF + el_offset][2] = NV1D*(iy+1) + ix + 1;
1087 QE[4*NEPF + el_offset][3] = NV1D*(iy+1) + ix;
1089 int z_off = NV1D*NV1D*(NV1D-1);
1090 QE[5*NEPF + el_offset][0] = z_off + NV1D*iy + ix;
1091 QE[5*NEPF + el_offset][1] = z_off + NV1D*iy + ix + 1;
1092 QE[5*NEPF + el_offset][2] = z_off + NV1D*(iy+1) + ix + 1;
1093 QE[5*NEPF + el_offset][3] = z_off + NV1D*(iy+1) + ix;
1097 QE[0*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1098 QE[0*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1100 QE[1*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1101 QE[1*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1103 QE[3*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1104 QE[3*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1106 QE[5*NEPF + 0 + 1*(NV1D-1)][0] = -1;
1107 QE[5*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1108 QE[5*NEPF + 2 + 1*(NV1D-1)][0] = -1;
1110 QE[4*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1111 QE[4*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1112 QE[4*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1114 QE[2*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1115 QE[2*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1116 for (
int j = 0; j < NV; j++) { AddVertex(QV[j]); }
1117 for (
int j = 0; j < NE; j++)
1119 if (QE[j][0] < 0) {
continue; }
1120 AddQuad(QE[j], j+1);
1122 RemoveUnusedVertices();
1123 FinalizeQuadMesh(
false, 0,
true);
1128 void Snap() { SnapNodesToUnitSphere(); }
1131 static int Problem0(Opt &opt)
1134 Surface *S =
nullptr;
1135 switch (opt.surface)
1137 case 0: S =
new MeshFromFile(opt);
break;
1138 case 1: S =
new Catenoid(opt);
break;
1139 case 2: S =
new Helicoid(opt);
break;
1140 case 3: S =
new Enneper(opt);
break;
1141 case 4: S =
new Hold(opt);
break;
1142 case 5: S =
new Costa(opt);
break;
1143 case 6: S =
new Shell(opt);
break;
1144 case 7: S =
new Scherk(opt);
break;
1145 case 8: S =
new FullPeach(opt);
break;
1146 case 9: S =
new QuarterPeach(opt);
break;
1147 case 10: S =
new SlottedSphere(opt);
break;
1148 default: MFEM_ABORT(
"Unknown surface (surface <= 10)!");
1158 static double u0(
const Vector &x) {
return sin(3.0 *
PI * (x[1] + x[0])); }
1162 static double qf(
const int order,
const int ker,
Mesh &m,
1169 const int NE(m.
GetNE());
1177 MFEM_VERIFY(ND == D1D*D1D,
"");
1178 MFEM_VERIFY(NQ == Q1D*Q1D,
"");
1180 Vector Eu(ND*NE), grad_u(
DIM*NQ*NE), sum(NE*NQ), one(NE*NQ);
1185 qi->Derivatives(Eu, grad_u);
1190 auto grdU =
Reshape(grad_u.Read(),
DIM, Q1D, Q1D, NE);
1191 auto S =
Reshape(sum.Write(), Q1D, Q1D, NE);
1193 MFEM_FORALL_2D(e, NE, Q1D, Q1D, 1,
1195 MFEM_FOREACH_THREAD(qy,y,Q1D)
1197 MFEM_FOREACH_THREAD(qx,x,Q1D)
1199 const double w = W(qx, qy);
1200 const double J11 = J(qx, qy, 0, 0, e);
1201 const double J12 = J(qx, qy, 1, 0, e);
1202 const double J21 = J(qx, qy, 0, 1, e);
1203 const double J22 = J(qx, qy, 1, 1, e);
1204 const double det = detJ(qx, qy, e);
1205 const double area = w * det;
1206 const double gu0 = grdU(0, qx, qy, e);
1207 const double gu1 = grdU(1, qx, qy, e);
1208 const double tgu0 = (J22*gu0 - J12*gu1)/det;
1209 const double tgu1 = (J11*gu1 - J21*gu0)/det;
1210 const double ngu = tgu0*tgu0 + tgu1*tgu1;
1211 const double s = (ker ==
AREA) ? sqrt(1.0 + ngu) :
1212 (ker ==
NORM) ? ngu : 0.0;
1213 S(qx, qy, e) = area *
s;
1221 static int Problem1(Opt &opt)
1223 const int order = opt.order;
1240 if (opt.vis) { Surface::Visualize(opt, &mesh,
GLVIZ_W,
GLVIZ_H, &u); }
1249 for (
int i = 0; i < opt.niters; i++)
1255 const double q_uold = qf(order,
AREA, mesh, fes, uold);
1256 MFEM_VERIFY(std::fabs(q_uold) > EPS,
"");
1260 a.FormLinearSystem(ess_tdof_list, u,
b, A, X, B);
1263 a.RecoverFEMSolution(X,
b, u);
1265 const double norm = sqrt(std::fabs(qf(order,
NORM, mesh, fes, eps)));
1266 const double area = qf(order,
AREA, mesh, fes, u);
1269 mfem::out <<
"Iteration " << i <<
", norm: " << norm
1270 <<
", area: " << area << std::endl;
1272 if (opt.vis) { Surface::Visualize(opt, &mesh, &u); }
1274 if (norm <
NRM) {
break; }
1282 opt.sz = opt.id = 0;
1285 args.
AddOption(&opt.pb,
"-p",
"--problem",
"Problem to solve.");
1286 args.
AddOption(&opt.mesh_file,
"-m",
"--mesh",
"Mesh file to use.");
1287 args.
AddOption(&opt.wait,
"-w",
"--wait",
"-no-w",
"--no-wait",
1288 "Enable or disable a GLVis pause.");
1289 args.
AddOption(&opt.radial,
"-rad",
"--radial",
"-no-rad",
"--no-radial",
1290 "Enable or disable radial constraints in solver.");
1291 args.
AddOption(&opt.nx,
"-x",
"--num-elements-x",
1292 "Number of elements in x-direction.");
1293 args.
AddOption(&opt.ny,
"-y",
"--num-elements-y",
1294 "Number of elements in y-direction.");
1295 args.
AddOption(&opt.order,
"-o",
"--order",
"Finite element order.");
1296 args.
AddOption(&opt.refine,
"-r",
"--ref-levels",
"Refinement");
1297 args.
AddOption(&opt.niters,
"-n",
"--niter-max",
"Max number of iterations");
1298 args.
AddOption(&opt.surface,
"-s",
"--surface",
"Choice of the surface.");
1299 args.
AddOption(&opt.pa,
"-pa",
"--partial-assembly",
"-no-pa",
1300 "--no-partial-assembly",
"Enable Partial Assembly.");
1301 args.
AddOption(&opt.tau,
"-t",
"--tau",
"Costa scale factor.");
1302 args.
AddOption(&opt.lambda,
"-l",
"--lambda",
"Lambda step toward solution.");
1303 args.
AddOption(&opt.amr,
"-a",
"--amr",
"-no-a",
"--no-amr",
"Enable AMR.");
1304 args.
AddOption(&opt.amr_threshold,
"-at",
"--amr-threshold",
"AMR threshold.");
1305 args.
AddOption(&opt.device_config,
"-d",
"--device",
1306 "Device configuration string, see Device::Configure().");
1307 args.
AddOption(&opt.keys,
"-k",
"--keys",
"GLVis configuration keys.");
1308 args.
AddOption(&opt.vis,
"-vis",
"--visualization",
"-no-vis",
1309 "--no-visualization",
"Enable or disable visualization.");
1310 args.
AddOption(&opt.vis_mesh,
"-vm",
"--vis-mesh",
"-no-vm",
1311 "--no-vis-mesh",
"Enable or disable mesh visualization.");
1312 args.
AddOption(&opt.by_vdim,
"-c",
"--solve-byvdim",
1313 "-no-c",
"--solve-bynodes",
1314 "Enable or disable the 'ByVdim' solver");
1315 args.
AddOption(&opt.print,
"-print",
"--print",
"-no-print",
"--no-print",
1316 "Enable or disable result output (files in mfem format).");
1317 args.
AddOption(&opt.snapshot,
"-ss",
"--snapshot",
"-no-ss",
"--no-snapshot",
1318 "Enable or disable GLVis snapshot.");
1321 MFEM_VERIFY(opt.lambda >= 0.0 && opt.lambda <= 1.0,
"");
1325 Device device(opt.device_config);
1326 if (!opt.id) { device.
Print(); }
1328 if (opt.pb == 0) { Problem0(opt); }
1330 if (opt.pb == 1) { Problem1(opt); }
int GetNPoints() const
Returns the number of the points in the integration rule.
Geometry::Type GetGeometryType() const
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
int Size() const
Return the logical size of the array.
virtual void Print(std::ostream &out=mfem::out) const
Conjugate gradient method.
int GetNDofs() const
Returns number of degrees of freedom.
Class for an integration rule - an Array of IntegrationPoint.
Class for grid function - Vector with associated FE space.
int DofToVDof(int dof, int vd, int ndofs=-1) const
cdouble LogEllipticTheta1Prime(const cdouble u, const cdouble q)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
A coefficient that is constant across space and time.
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.
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
cdouble WeierstrassZeta(const cdouble z, const cdouble w1=0.5, const cdouble w3=0.5 *I)
void SetSize(int s)
Resize the vector to size s.
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
double Norml2() const
Returns the l2 norm of the vector.
Pointer to an Operator of a specified type.
std::complex< double > cdouble
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
void Transform(void(*f)(const Vector &, Vector &))
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
int GetNE() const
Returns number of elements.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
constexpr Element::Type QUAD
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
Geometry::Type GetElementBaseGeometry(int i) const
void add(const Vector &v1, const Vector &v2, Vector &v)
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
int close()
Close the socketstream.
Data type for Gauss-Seidel smoother of sparse matrix.
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Vector J
Jacobians of the element transformations at all quadrature points.
void SetPrintLevel(int print_lvl)
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
cdouble WeierstrassP(const cdouble z, const cdouble w1=0.5, const cdouble w3=0.5 *I)
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Vector detJ
Determinants of the Jacobians at all quadrature points.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void SetMaxIter(int max_it)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Type
Constants for the classes derived from Element.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
cdouble EllipticTheta(const int a, const cdouble u, const cdouble q)
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
const Element * GetElement(int i) const
FiniteElementSpace * FESpace()
void PrintUsage(std::ostream &out) const
Print the usage message.
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double p(const Vector &x, double t)
void SetRelTol(double rtol)
void Transpose()
(*this) = (*this)^t
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
void subtract(const Vector &x, const Vector &y, Vector &z)
int GetOrder() const
Returns the order of the integration rule.
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...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, double sx=1.0, double sy=1.0, bool sfc_ordering=true)
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
void PrintOptions(std::ostream &out) const
Print the options.
Lexicographic ordering for tensor-product FiniteElements.
virtual int GetNVertices() const =0
void SetAttribute(const int attr)
Set element's attribute.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A general function coefficient.
void Mult(const double *x, double *y) const
Matrix vector multiplication.
void GetNodes(Vector &node_coord) const
Arbitrary order H1-conforming (continuous) finite elements.
double u(const Vector &xvec)
cdouble EllipticTheta1Prime(const int k, const cdouble u, const cdouble q)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Abstract data type element.
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
bool Good() const
Return true if the command line options were parsed successfully.