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