MFEM  v4.1.0 Finite element discretization library
lissajous.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // ---------------------------------------------
13 // 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
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);
54  "Number of elements in x-direction.");
56  "Number of elements in y-direction.");
58  "Order (polynomial degree) of the mesh elements.");
60  "Frequency of the x-component.");
62  "Frequency of the y-component.");
64  "Phase angle of the x-component.");
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(nx, ny, Element::QUADRILATERAL, 1, 2*M_PI, 2*M_PI);
84  mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
86
87  H1_FECollection fec(order, 3);
88  FiniteElementSpace fes(&mesh, &fec);
89  GridFunction u(&fes);
91  u.ProjectCoefficient(ufc);
92
93  if (visualization)
94  {
95  soutv.open(vishost, visport);
96  soutv << "solution\n" << mesh << u;
97  soutv << "keys 'ARRj" << std::string(90, '7') << "'\n";
98  soutv << "palette 17 zoom 1.65 subdivisions 32 0\n";
99  soutv << "window_title 'V' window_geometry 0 0 500 500\n";
100  soutv << flush;
101  }
102  }
103
104  {
105  Mesh mesh(nx, ny, Element::QUADRILATERAL, 1, 2*M_PI, 2*M_PI);
106  mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
108
109  H1_FECollection fec(order, 3);
110  FiniteElementSpace fes(&mesh, &fec);
111  GridFunction u(&fes);
113  u.ProjectCoefficient(ufc);
114
115  if (visualization)
116  {
117  south.open(vishost, visport);
118  south << "solution\n" << mesh << u;
119  south << "keys 'ARRj'\n";
120  south << "palette 17 zoom 1.65 subdivisions 32 0\n";
121  south << "window_title 'H' window_geometry 500 0 500 500\n";
122  south << flush;
123  }
124
125  ofstream mesh_ofs("lissajous.mesh");
126  mesh_ofs.precision(8);
127  mesh.Print(mesh_ofs);
128  ofstream sol_ofs("lissajous.gf");
129  sol_ofs.precision(8);
130  u.Save(sol_ofs);
131  }
132
133  soutv << "keys '.0" << std::string(b, '0') << "'\n" << flush;
134  south << "keys '.0" << std::string(a, '0') << "'\n" << flush;
135
136  cout << "Which direction(s) are the two curves spinning in?\n";
137
138  return 0;
139 }
140
141 // Simple function to project to help identify the spinning
142 double u_function(const Vector &x)
143 {
144  return x[2];
145 }
146
147 // Tubular Lissajous curve with the given parameters (a, b, theta)
148 void lissajous_trans(const Vector &x, Vector &p,
149  double a, double b, double delta)
150 {
151  p.SetSize(3);
152
153  double phi = x[0];
154  double theta = x[1];
155  double t = phi;
156
157  double A = b; // Scaling of the curve along the x-axis
158  double B = a; // Scaling of the curve along the y-axis
159
160  // Lissajous curve on a 3D cylinder
161  p[0] = B*cos(b*t);
162  p[1] = B*sin(b*t); // Y
163  p[2] = A*sin(a*t + delta); // X
164
165  // Turn the curve into a tubular surface
166  {
168  double R = 0.02*(A+B);
169
170  // normal to the cylinder at p(t)
171  double normal[3] = { cos(b*t), sin(b*t), 0 };
172
173  // tangent to the curve, dp/dt(t)
174  // double tangent[3] = { -b*B*sin(b*t), b*B*cos(b*t), A*a*cos(a*t+delta) };
175
176  // normalized cross product of tangent and normal at p(t)
177  double cn = 1e-128;
178  double cross[3] = { A*a*sin(b*t)*cos(a*t+delta), -A*a*cos(b*t)*cos(a*t+delta), b*B };
179  for (int i = 0; i < 3; i++) { cn += cross[i]*cross[i]; }
180  for (int i = 0; i < 3; i++) { cross[i] /= sqrt(cn); }
181
182  // create a tubular surface of radius R around the curve p(t), in the plane
183  // orthogonal to the tangent (with basis given by normal and cross)
184  for (int i = 0; i < 3; i++)
185  {
186  p[i] += R * (cos(theta)*normal[i] + sin(theta)*cross[i]);
187  }
188  }
189 }
190
191 // Vertically spinning curve
192 void lissajous_trans_v(const Vector &x, Vector &p)
193 {
194  return lissajous_trans(x, p, a, b, delta);
195 }
196
197 // Horizontally spinning curve
198 void lissajous_trans_h(const Vector &x, Vector &p)
199 {
200  return lissajous_trans(x, p, b, a, delta);
201 }
void lissajous_trans_h(const Vector &x, Vector &p)
Definition: lissajous.cpp:198
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:9864
void lissajous_trans_v(const Vector &x, Vector &p)
Definition: lissajous.cpp:192
int main(int argc, char *argv[])
Definition: ex1.cpp:62
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2731
double u_function(const Vector &x)
Definition: lissajous.cpp:142
double b
Definition: lissajous.cpp:42
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
double delta
Definition: lissajous.cpp:43
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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)
Definition: optparser.hpp:76
void lissajous_trans(const Vector &x, Vector &p, double a, double b, double delta)
Definition: lissajous.cpp:148
double a
Definition: lissajous.cpp:41
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1685
int open(const char hostname[], int port)
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
bool Good() const
Definition: optparser.hpp:125