MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
lissajous.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
33using namespace std;
34using namespace mfem;
35
36real_t u_function(const Vector &x);
37void lissajous_trans_v(const Vector &x, Vector &p);
38void lissajous_trans_h(const Vector &x, Vector &p);
39
40// Default Lissajous curve parameters
41real_t a = 3.0;
42real_t b = 2.0;
44
45int 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 {
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 {
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
145{
146 return x[2];
147}
148
149// Tubular Lissajous curve with the given parameters (a, b, theta)
151 real_t a_, real_t b_, real_t delta_)
152{
153 p.SetSize(3);
154
155 real_t phi = x[0];
156 real_t theta = x[1];
157 real_t t = phi;
158
159 real_t A = b_; // Scaling of the curve along the x-axis
160 real_t 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 {
170 real_t R = 0.02*(A+B);
171
172 // normal to the cylinder at p(t)
173 real_t normal[3] = { cos(b_*t), sin(b_*t), 0 };
174
175 // tangent to the curve, dp/dt(t)
176 // real_t 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#ifdef MFEM_USE_SINGLE
180 real_t cn = 1e-32;
181#elif defined MFEM_USE_DOUBLE
182 real_t cn = 1e-128;
183#else
184 real_t cn = 0;
185 MFEM_ABORT("Floating point type undefined");
186#endif
187 real_t cross[3] = { A*a_*sin(b_*t)*cos(a_*t+delta_), -A*a_*cos(b_*t)*cos(a_*t+delta_), b_*B };
188 for (int i = 0; i < 3; i++) { cn += cross[i]*cross[i]; }
189 for (int i = 0; i < 3; i++) { cross[i] /= sqrt(cn); }
190
191 // create a tubular surface of radius R around the curve p(t), in the plane
192 // orthogonal to the tangent (with basis given by normal and cross)
193 for (int i = 0; i < 3; i++)
194 {
195 p[i] += R * (cos(theta)*normal[i] + sin(theta)*cross[i]);
196 }
197 }
198}
199
200// Vertically spinning curve
202{
203 return lissajous_trans(x, p, a, b, delta);
204}
205
206// Horizontally spinning curve
208{
209 return lissajous_trans(x, p, b, a, delta);
210}
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Mesh data type.
Definition mesh.hpp:56
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
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:4243
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:6211
void Transform(void(*f)(const Vector &, Vector &))
Definition mesh.cpp:12821
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.
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int main()
real_t b
Definition lissajous.cpp:42
void lissajous_trans_v(const Vector &x, Vector &p)
void lissajous_trans(const Vector &x, Vector &p, real_t a_, real_t b_, real_t delta_)
real_t delta
Definition lissajous.cpp:43
void lissajous_trans_h(const Vector &x, Vector &p)
real_t u_function(const Vector &x)
real_t a
Definition lissajous.cpp:41
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
const char vishost[]
real_t p(const Vector &x, real_t t)
RefCoord t[3]