MFEM  v4.6.0
Finite element discretization library
lissajous.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 
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);
87 
88  H1_FECollection fec(order, 3);
89  FiniteElementSpace fes(&mesh, &fec);
90  GridFunction u(&fes);
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);
110 
111  H1_FECollection fec(order, 3);
112  FiniteElementSpace fes(&mesh, &fec);
113  GridFunction u(&fes);
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((int)b, '0') << "'\n" << flush;
136  south << "keys '.0" << std::string((int)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
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)
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
195 {
196  return lissajous_trans(x, p, a, b, delta);
197 }
198 
199 // Horizontally spinning curve
201 {
202  return lissajous_trans(x, p, b, a, delta);
203 }
void lissajous_trans_h(const Vector &x, Vector &p)
Definition: lissajous.cpp:200
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:12045
int main(int argc, char *argv[])
Definition: lissajous.cpp:45
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
void lissajous_trans_v(const Vector &x, Vector &p)
Definition: lissajous.cpp:194
double u_function(const Vector &x)
Definition: lissajous.cpp:144
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
double b
Definition: lissajous.cpp:42
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:5635
double delta
Definition: lissajous.cpp:43
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
void lissajous_trans(const Vector &x, Vector &p, double a_, double b_, double delta_)
Definition: lissajous.cpp:150
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
double a
Definition: lissajous.cpp:41
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
RefCoord t[3]
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
double u(const Vector &xvec)
Definition: lor_mms.hpp:22