MFEM  v4.6.0 Finite element discretization library
pminimal-surface.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // --------------------------------------------------------------------
13 // Minimal Surface Miniapp: Compute minimal surfaces - Parallel Version
14 // --------------------------------------------------------------------
15 //
16 // This miniapp solves Plateau's problem: the Dirichlet problem for the minimal
17 // surface equation.
18 //
19 // Two problems can be run:
20 //
21 // - Problem 0 solves the minimal surface equation of parametric surfaces.
22 // The surface (-s) option allow the selection of different
23 // parametrization:
24 // s=0: Uses the given mesh from command line options
25 // s=1: Catenoid
26 // s=2: Helicoid
27 // s=3: Enneper
28 // s=4: Hold
29 // s=5: Costa
30 // s=6: Shell
31 // s=7: Scherk
32 // s=8: FullPeach
33 // s=9: QuarterPeach
34 // s=10: SlottedSphere
35 //
36 // - Problem 1 solves the minimal surface equation of the form z=f(x,y),
37 // for the Dirichlet problem, using Picard iterations:
38 // -div( q grad u) = 0, with q(u) = (1 + |∇u|²)^{-1/2}
39 //
40 // Compile with: make pminimal-surface
41 //
42 // Sample runs: mpirun -np 4 pminimal-surface
43 // mpirun -np 4 pminimal-surface -a
44 // mpirun -np 4 pminimal-surface -c
45 // mpirun -np 4 pminimal-surface -c -a
46 // mpirun -np 4 pminimal-surface -no-pa
47 // mpirun -np 4 pminimal-surface -no-pa -a
48 // mpirun -np 4 pminimal-surface -no-pa -a -c
49 // mpirun -np 4 pminimal-surface -p 1
50 //
51 // Device sample runs:
52 // mpirun -np 4 pminimal-surface -d debug
53 // mpirun -np 4 pminimal-surface -d debug -a
54 // mpirun -np 4 pminimal-surface -d debug -c
55 // mpirun -np 4 pminimal-surface -d debug -c -a
56 // mpirun -np 4 pminimal-surface -d cuda
57 // mpirun -np 4 pminimal-surface -d cuda -a
58 // mpirun -np 4 pminimal-surface -d cuda -c
59 // mpirun -np 4 pminimal-surface -d cuda -c -a
60 // mpirun -np 4 pminimal-surface -d cuda -no-pa
61 // mpirun -np 4 pminimal-surface -d cuda -no-pa -a
62 // mpirun -np 4 pminimal-surface -d cuda -no-pa -c
63 // mpirun -np 4 pminimal-surface -d cuda -no-pa -c -a
64
65 #include "mfem.hpp"
66 #include "../../general/forall.hpp"
67
68 using namespace mfem;
69
70 // Constant variables
71 constexpr int DIM = 2;
72 constexpr int SDIM = 3;
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();
78
79 // Static variables for GLVis
80 constexpr int GLVIZ_W = 1024;
81 constexpr int GLVIZ_H = 1024;
82 constexpr int visport = 19916;
83 constexpr char vishost[] = "localhost";
84
85 // Context/Options for the solver
86 struct Opt
87 {
88  int sz, id;
89  int pb = 0;
90  int nx = 6;
91  int ny = 6;
92  int order = 3;
93  int refine = 2;
94  int niters = 8;
95  int surface = 5;
96  bool pa = true;
97  bool vis = true;
98  bool amr = false;
99  bool wait = false;
100  bool print = false;
102  bool by_vdim = false;
103  bool snapshot = false;
104  // bool vis_mesh = false; // Not supported by GLVis
105  double tau = 1.0;
106  double lambda = 0.1;
107  double amr_threshold = 0.6;
108  const char *keys = "gAm";
109  const char *device_config = "cpu";
110  const char *mesh_file = "../../data/mobius-strip.mesh";
111  void (*Tptr)(const Vector&, Vector&) = nullptr;
112 };
113
114 class Surface: public Mesh
115 {
116 protected:
117  Opt &opt;
118  ParMesh *mesh;
119  Array<int> bc;
120  socketstream glvis;
121  H1_FECollection *fec;
123 public:
124  // Reading from mesh file
125  Surface(Opt &opt, const char *file): Mesh(file, true), opt(opt) { }
126
127  // Generate 2D empty surface mesh
128  Surface(Opt &opt, bool): Mesh(), opt(opt) { }
129
130  // Generate 2D quad surface mesh
131  Surface(Opt &opt)
132  : Mesh(Mesh::MakeCartesian2D(opt.nx, opt.ny, QUAD, true)), opt(opt) { }
133
134  // Generate 2D generic surface mesh
135  Surface(Opt &opt, int nv, int ne, int nbe):
136  Mesh(DIM, nv, ne, nbe, SDIM), opt(opt) { }
137
138  void Create()
139  {
140  if (opt.surface > 0)
141  {
142  Prefix();
143  Transform();
144  }
145  Postfix();
146  Refine();
147  Snap();
148  fec = new H1_FECollection(opt.order, DIM);
149  if (opt.amr) { EnsureNCMesh(); }
150  mesh = new ParMesh(MPI_COMM_WORLD, *this);
151  fes = new ParFiniteElementSpace(mesh, fec, opt.by_vdim ? 1 : SDIM);
152  BoundaryConditions();
153  }
154
155  int Solve()
156  {
157  // Initialize GLVis server if 'visualization' is set
158  if (opt.vis) { opt.vis = glvis.open(vishost, visport) == 0; }
159  // Send to GLVis the first mesh
160  if (opt.vis) { Visualize(glvis, opt, mesh, GLVIZ_W, GLVIZ_H); }
161  // Create and launch the surface solver
162  if (opt.by_vdim)
163  {
164  ByVDim(*this, opt).Solve();
165  }
166  else
167  {
168  ByNodes(*this, opt).Solve();
169  }
170  if (opt.vis && opt.snapshot)
171  {
172  opt.keys = "Sq";
173  Visualize(glvis, opt, mesh, mesh->GetNodes());
174  }
175  return 0;
176  }
177
178  ~Surface()
179  {
180  if (opt.vis) { glvis.close(); }
181  delete mesh; delete fec; delete fes;
182  }
183
184  virtual void Prefix()
185  {
186  SetCurvature(opt.order, false, SDIM, Ordering::byNODES);
187  }
188
189  virtual void Transform() { if (opt.Tptr) { Mesh::Transform(opt.Tptr); } }
190
191  virtual void Postfix()
192  {
193  SetCurvature(opt.order, false, SDIM, Ordering::byNODES);
194  }
195
196  virtual void Refine()
197  {
198  for (int l = 0; l < opt.refine; l++)
199  {
200  UniformRefinement();
201  }
202  }
203
204  virtual void Snap()
205  {
206  GridFunction &nodes = *GetNodes();
207  for (int i = 0; i < nodes.Size(); i++)
208  {
209  if (std::abs(nodes(i)) < EPS)
210  {
211  nodes(i) = 0.0;
212  }
213  }
214  }
215
216  void SnapNodesToUnitSphere()
217  {
218  Vector node(SDIM);
219  GridFunction &nodes = *GetNodes();
220  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
221  {
222  for (int d = 0; d < SDIM; d++)
223  {
224  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
225  }
226  node /= node.Norml2();
227  for (int d = 0; d < SDIM; d++)
228  {
229  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
230  }
231  }
232  }
233
234  virtual void BoundaryConditions()
235  {
236  if (bdr_attributes.Size())
237  {
238  Array<int> ess_bdr(bdr_attributes.Max());
239  ess_bdr = 1;
241  fes->GetEssentialTrueDofs(ess_bdr, bc);
242  }
243  }
244
245  // Initialize visualization of some given mesh
246  static void Visualize(socketstream &glvis,
247  Opt &opt, const Mesh *mesh,
248  const int w, const int h,
249  const GridFunction *sol = nullptr)
250  {
251  const GridFunction &solution = sol ? *sol : *mesh->GetNodes();
252  glvis << "parallel " << opt.sz << " " << opt.id << "\n";
253  glvis << "solution\n" << *mesh << solution;
254  glvis.precision(8);
255  glvis << "window_size " << w << " " << h << "\n";
256  glvis << "keys " << opt.keys << "\n";
257  opt.keys = nullptr;
258  if (opt.wait) { glvis << "pause\n"; }
259  glvis << std::flush;
260  }
261
262  // Visualize some solution on the given mesh
263  static void Visualize(socketstream &glvis,
264  const Opt &opt, const Mesh *mesh,
265  const GridFunction *sol = nullptr)
266  {
267  glvis << "parallel " << opt.sz << " " << opt.id << "\n";
268  const GridFunction &solution = sol ? *sol : *mesh->GetNodes();
269  glvis << "solution\n" << *mesh << solution;
270  if (opt.wait) { glvis << "pause\n"; }
271  if (opt.snapshot && opt.keys) { glvis << "keys " << opt.keys << "\n"; }
272  glvis << std::flush;
273  }
274
275  using Mesh::Print;
276  static void Print(const Opt &opt, ParMesh *mesh, const GridFunction *sol)
277  {
278  const char *mesh_file = "surface.mesh";
279  const char *sol_file = "sol.gf";
280  if (!opt.id)
281  {
282  mfem::out << "Printing " << mesh_file << ", " << sol_file << std::endl;
283  }
284
285  std::ostringstream mesh_name;
286  mesh_name << mesh_file << "." << std::setfill('0') << std::setw(6) << opt.id;
287  std::ofstream mesh_ofs(mesh_name.str().c_str());
288  mesh_ofs.precision(8);
289  mesh->Print(mesh_ofs);
290  mesh_ofs.close();
291
292  std::ostringstream sol_name;
293  sol_name << sol_file << "." << std::setfill('0') << std::setw(6) << opt.id;
294  std::ofstream sol_ofs(sol_name.str().c_str());
295  sol_ofs.precision(8);
296  sol->Save(sol_ofs);
297  sol_ofs.close();
298  }
299
300  // Surface Solver class
301  class Solver
302  {
303  protected:
304  Opt &opt;
305  Surface &S;
306  CGSolver cg;
307  OperatorPtr A;
309  ParGridFunction x, x0, b;
311  mfem::Solver *M = nullptr;
312  const int print_iter = -1, max_num_iter = 2000;
313  const double RTOLERANCE = EPS, ATOLERANCE = EPS*EPS;
314  public:
315  Solver(Surface &S, Opt &opt): opt(opt), S(S), cg(MPI_COMM_WORLD),
316  a(S.fes), x(S.fes), x0(S.fes), b(S.fes), one(1.0)
317  {
318  cg.SetRelTol(RTOLERANCE);
319  cg.SetAbsTol(ATOLERANCE);
320  cg.SetMaxIter(max_num_iter);
321  cg.SetPrintLevel(print_iter);
322  }
323
324  ~Solver() { delete M; }
325
326  void Solve()
327  {
328  constexpr bool converged = true;
329  if (opt.pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL);}
330  for (int i=0; i < opt.niters; ++i)
331  {
332  if (opt.amr) { Amr(); }
333  if (opt.vis) { Surface::Visualize(S.glvis, opt, S.mesh); }
334  if (!opt.id) { mfem::out << "Iteration " << i << ": "; }
335  S.mesh->NodesUpdated();
336  a.Update();
337  a.Assemble();
338  if (Step() == converged) { break; }
339  }
340  if (opt.print) { Surface::Print(opt, S.mesh, S.mesh->GetNodes()); }
341  }
342
343  virtual bool Step() = 0;
344
345  protected:
346  bool Converged(const double rnorm) { return rnorm < NRM; }
347
348  bool ParAXeqB()
349  {
350  b = 0.0;
351  Vector X, B;
352  a.FormLinearSystem(S.bc, x, b, A, X, B);
353  if (!opt.pa) { M = new HypreBoomerAMG; }
354  if (M) { cg.SetPreconditioner(*M); }
355  cg.SetOperator(*A);
356  cg.Mult(B, X);
357  a.RecoverFEMSolution(X, b, x);
358  const bool by_vdim = opt.by_vdim;
359  GridFunction *nodes = by_vdim ? &x0 : S.fes->GetMesh()->GetNodes();
362  double rnorm = nodes->DistanceTo(x) / nodes->Norml2();
363  double glob_norm;
364  MPI_Allreduce(&rnorm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
365  rnorm = glob_norm;
366  if (!opt.id) { mfem::out << "rnorm = " << rnorm << std::endl; }
367  const double lambda = opt.lambda;
368  if (by_vdim)
369  {
371  return Converged(rnorm);
372  }
374  {
375  GridFunction delta(S.fes);
376  subtract(x, *nodes, delta); // delta = x - nodes
377  // position and Δ vectors at point i
378  Vector ni(SDIM), di(SDIM);
379  for (int i = 0; i < delta.Size()/SDIM; i++)
380  {
381  // extract local vectors
382  const int ndof = S.fes->GetNDofs();
383  for (int d = 0; d < SDIM; d++)
384  {
385  ni(d) = (*nodes)(d*ndof + i);
386  di(d) = delta(d*ndof + i);
387  }
388  // project the delta vector in radial direction
389  const double ndotd = (ni*di) / (ni*ni);
390  di.Set(ndotd,ni);
391  // set global vectors
392  for (int d = 0; d < SDIM; d++) { delta(d*ndof + i) = di(d); }
393  }
395  }
396  // x = lambda*nodes + (1-lambda)*x
397  add(lambda, *nodes, (1.0-lambda), x, x);
398  return Converged(rnorm);
399  }
400
401  void Amr()
402  {
403  MFEM_VERIFY(opt.amr_threshold >= 0.0 && opt.amr_threshold <= 1.0, "");
404  Mesh *smesh = S.mesh;
405  Array<Refinement> amr;
406  const int NE = smesh->GetNE();
408  for (int e = 0; e < NE; e++)
409  {
410  double minW = +NL_DMAX;
411  double maxW = -NL_DMAX;
413  const Geometry::Type &type =
414  smesh->GetElement(e)->GetGeometryType();
415
416  const IntegrationRule *ir = &IntRules.Get(type, opt.order);
417  const int NQ = ir->GetNPoints();
418  for (int q = 0; q < NQ; q++)
419  {
420  eTr->SetIntPoint(&ir->IntPoint(q));
421  const DenseMatrix &J = eTr->Jacobian();
425  const double w = Jadjt.Weight();
426  minW = std::fmin(minW, w);
427  maxW = std::fmax(maxW, w);
428  }
429  if (std::fabs(maxW) != 0.0)
430  {
431  const double rho = minW / maxW;
432  MFEM_VERIFY(rho <= 1.0, "");
433  if (rho < opt.amr_threshold) { amr.Append(Refinement(e)); }
434  }
435  }
436  if (amr.Size()>0)
437  {
439  smesh->GeneralRefinement(amr);
440  S.fes->Update();
442  x.Update();
443  a.Update();
445  b.Update();
446  S.BoundaryConditions();
447  }
448  }
449  };
450
451  // Surface solver 'by vector'
452  class ByNodes: public Solver
453  {
454  public:
455  ByNodes(Surface &S, Opt &opt): Solver(S, opt)
457
458  bool Step()
459  {
460  x = *S.fes->GetMesh()->GetNodes();
461  bool converge = ParAXeqB();
462  S.mesh->SetNodes(x);
463  return converge ? true : false;
464  }
465  };
466
467  // Surface solver 'by ByVDim'
468  class ByVDim: public Solver
469  {
470  public:
471  void SetNodes(const GridFunction &Xi, const int c)
472  {
474  auto d_nodes = S.fes->GetMesh()->GetNodes()->Write();
475  const int ndof = S.fes->GetNDofs();
476  mfem::forall(ndof, [=] MFEM_HOST_DEVICE (int i)
477  {
478  d_nodes[c*ndof + i] = d_Xi[i];
479  });
480  }
481
482  void GetNodes(GridFunction &Xi, const int c)
483  {
484  auto d_Xi = Xi.Write();
485  const int ndof = S.fes->GetNDofs();
487  mfem::forall(ndof, [=] MFEM_HOST_DEVICE (int i)
488  {
489  d_Xi[i] = d_nodes[c*ndof + i];
490  });
491  }
492
493  ByVDim(Surface &S, Opt &opt): Solver(S, opt)
495
496  bool Step()
497  {
498  bool cvg[SDIM] {false};
499  for (int c=0; c < SDIM; ++c)
500  {
501  GetNodes(x, c);
502  x0 = x;
503  cvg[c] = ParAXeqB();
504  SetNodes(x, c);
505  }
506  const bool converged = cvg[0] && cvg[1] && cvg[2];
507  return converged ? true : false;
508  }
509  };
510 };
511
512 // #0: Default surface mesh file
513 struct MeshFromFile: public Surface
514 {
515  MeshFromFile(Opt &opt): Surface(opt, opt.mesh_file) { }
516 };
517
518 // #1: Catenoid surface
519 struct Catenoid: public Surface
520 {
521  Catenoid(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
522
523  void Prefix()
524  {
525  SetCurvature(opt.order, false, SDIM, Ordering::byNODES);
526  Array<int> v2v(GetNV());
527  for (int i = 0; i < v2v.Size(); i++) { v2v[i] = i; }
528  // identify vertices on vertical lines
529  for (int j = 0; j <= opt.ny; j++)
530  {
531  const int v_old = opt.nx + j * (opt.nx + 1);
532  const int v_new = j * (opt.nx + 1);
533  v2v[v_old] = v_new;
534  }
535  // renumber elements
536  for (int i = 0; i < GetNE(); i++)
537  {
538  Element *el = GetElement(i);
539  int *v = el->GetVertices();
540  const int nv = el->GetNVertices();
541  for (int j = 0; j < nv; j++)
542  {
543  v[j] = v2v[v[j]];
544  }
545  }
546  // renumber boundary elements
547  for (int i = 0; i < GetNBE(); i++)
548  {
549  Element *el = GetBdrElement(i);
550  int *v = el->GetVertices();
551  const int nv = el->GetNVertices();
552  for (int j = 0; j < nv; j++)
553  {
554  v[j] = v2v[v[j]];
555  }
556  }
557  RemoveUnusedVertices();
558  RemoveInternalBoundaries();
559  }
560
561  static void Parametrization(const Vector &x, Vector &p)
562  {
563  p.SetSize(SDIM);
564  // u in [0,2π] and v in [-π/6,π/6]
565  const double u = 2.0*PI*x[0];
566  const double v = PI*(x[1]-0.5)/3.;
567  p[0] = cos(u);
568  p[1] = sin(u);
569  p[2] = v;
570  }
571 };
572
573 // #2: Helicoid surface
574 struct Helicoid: public Surface
575 {
576  Helicoid(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
577
578  static void Parametrization(const Vector &x, Vector &p)
579  {
580  p.SetSize(SDIM);
581  // u in [0,2π] and v in [-2π/3,2π/3]
582  const double u = 2.0*PI*x[0];
583  const double v = 2.0*PI*(2.0*x[1]-1.0)/3.0;
584  p(0) = sin(u)*v;
585  p(1) = cos(u)*v;
586  p(2) = u;
587  }
588 };
589
590 // #3: Enneper's surface
591 struct Enneper: public Surface
592 {
593  Enneper(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
594
595  static void Parametrization(const Vector &x, Vector &p)
596  {
597  p.SetSize(SDIM);
598  // (u,v) in [-2, +2]
599  const double u = 4.0*(x[0]-0.5);
600  const double v = 4.0*(x[1]-0.5);
601  p[0] = +u - u*u*u/3.0 + u*v*v;
602  p[1] = -v - u*u*v + v*v*v/3.0;
603  p[2] = u*u - v*v;
604  }
605 };
606
607 // #4: Hold surface
608 struct Hold: public Surface
609 {
610  Hold(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
611
612  void Prefix()
613  {
614  SetCurvature(opt.order, false, SDIM, Ordering::byNODES);
615  Array<int> v2v(GetNV());
616  for (int i = 0; i < v2v.Size(); i++) { v2v[i] = i; }
617  // identify vertices on vertical lines
618  for (int j = 0; j <= opt.ny; j++)
619  {
620  const int v_old = opt.nx + j * (opt.nx + 1);
621  const int v_new = j * (opt.nx + 1);
622  v2v[v_old] = v_new;
623  }
624  // renumber elements
625  for (int i = 0; i < GetNE(); i++)
626  {
627  Element *el = GetElement(i);
628  int *v = el->GetVertices();
629  const int nv = el->GetNVertices();
630  for (int j = 0; j < nv; j++)
631  {
632  v[j] = v2v[v[j]];
633  }
634  }
635  // renumber boundary elements
636  for (int i = 0; i < GetNBE(); i++)
637  {
638  Element *el = GetBdrElement(i);
639  int *v = el->GetVertices();
640  const int nv = el->GetNVertices();
641  for (int j = 0; j < nv; j++)
642  {
643  v[j] = v2v[v[j]];
644  }
645  }
646  RemoveUnusedVertices();
647  RemoveInternalBoundaries();
648  }
649
650  static void Parametrization(const Vector &x, Vector &p)
651  {
652  p.SetSize(SDIM);
653  // u in [0,2π] and v in [0,1]
654  const double u = 2.0*PI*x[0];
655  const double v = x[1];
656  p[0] = cos(u)*(1.0 + 0.3*sin(3.*u + PI*v));
657  p[1] = sin(u)*(1.0 + 0.3*sin(3.*u + PI*v));
658  p[2] = v;
659  }
660 };
661
662 // #5: Costa minimal surface
663 #include <complex>
664 using cdouble = std::complex<double>;
665 #define I cdouble(0.0, 1.0)
666
667 // https://dlmf.nist.gov/20.2
668 cdouble EllipticTheta(const int a, const cdouble u, const cdouble q)
669 {
670  cdouble J = 0.0;
671  double delta = std::numeric_limits<double>::max();
672  switch (a)
673  {
674  case 1:
675  for (int n=0; delta > EPS; n+=1)
676  {
677  const cdouble j = pow(-1,n)*pow(q,n*(n+1.0))*sin((2.0*n+1.0)*u);
678  delta = abs(j);
679  J += j;
680  }
681  return 2.0*pow(q,0.25)*J;
682
683  case 2:
684  for (int n=0; delta > EPS; n+=1)
685  {
686  const cdouble j = pow(q,n*(n+1))*cos((2.0*n+1.0)*u);
687  delta = abs(j);
688  J += j;
689  }
690  return 2.0*pow(q,0.25)*J;
691  case 3:
692  for (int n=1; delta > EPS; n+=1)
693  {
694  const cdouble j = pow(q,n*n)*cos(2.0*n*u);
695  delta = abs(j);
696  J += j;
697  }
698  return 1.0 + 2.0*J;
699  case 4:
700  for (int n=1; delta > EPS; n+=1)
701  {
702  const cdouble j =pow(-1,n)*pow(q,n*n)*cos(2.0*n*u);
703  delta = abs(j);
704  J += j;
705  }
706  return 1.0 + 2.0*J;
707  }
708  return J;
709 }
710
711 // https://dlmf.nist.gov/23.6#E5
713  const cdouble w1 = 0.5,
714  const cdouble w3 = 0.5*I)
715 {
716  const cdouble tau = w3/w1;
717  const cdouble q = exp(I*M_PI*tau);
718  const cdouble e1 = M_PI*M_PI/(12.0*w1*w1)*
719  (1.0*pow(EllipticTheta(2,0,q),4) +
720  2.0*pow(EllipticTheta(4,0,q),4));
721  const cdouble u = M_PI*z / (2.0*w1);
722  const cdouble P = M_PI * EllipticTheta(3,0,q)*EllipticTheta(4,0,q) *
723  EllipticTheta(2,u,q)/(2.0*w1*EllipticTheta(1,u,q));
724  return P*P + e1;
725 }
726
727 cdouble EllipticTheta1Prime(const int k, const cdouble u, const cdouble q)
728 {
729  cdouble J = 0.0;
730  double delta = std::numeric_limits<double>::max();
731  for (int n=0; delta > EPS; n+=1)
732  {
733  const double alpha = 2.0*n+1.0;
734  const cdouble Dcosine = pow(alpha,k)*sin(k*M_PI/2.0 + alpha*u);
735  const cdouble j = pow(-1,n)*pow(q,n*(n+1.0))*Dcosine;
736  delta = abs(j);
737  J += j;
738  }
739  return 2.0*pow(q,0.25)*J;
740 }
741
742 // Logarithmic Derivative of Theta Function 1
744 {
745  cdouble J = 0.0;
746  double delta = std::numeric_limits<double>::max();
747  for (int n=1; delta > EPS; n+=1)
748  {
749  cdouble q2n = pow(q, 2*n);
750  if (abs(q2n) < EPS) { q2n = 0.0; }
751  const cdouble j = q2n/(1.0-q2n)*sin(2.0*n*u);
752  delta = abs(j);
753  J += j;
754  }
755  return 1.0/tan(u) + 4.0*J;
756 }
757
758 // https://dlmf.nist.gov/23.6#E13
760  const cdouble w1 = 0.5,
761  const cdouble w3 = 0.5*I)
762 {
763  const cdouble tau = w3/w1;
764  const cdouble q = exp(I*M_PI*tau);
765  const cdouble n1 = -M_PI*M_PI/(12.0*w1) *
766  (EllipticTheta1Prime(3,0,q)/
767  EllipticTheta1Prime(1,0,q));
768  const cdouble u = M_PI*z / (2.0*w1);
769  return z*n1/w1 + M_PI/(2.0*w1)*LogEllipticTheta1Prime(u,q);
770 }
771
772 // https://www.mathcurve.com/surfaces.gb/costa/costa.shtml
773 static double ALPHA[4] {0.0};
774 struct Costa: public Surface
775 {
776  Costa(Opt &opt): Surface((opt.Tptr = Parametrization, opt), false) { }
777
778  void Prefix()
779  {
780  ALPHA[3] = opt.tau;
781  const int nx = opt.nx, ny = opt.ny;
782  MFEM_VERIFY(nx>2 && ny>2, "");
783  const int nXhalf = (nx%2)==0 ? 4 : 2;
784  const int nYhalf = (ny%2)==0 ? 4 : 2;
785  const int nxh = nXhalf + nYhalf;
786  const int NVert = (nx+1) * (ny+1);
787  const int NElem = nx*ny - 4 - nxh;
788  const int NBdrElem = 0;
789  InitMesh(DIM, SDIM, NVert, NElem, NBdrElem);
790  // Sets vertices and the corresponding coordinates
791  for (int j = 0; j <= ny; j++)
792  {
793  const double cy = ((double) j / ny) ;
794  for (int i = 0; i <= nx; i++)
795  {
796  const double cx = ((double) i / nx);
797  const double coords[SDIM] = {cx, cy, 0.0};
799  }
800  }
801  // Sets elements and the corresponding indices of vertices
802  for (int j = 0; j < ny; j++)
803  {
804  for (int i = 0; i < nx; i++)
805  {
806  if (i == 0 && j == 0) { continue; }
807  if (i+1 == nx && j == 0) { continue; }
808  if (i == 0 && j+1 == ny) { continue; }
809  if (i+1 == nx && j+1 == ny) { continue; }
810  if ((j == 0 || j+1 == ny) && (abs(nx-(i<<1)-1)<=1)) { continue; }
811  if ((i == 0 || i+1 == nx) && (abs(ny-(j<<1)-1)<=1)) { continue; }
812  const int i0 = i + j*(nx+1);
813  const int i1 = i+1 + j*(nx+1);
814  const int i2 = i+1 + (j+1)*(nx+1);
815  const int i3 = i + (j+1)*(nx+1);
816  const int ind[4] = {i0, i1, i2, i3};
818  }
819  }
820  RemoveUnusedVertices();
822  FinalizeTopology();
823  SetCurvature(opt.order, false, SDIM, Ordering::byNODES);
824  }
825
826  static void Parametrization(const Vector &X, Vector &p)
827  {
828  const double tau = ALPHA[3];
829  Vector x = X;
830  x -= +0.5;
831  x *= tau;
832  x -= -0.5;
833
834  p.SetSize(3);
835  const bool y_top = x[1] > 0.5;
836  const bool x_top = x[0] > 0.5;
837  double u = x[0];
838  double v = x[1];
839  if (y_top) { v = 1.0 - x[1]; }
840  if (x_top) { u = 1.0 - x[0]; }
841  const cdouble w = u + I*v;
842  const cdouble w3 = I/2.;
843  const cdouble w1 = 1./2.;
844  const cdouble pw = WeierstrassP(w);
845  const cdouble e1 = WeierstrassP(0.5);
846  const cdouble zw = WeierstrassZeta(w);
847  const cdouble dw = WeierstrassZeta(w-w1) - WeierstrassZeta(w-w3);
848  p[0] = real(PI*(u+PI/(4.*e1))- zw +PI/(2.*e1)*(dw));
849  p[1] = real(PI*(v+PI/(4.*e1))-I*zw-PI*I/(2.*e1)*(dw));
850  p[2] = sqrt(PI/2.)*log(abs((pw-e1)/(pw+e1)));
851  if (y_top) { p[1] *= -1.0; }
852  if (x_top) { p[0] *= -1.0; }
853  const bool nan = std::isnan(p[0]) || std::isnan(p[1]) || std::isnan(p[2]);
854  MFEM_VERIFY(!nan, "nan");
855  ALPHA[0] = std::fmax(p[0], ALPHA[0]);
856  ALPHA[1] = std::fmax(p[1], ALPHA[1]);
857  ALPHA[2] = std::fmax(p[2], ALPHA[2]);
858  }
859
860  void Snap()
861  {
862  Vector node(SDIM);
863  MFEM_VERIFY(ALPHA[0] > 0.0,"");
864  MFEM_VERIFY(ALPHA[1] > 0.0,"");
865  MFEM_VERIFY(ALPHA[2] > 0.0,"");
866  GridFunction &nodes = *GetNodes();
867  const double phi = (1.0 + sqrt(5.0)) / 2.0;
868  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
869  {
870  for (int d = 0; d < SDIM; d++)
871  {
872  const double alpha = d==2 ? phi : 1.0;
873  const int vdof = nodes.FESpace()->DofToVDof(i, d);
874  nodes(vdof) /= alpha * ALPHA[d];
875  }
876  }
877  }
878 };
879
880 // #6: Shell surface model
881 struct Shell: public Surface
882 {
883  Shell(Opt &opt): Surface((opt.niters = 1, opt.Tptr = Parametrization, opt)) { }
884
885  static void Parametrization(const Vector &x, Vector &p)
886  {
887  p.SetSize(3);
888  // u in [0,2π] and v in [-15, 6]
889  const double u = 2.0*PI*x[0];
890  const double v = 21.0*x[1]-15.0;
891  p[0] = +1.0*pow(1.16,v)*cos(v)*(1.0+cos(u));
892  p[1] = -1.0*pow(1.16,v)*sin(v)*(1.0+cos(u));
893  p[2] = -2.0*pow(1.16,v)*(1.0+sin(u));
894  }
895 };
896
897 // #7: Scherk's doubly periodic surface
898 struct Scherk: public Surface
899 {
900  static void Parametrization(const Vector &x, Vector &p)
901  {
902  p.SetSize(SDIM);
903  const double alpha = 0.49;
904  // (u,v) in [-απ, +απ]
905  const double u = alpha*PI*(2.0*x[0]-1.0);
906  const double v = alpha*PI*(2.0*x[1]-1.0);
907  p[0] = u;
908  p[1] = v;
909  p[2] = log(cos(v)/cos(u));
910  }
911
912  Scherk(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
913 };
914
915 // #8: Full Peach street model
916 struct FullPeach: public Surface
917 {
918  static constexpr int NV = 8;
919  static constexpr int NE = 6;
920
921  FullPeach(Opt &opt):
922  Surface((opt.niters = std::min(4, opt.niters), opt), NV, NE, 0) { }
923
924  void Prefix()
925  {
927  {
928  {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
929  {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
930  };
932  {
933  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
934  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
935
936  };
937  for (int j = 0; j < NV; j++) { AddVertex(quad_v[j]); }
939
941  FinalizeTopology(false);
942  UniformRefinement();
943  SetCurvature(opt.order, false, SDIM, Ordering::byNODES);
944  }
945
946  void Snap() { SnapNodesToUnitSphere(); }
947
948  void BoundaryConditions()
949  {
950  Vector X(SDIM);
951  Array<int> dofs;
952  Array<int> ess_vdofs, ess_tdofs;
953  ess_vdofs.SetSize(fes->GetVSize());
954  MFEM_VERIFY(fes->GetVSize() >= fes->GetTrueVSize(),"");
955  ess_vdofs = 0;
956  DenseMatrix PointMat;
958  for (int e = 0; e < fes->GetNE(); e++)
959  {
960  fes->GetElementDofs(e, dofs);
961  const IntegrationRule &ir = fes->GetFE(e)->GetNodes();
963  eTr->Transform(ir, PointMat);
964  Vector one(dofs.Size());
965  for (int dof = 0; dof < dofs.Size(); dof++)
966  {
967  one = 0.0;
968  one[dof] = 1.0;
969  const int k = dofs[dof];
970  MFEM_ASSERT(k >= 0, "");
971  PointMat.Mult(one, X);
972  const bool halfX = std::fabs(X[0]) < EPS && X[1] <= 0.0;
973  const bool halfY = std::fabs(X[2]) < EPS && X[1] >= 0.0;
974  const bool is_on_bc = halfX || halfY;
975  for (int c = 0; c < SDIM; c++)
976  { ess_vdofs[fes->DofToVDof(k, c)] = is_on_bc; }
977  }
978  }
979  const SparseMatrix *R = fes->GetRestrictionMatrix();
980  if (!R)
981  {
982  ess_tdofs.MakeRef(ess_vdofs);
983  }
984  else
985  {
986  R->BooleanMult(ess_vdofs, ess_tdofs);
987  }
990  }
991 };
992
993 // #9: 1/4th Peach street model
994 struct QuarterPeach: public Surface
995 {
996  QuarterPeach(Opt &opt): Surface((opt.Tptr = Parametrization, opt)) { }
997
998  static void Parametrization(const Vector &X, Vector &p)
999  {
1000  p = X;
1001  const double x = 2.0*X[0]-1.0;
1002  const double y = X[1];
1003  const double r = sqrt(x*x + y*y);
1004  const double t = (x==0.0) ? PI/2.0 :
1005  (y==0.0 && x>0.0) ? 0. :
1006  (y==0.0 && x<0.0) ? PI : acos(x/r);
1007  const double sqrtx = sqrt(1.0 + x*x);
1008  const double sqrty = sqrt(1.0 + y*y);
1009  const bool yaxis = PI/4.0<t && t < 3.0*PI/4.0;
1010  const double R = yaxis?sqrtx:sqrty;
1011  const double gamma = r/R;
1012  p[0] = gamma * cos(t);
1013  p[1] = gamma * sin(t);
1014  p[2] = 1.0 - gamma;
1015  }
1016
1017  void Postfix()
1018  {
1019  for (int i = 0; i < GetNBE(); i++)
1020  {
1021  Element *el = GetBdrElement(i);
1022  const int fn = GetBdrElementEdgeIndex(i);
1023  MFEM_VERIFY(!FaceIsTrueInterior(fn),"");
1024  Array<int> vertices;
1025  GetFaceVertices(fn, vertices);
1026  const GridFunction *nodes = GetNodes();
1027  Vector nval;
1028  double R[2], X[2][SDIM];
1029  for (int v = 0; v < 2; v++)
1030  {
1031  R[v] = 0.0;
1032  const int iv = vertices[v];
1033  for (int d = 0; d < SDIM; d++)
1034  {
1035  nodes->GetNodalValues(nval, d+1);
1036  const double x = X[v][d] = nval[iv];
1037  if (d < 2) { R[v] += x*x; }
1038  }
1039  }
1040  if (std::fabs(X[0][1])<=EPS && std::fabs(X[1][1])<=EPS &&
1041  (R[0]>0.1 || R[1]>0.1))
1042  { el->SetAttribute(1); }
1043  else { el->SetAttribute(2); }
1044  }
1045  }
1046 };
1047
1048 // #10: Slotted sphere mesh
1049 struct SlottedSphere: public Surface
1050 {
1051  SlottedSphere(Opt &opt): Surface((opt.niters = 4, opt), 64, 40, 0) { }
1052
1053  void Prefix()
1054  {
1055  constexpr double delta = 0.15;
1056  constexpr int NV1D = 4;
1057  constexpr int NV = NV1D*NV1D*NV1D;
1058  constexpr int NEPF = (NV1D-1)*(NV1D-1);
1059  constexpr int NE = NEPF*6;
1060  const double V1D[NV1D] = {-1.0, -delta, delta, 1.0};
1061  double QV[NV][3];
1062  for (int iv=0; iv<NV; ++iv)
1063  {
1064  int ix = iv % NV1D;
1065  int iy = (iv / NV1D) % NV1D;
1066  int iz = (iv / NV1D) / NV1D;
1067
1068  QV[iv][0] = V1D[ix];
1069  QV[iv][1] = V1D[iy];
1070  QV[iv][2] = V1D[iz];
1071  }
1072  int QE[NE][4];
1073  for (int ix=0; ix<NV1D-1; ++ix)
1074  {
1075  for (int iy=0; iy<NV1D-1; ++iy)
1076  {
1077  int el_offset = ix + iy*(NV1D-1);
1078  // x = 0
1079  QE[0*NEPF + el_offset][0] = NV1D*ix + NV1D*NV1D*iy;
1080  QE[0*NEPF + el_offset][1] = NV1D*(ix+1) + NV1D*NV1D*iy;
1081  QE[0*NEPF + el_offset][2] = NV1D*(ix+1) + NV1D*NV1D*(iy+1);
1082  QE[0*NEPF + el_offset][3] = NV1D*ix + NV1D*NV1D*(iy+1);
1083  // x = 1
1084  int x_off = NV1D-1;
1085  QE[1*NEPF + el_offset][3] = x_off + NV1D*ix + NV1D*NV1D*iy;
1086  QE[1*NEPF + el_offset][2] = x_off + NV1D*(ix+1) + NV1D*NV1D*iy;
1087  QE[1*NEPF + el_offset][1] = x_off + NV1D*(ix+1) + NV1D*NV1D*(iy+1);
1088  QE[1*NEPF + el_offset][0] = x_off + NV1D*ix + NV1D*NV1D*(iy+1);
1089  // y = 0
1090  QE[2*NEPF + el_offset][0] = NV1D*NV1D*iy + ix;
1091  QE[2*NEPF + el_offset][1] = NV1D*NV1D*iy + ix + 1;
1092  QE[2*NEPF + el_offset][2] = NV1D*NV1D*(iy+1) + ix + 1;
1093  QE[2*NEPF + el_offset][3] = NV1D*NV1D*(iy+1) + ix;
1094  // y = 1
1095  int y_off = NV1D*(NV1D-1);
1096  QE[3*NEPF + el_offset][0] = y_off + NV1D*NV1D*iy + ix;
1097  QE[3*NEPF + el_offset][1] = y_off + NV1D*NV1D*iy + ix + 1;
1098  QE[3*NEPF + el_offset][2] = y_off + NV1D*NV1D*(iy+1) + ix + 1;
1099  QE[3*NEPF + el_offset][3] = y_off + NV1D*NV1D*(iy+1) + ix;
1100  // z = 0
1101  QE[4*NEPF + el_offset][0] = NV1D*iy + ix;
1102  QE[4*NEPF + el_offset][1] = NV1D*iy + ix + 1;
1103  QE[4*NEPF + el_offset][2] = NV1D*(iy+1) + ix + 1;
1104  QE[4*NEPF + el_offset][3] = NV1D*(iy+1) + ix;
1105  // z = 1
1106  int z_off = NV1D*NV1D*(NV1D-1);
1107  QE[5*NEPF + el_offset][0] = z_off + NV1D*iy + ix;
1108  QE[5*NEPF + el_offset][1] = z_off + NV1D*iy + ix + 1;
1109  QE[5*NEPF + el_offset][2] = z_off + NV1D*(iy+1) + ix + 1;
1110  QE[5*NEPF + el_offset][3] = z_off + NV1D*(iy+1) + ix;
1111  }
1112  }
1113  // Delete on x = 0 face
1114  QE[0*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1115  QE[0*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1116  // Delete on x = 1 face
1117  QE[1*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1118  QE[1*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1119  // Delete on y = 1 face
1120  QE[3*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1121  QE[3*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1122  // Delete on z = 1 face
1123  QE[5*NEPF + 0 + 1*(NV1D-1)][0] = -1;
1124  QE[5*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1125  QE[5*NEPF + 2 + 1*(NV1D-1)][0] = -1;
1126  // Delete on z = 0 face
1127  QE[4*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1128  QE[4*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1129  QE[4*NEPF + 1 + 2*(NV1D-1)][0] = -1;
1130  // Delete on y = 0 face
1131  QE[2*NEPF + 1 + 0*(NV1D-1)][0] = -1;
1132  QE[2*NEPF + 1 + 1*(NV1D-1)][0] = -1;
1133  for (int j = 0; j < NV; j++) { AddVertex(QV[j]); }
1134  for (int j = 0; j < NE; j++)
1135  {
1136  if (QE[j][0] < 0) { continue; }
1138  }
1139  RemoveUnusedVertices();
1141  EnsureNodes();
1142  FinalizeTopology();
1143  }
1144
1145  void Snap() { SnapNodesToUnitSphere(); }
1146 };
1147
1148 static int Problem0(Opt &opt)
1149 {
1150  // Create our surface mesh from command line options
1151  Surface *S = nullptr;
1152  switch (opt.surface)
1153  {
1154  case 0: S = new MeshFromFile(opt); break;
1155  case 1: S = new Catenoid(opt); break;
1156  case 2: S = new Helicoid(opt); break;
1157  case 3: S = new Enneper(opt); break;
1158  case 4: S = new Hold(opt); break;
1159  case 5: S = new Costa(opt); break;
1160  case 6: S = new Shell(opt); break;
1161  case 7: S = new Scherk(opt); break;
1162  case 8: S = new FullPeach(opt); break;
1163  case 9: S = new QuarterPeach(opt); break;
1164  case 10: S = new SlottedSphere(opt); break;
1165  default: MFEM_ABORT("Unknown surface (surface <= 10)!");
1166  }
1167  S->Create();
1168  S->Solve();
1169  delete S;
1170  return 0;
1171 }
1172
1173 // Problem 1: solve the Dirichlet problem for the minimal surface equation
1174 // of the form z=f(x,y), using Picard iterations.
1175 static double u0(const Vector &x) { return sin(3.0 * PI * (x[1] + x[0])); }
1176
1177 enum {NORM, AREA};
1178
1179 static double qf(const int order, const int ker, Mesh &m,
1181 {
1182  const Geometry::Type type = m.GetElementBaseGeometry(0);
1183  const IntegrationRule &ir(IntRules.Get(type, order));
1185
1186  const int NE(m.GetNE());
1187  const int ND(fes.GetFE(0)->GetDof());
1188  const int NQ(ir.GetNPoints());
1190  const GeometricFactors *geom = m.GetGeometricFactors(ir, flags);
1191
1192  const int D1D = fes.GetFE(0)->GetOrder() + 1;
1193  const int Q1D = IntRules.Get(Geometry::SEGMENT, ir.GetOrder()).GetNPoints();
1194  MFEM_VERIFY(ND == D1D*D1D, "");
1195  MFEM_VERIFY(NQ == Q1D*Q1D, "");
1196
1197  Vector Eu(ND*NE), grad_u(DIM*NQ*NE), sum(NE*NQ), one(NE*NQ);
1198  qi->SetOutputLayout(QVectorLayout::byVDIM);
1200  const Operator *G(fes.GetElementRestriction(e_ordering));
1201  G->Mult(u, Eu);
1203
1204  auto W = Reshape(ir.GetWeights().Read(), Q1D, Q1D);
1205  auto J = Reshape(geom->J.Read(), Q1D, Q1D, DIM, DIM, NE);
1206  auto detJ = Reshape(geom->detJ.Read(), Q1D, Q1D, NE);
1208  auto S = Reshape(sum.Write(), Q1D, Q1D, NE);
1209
1210  mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
1211  {
1213  {
1215  {
1216  const double w = W(qx, qy);
1217  const double J11 = J(qx, qy, 0, 0, e);
1218  const double J12 = J(qx, qy, 1, 0, e);
1219  const double J21 = J(qx, qy, 0, 1, e);
1220  const double J22 = J(qx, qy, 1, 1, e);
1221  const double det = detJ(qx, qy, e);
1222  const double area = w * det;
1223  const double gu0 = grdU(0, qx, qy, e);
1224  const double gu1 = grdU(1, qx, qy, e);
1225  const double tgu0 = (J22*gu0 - J12*gu1)/det;
1226  const double tgu1 = (J11*gu1 - J21*gu0)/det;
1227  const double ngu = tgu0*tgu0 + tgu1*tgu1;
1228  const double s = (ker == AREA) ? sqrt(1.0 + ngu) :
1229  (ker == NORM) ? ngu : 0.0;
1230  S(qx, qy, e) = area * s;
1231  }
1232  }
1233  });
1234  one = 1.0;
1235  return sum * one;
1236 }
1237
1238 static int Problem1(Opt &opt)
1239 {
1240  const int order = opt.order;
1241  Mesh smesh = Mesh::MakeCartesian2D(opt.nx, opt.ny, QUAD);
1242  smesh.SetCurvature(opt.order, false, DIM, Ordering::byNODES);
1243  for (int l = 0; l < opt.refine; l++) { smesh.UniformRefinement(); }
1244  ParMesh mesh(MPI_COMM_WORLD, smesh);
1245  const H1_FECollection fec(order, DIM);
1246  ParFiniteElementSpace fes(&mesh, &fec);
1247  Array<int> ess_tdof_list;
1248  if (mesh.bdr_attributes.Size())
1249  {
1250  Array<int> ess_bdr(mesh.bdr_attributes.Max());
1251  ess_bdr = 1;
1252  fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
1253  }
1254  ParGridFunction uold(&fes), u(&fes), b(&fes);
1255  FunctionCoefficient u0_fc(u0);
1256  u.ProjectCoefficient(u0_fc);
1257  socketstream glvis;
1258  if (opt.vis) { opt.vis = glvis.open(vishost, visport) == 0; }
1259  if (opt.vis) { Surface::Visualize(glvis, opt, &mesh, GLVIZ_W, GLVIZ_H, &u); }
1260  Vector B, X;
1261  OperatorPtr A;
1262  CGSolver cg(MPI_COMM_WORLD);
1263  cg.SetRelTol(EPS);
1264  cg.SetAbsTol(EPS*EPS);
1265  cg.SetMaxIter(400);
1266  cg.SetPrintLevel(0);
1267  ParGridFunction eps(&fes);
1268  for (int i = 0; i < opt.niters; i++)
1269  {
1270  b = 0.0;
1271  uold = u;
1272  ParBilinearForm a(&fes);
1273  if (opt.pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
1274  const double q_uold = qf(order, AREA, mesh, fes, uold);
1275  MFEM_VERIFY(std::fabs(q_uold) > EPS,"");
1276  ConstantCoefficient q_uold_cc(1.0/sqrt(q_uold));
1278  a.Assemble();
1279  a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
1280  cg.SetOperator(*A);
1281  cg.Mult(B, X);
1282  a.RecoverFEMSolution(X, b, u);
1283  subtract(u, uold, eps);
1284  const double norm = sqrt(std::fabs(qf(order, NORM, mesh, fes, eps)));
1285  const double area = qf(order, AREA, mesh, fes, u);
1286  if (!opt.id)
1287  {
1288  mfem::out << "Iteration " << i << ", norm: " << norm
1289  << ", area: " << area << std::endl;
1290  }
1291  if (opt.vis) { Surface::Visualize(glvis, opt, &mesh, &u); }
1292  if (opt.print) { Surface::Print(opt, &mesh, &u); }
1293  if (norm < NRM) { break; }
1294  }
1295  return 0;
1296 }
1297
1298 int main(int argc, char *argv[])
1299 {
1300  Opt opt;
1301  Mpi::Init(argc, argv);
1302  opt.id = Mpi::WorldRank();
1303  opt.sz = Mpi::WorldSize();
1304  Hypre::Init();
1305
1306  // Parse command-line options.
1307  OptionsParser args(argc, argv);
1308  args.AddOption(&opt.pb, "-p", "--problem", "Problem to solve.");
1309  args.AddOption(&opt.mesh_file, "-m", "--mesh", "Mesh file to use.");
1310  args.AddOption(&opt.wait, "-w", "--wait", "-no-w", "--no-wait",
1311  "Enable or disable a GLVis pause.");
1313  "Enable or disable radial constraints in solver.");
1315  "Number of elements in x-direction.");
1317  "Number of elements in y-direction.");
1318  args.AddOption(&opt.order, "-o", "--order", "Finite element order.");
1320  args.AddOption(&opt.niters, "-n", "--niter-max", "Max number of iterations");
1321  args.AddOption(&opt.surface, "-s", "--surface", "Choice of the surface.");
1323  "--no-partial-assembly", "Enable Partial Assembly.");
1324  args.AddOption(&opt.tau, "-t", "--tau", "Costa scale factor.");
1325  args.AddOption(&opt.lambda, "-l", "--lambda", "Lambda step toward solution.");
1326  args.AddOption(&opt.amr, "-a", "--amr", "-no-a", "--no-amr", "Enable AMR.");
1327  args.AddOption(&opt.amr_threshold, "-at", "--amr-threshold", "AMR threshold.");
1329  "Device configuration string, see Device::Configure().");
1330  args.AddOption(&opt.keys, "-k", "--keys", "GLVis configuration keys.");
1332  "--no-visualization", "Enable or disable visualization.");
1334  "-no-c", "--solve-bynodes",
1335  "Enable or disable the 'ByVdim' solver");
1336  args.AddOption(&opt.print, "-print", "--print", "-no-print", "--no-print",
1337  "Enable or disable result output (files in mfem format).");
1338  args.AddOption(&opt.snapshot, "-ss", "--snapshot", "-no-ss", "--no-snapshot",
1339  "Enable or disable GLVis snapshot.");
1340  args.Parse();
1341  if (!args.Good()) { args.PrintUsage(mfem::out); return 1; }
1342  MFEM_VERIFY(opt.lambda >= 0.0 && opt.lambda <= 1.0,"");
1343  if (!opt.id) { args.PrintOptions(mfem::out); }
1344
1345  // Initialize hardware devices
1346  Device device(opt.device_config);
1347  if (!opt.id) { device.Print(); }
1348
1349  if (opt.pb == 0) { Problem0(opt); }
1350
1351  if (opt.pb == 1) { Problem1(opt); }
1352
1353  return 0;
1354 }
const T * Read(bool on_dev=true) const
Definition: array.hpp:307
const char vishost[]
Definition: array.hpp:327
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:80
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
Definition: solvers.hpp:493
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
constexpr double NRM
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
void forall_2D(int N, int X, int Y, lambda &&body)
Definition: forall.hpp:751
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
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.
Definition: mesh.cpp:847
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
constexpr double PI
virtual const double * HostRead() const
Definition: vector.hpp:457
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
std::complex< double > cdouble
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
Definition: densemat.cpp:2478
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:12045
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:706
virtual const double * Read(bool on_dev=true) const
Definition: vector.hpp:453
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 ...
Definition: fespace.cpp:587
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.
Definition: device.cpp:279
const Element * GetElement(int i) const
Return pointer to the i&#39;th element object.
Definition: mesh.hpp:1143
cdouble LogEllipticTheta1Prime(const cdouble u, const cdouble q)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition: mesh.hpp:2183
STL namespace.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
double Weight() const
Definition: densemat.cpp:544
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:379
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition: fespace.cpp:621
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
int close()
Close the socketstream.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1302
constexpr int GLVIZ_W
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
Vector J
Jacobians of the element transformations at all quadrature points.
Definition: mesh.hpp:2224
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:86
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
Definition: vector.hpp:669
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
constexpr int DIM
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
Vector detJ
Determinants of the Jacobians at all quadrature points.
Definition: mesh.hpp:2230
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
int GetOrder() const
Returns the order of the integration rule.
Definition: intrules.hpp:246
int main(int argc, char *argv[])
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
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.
Definition: mesh.cpp:5635
double delta
Definition: lissajous.cpp:43
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
static void Init()
Singleton creation with Mpi::Init();.
const int visport
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:172
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i&#39;th element for dimension vdim.
Definition: gridfunc.cpp:395
virtual DofTransformation * GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:470
cdouble WeierstrassP(const cdouble z, const cdouble w1=0.5, const cdouble w3=0.5 *I)
void SetAbsTol(double atol)
Definition: solvers.hpp:200
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1419
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
cdouble EllipticTheta(const int a, 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...
Definition: globals.hpp:66
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:473
constexpr double EPS
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
constexpr int SDIM
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)
Definition: mesh.cpp:3773
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1026
constexpr double NL_DMAX
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
double a
Definition: lissajous.cpp:41
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition: fespace.cpp:1370
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:355
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
Lexicographic ordering for tensor-product FiniteElements.
virtual int GetNVertices() const =0
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
Class for parallel bilinear form.
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:835
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
cdouble EllipticTheta1Prime(const int k, const cdouble u, const cdouble q)
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:251
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
constexpr int GLVIZ_H
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
RefCoord s[3]
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4825
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Abstract operator.
Definition: operator.hpp:24
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:389
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9820
Class for parallel meshes.
Definition: pmesh.hpp:32
Abstract data type element.
Definition: element.hpp:28
cdouble WeierstrassZeta(const cdouble z, const cdouble w1=0.5, const cdouble w3=0.5 *I)
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327