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