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