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 // Lissajous Miniapp: Spinning optical illusion
14 // ---------------------------------------------
15 //
16 // This miniapp generates two different Lissajous curves in 3D which appear to
17 // spin vertically and/or horizontally, even though the net motion is the same.
18 // Based on the 2019 Illusion of the year "Dual Axis Illusion" by Frank Force,
19 // see http://illusionoftheyear.com/2019/12/dual-axis-illusion.
20 //
21 // Compile with: make lissajous
22 //
23 // Sample runs: lissajous
24 // lissajous -a 5 -b 4
25 // lissajous -a 4 -b 3 -delta -90
26 // lissajous -o 8 -nx 3 -ny 3
27 // lissajous -a 11 -b 10 -o 4
28
29 #include "mfem.hpp"
30 #include <fstream>
31 #include <iostream>
32
33 using namespace std;
34 using namespace mfem;
35
36 double u_function(const Vector &x);
37 void lissajous_trans_v(const Vector &x, Vector &p);
38 void lissajous_trans_h(const Vector &x, Vector &p);
39
40 // Default Lissajous curve parameters
41 double a = 3.0;
42 double b = 2.0;
43 double delta = 90;
44
main(int argc,char * argv[])45 int main(int argc, char *argv[])
46 {
47 int nx = 32;
48 int ny = 3;
49 int order = 2;
50 bool visualization = true;
51
52 OptionsParser args(argc, argv);
53 args.AddOption(&nx, "-nx", "--num-elements-x",
54 "Number of elements in x-direction.");
55 args.AddOption(&ny, "-ny", "--num-elements-y",
56 "Number of elements in y-direction.");
57 args.AddOption(&order, "-o", "--mesh-order",
58 "Order (polynomial degree) of the mesh elements.");
59 args.AddOption(&a, "-a", "--x-frequency",
60 "Frequency of the x-component.");
61 args.AddOption(&b, "-b", "--y-frequency",
62 "Frequency of the y-component.");
63 args.AddOption(&delta, "-delta", "--x-phase",
64 "Phase angle of the x-component.");
65 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
66 "--no-visualization",
67 "Enable or disable GLVis visualization.");
68 args.Parse();
69 if (!args.Good())
70 {
71 args.PrintUsage(cout);
72 return 1;
73 }
74 args.PrintOptions(cout);
75
76 delta *= M_PI / 180.0; // convert to radians
77
78 char vishost[] = "localhost";
79 int visport = 19916;
80 socketstream soutv, south;
81
82 {
83 Mesh mesh = Mesh::MakeCartesian2D(
84 nx, ny, Element::QUADRILATERAL, 1, 2*M_PI, 2*M_PI);
85 mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
86 mesh.Transform(lissajous_trans_v);
87
88 H1_FECollection fec(order, 3);
89 FiniteElementSpace fes(&mesh, &fec);
90 GridFunction u(&fes);
91 FunctionCoefficient ufc(u_function);
92 u.ProjectCoefficient(ufc);
93
94 if (visualization)
95 {
96 soutv.open(vishost, visport);
97 soutv << "solution\n" << mesh << u;
98 soutv << "keys 'ARRj" << std::string(90, '7') << "'\n";
99 soutv << "palette 17 zoom 1.65 subdivisions 32 0\n";
100 soutv << "window_title 'V' window_geometry 0 0 500 500\n";
101 soutv << flush;
102 }
103 }
104
105 {
106 Mesh mesh = Mesh::MakeCartesian2D(
107 nx, ny, Element::QUADRILATERAL, 1, 2*M_PI, 2*M_PI);
108 mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
109 mesh.Transform(lissajous_trans_h);
110
111 H1_FECollection fec(order, 3);
112 FiniteElementSpace fes(&mesh, &fec);
113 GridFunction u(&fes);
114 FunctionCoefficient ufc(u_function);
115 u.ProjectCoefficient(ufc);
116
117 if (visualization)
118 {
119 south.open(vishost, visport);
120 south << "solution\n" << mesh << u;
121 south << "keys 'ARRj'\n";
122 south << "palette 17 zoom 1.65 subdivisions 32 0\n";
123 south << "window_title 'H' window_geometry 500 0 500 500\n";
124 south << flush;
125 }
126
127 ofstream mesh_ofs("lissajous.mesh");
128 mesh_ofs.precision(8);
129 mesh.Print(mesh_ofs);
130 ofstream sol_ofs("lissajous.gf");
131 sol_ofs.precision(8);
132 u.Save(sol_ofs);
133 }
134
135 soutv << "keys '.0" << std::string(b, '0') << "'\n" << flush;
136 south << "keys '.0" << std::string(a, '0') << "'\n" << flush;
137
138 cout << "Which direction(s) are the two curves spinning in?\n";
139
140 return 0;
141 }
142
143 // Simple function to project to help identify the spinning
u_function(const Vector & x)144 double u_function(const Vector &x)
145 {
146 return x[2];
147 }
148
149 // Tubular Lissajous curve with the given parameters (a, b, theta)
lissajous_trans(const Vector & x,Vector & p,double a,double b,double delta)150 void lissajous_trans(const Vector &x, Vector &p,
151 double a, double b, double delta)
152 {
153 p.SetSize(3);
154
155 double phi = x[0];
156 double theta = x[1];
157 double t = phi;
158
159 double A = b; // Scaling of the curve along the x-axis
160 double B = a; // Scaling of the curve along the y-axis
161
162 // Lissajous curve on a 3D cylinder
163 p[0] = B*cos(b*t);
164 p[1] = B*sin(b*t); // Y
165 p[2] = A*sin(a*t + delta); // X
166
167 // Turn the curve into a tubular surface
168 {
169 // tubular radius
170 double R = 0.02*(A+B);
171
172 // normal to the cylinder at p(t)
173 double normal[3] = { cos(b*t), sin(b*t), 0 };
174
175 // tangent to the curve, dp/dt(t)
176 // double tangent[3] = { -b*B*sin(b*t), b*B*cos(b*t), A*a*cos(a*t+delta) };
177
178 // normalized cross product of tangent and normal at p(t)
179 double cn = 1e-128;
180 double cross[3] = { A*a*sin(b*t)*cos(a*t+delta), -A*a*cos(b*t)*cos(a*t+delta), b*B };
181 for (int i = 0; i < 3; i++) { cn += cross[i]*cross[i]; }
182 for (int i = 0; i < 3; i++) { cross[i] /= sqrt(cn); }
183
184 // create a tubular surface of radius R around the curve p(t), in the plane
185 // orthogonal to the tangent (with basis given by normal and cross)
186 for (int i = 0; i < 3; i++)
187 {
188 p[i] += R * (cos(theta)*normal[i] + sin(theta)*cross[i]);
189 }
190 }
191 }
192
193 // Vertically spinning curve
lissajous_trans_v(const Vector & x,Vector & p)194 void lissajous_trans_v(const Vector &x, Vector &p)
195 {
196 return lissajous_trans(x, p, a, b, delta);
197 }
198
199 // Horizontally spinning curve
lissajous_trans_h(const Vector & x,Vector & p)200 void lissajous_trans_h(const Vector &x, Vector &p)
201 {
202 return lissajous_trans(x, p, b, a, delta);
203 }
204