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