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