MFEM  v4.4.0 Finite element discretization library
snake.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 // Snake Miniapp: Model of the Rubik's Snake_TM Puzzle
14 // ----------------------------------------------------
15 //
16 // This miniapp provides a light-hearted example of mesh manipulation and
17 // GLVis integration.
18 //
19 // The Rubik's Snake a.k.a. Twist is a simple tool for experimenting with
20 // geometric shapes in 3D. It consists of 24 triangular prisms attached in
21 // a row so that neighboring wedges can rotate against each other but cannot
22 // be separated. An astonishing variety of different configurations can be
23 // reached. Enjoy!
24 //
25 // Compile with: make snake
26 //
27 // Sample runs: snake
28 // snake -c 6
29
30 #include "mfem.hpp"
31 #include "../common/mesh_extras.hpp"
32 #include <fstream>
33 #include <iostream>
34
35 using namespace std;
36 using namespace mfem;
37 using namespace mfem::common;
38
39 static int joint_ = 0;
40 static int notch_ = 0;
41 static int step_ = 0;
42 static int nstep_ = 6;
43
44 static double cosa_ = cos(0.5 * M_PI / nstep_);
45 static double sina_ = sin(0.5 * M_PI / nstep_);
46
48
49  Each configuration must be 23 integers long corresponding to the 23 joints
50  making up the Snake_TM puzzle. The values can be 0-3 indicating how far to
51  rotate the joint in the clockwise direction when looking along the snake
52  from the starting (lower) end. The values 0, 1, 2, and 3 correspond to
53  angles of 0, 90, 180, and 270 degrees respectively.
54 */
55 static int conf[][23] =
56 {
57  /* 0 - Ball */ {
58  3,1,3,3,1,3,1,1,
59  3,1,3,3,1,3,1,1,
60  3,1,3,3,1,3,1
61  },
62  /* 1 - Triangle */ {
63  1,0,0,0,0,0,0,3,
64  1,0,0,0,0,0,0,3,
65  1,0,0,0,0,0,0
66  },
67  /* 2 - Hexagon */ {
68  0,0,1,0,0,0,0,3,
69  0,0,1,0,0,0,0,3,
70  0,0,1,0,0,0,0
71  },
72  /* 3 - Snow Flake */ {
73  1,1,1,1,3,3,3,3,
74  1,1,1,1,3,3,3,3,
75  1,1,1,1,3,3,3
76  },
77  /* 4 - Spiral */ {
78  2,1,2,1,2,1,2,1,
79  2,1,2,1,2,1,2,1,
80  2,1,2,1,2,1,2
81  },
82  /* 5 - Zig-zag */ {
83  3,3,3,1,1,1,3,3,
84  3,1,1,1,3,3,3,1,
85  1,1,3,3,3,1,1
86  },
87  /* 6 - Cobra */ {
88  2,0,0,2,1,3,0,2,
89  0,2,3,0,1,3,2,1,
90  1,2,3,1,0,3,0
91  },
92  /* 7 - Serenity */ {
93  3,2,3,2,1,2,0,2,
94  0,2,3,2,1,2,1,2,
95  0,2,3,2,1,2,0
96  },
97  /* 8 - Pinwheel */ {
98  3,2,1,0,2,3,3,2,
99  1,0,2,3,3,2,1,0,
100  2,3,3,2,1,0,2
101  },
102  /* 9 - Crane */ {
103  0,0,3,2,0,2,2,0,
104  0,2,0,0,0,2,2,0,
105  0,0,3,0,0,0,2
106  },
107  /* 10 - Snake */ {
108  0,1,0,1,0,1,0,1,
109  0,1,3,1,0,1,0,1,
110  0,1,0,1,0,1,2
111  },
112  /* 11 - Sculpture */ {
113  0,2,0,2,2,0,3,0,
114  2,2,0,1,0,2,2,0,
115  3,0,2,2,0,2,0
116  },
117  /* 12 - Angles */ {
118  0,2,0,2,2,0,3,0,
119  2,2,0,0,0,2,2,0,
120  3,0,2,2,0,2,0
121  },
122 };
123 static int NUM_CONFIGURATIONS = 13;
124
125 void trans(const int * conf, Mesh & mesh);
126
127 bool anim_step(const int * conf, Mesh & mesh);
128
129 int main(int argc, char *argv[])
130 {
131  int cfg = -1;
132  bool anim = true;
133  bool user = false;
134  bool visualization = true;
135
136  Array<int> myConf(0);
137
138  OptionsParser args(argc, argv);
140  "Select one of 13 pre-programmed configurations: 0-12");
142  "User defined configuration consisting of "
143  "23 joint positions defined by the integers 0, 1, 2, or 3.");
145  "--no-animation",
146  "Enable or disable GLVis animation.");
148  "--no-visualization",
149  "Enable or disable GLVis visualization.");
150  args.Parse();
151  if (!args.Good())
152  {
153  args.PrintUsage(cout);
154  return 1;
155  }
156  args.PrintOptions(cout);
157
158  // Test for a user supplied configuration
159  if (myConf.Size() > 0)
160  {
161  user = true;
162
163  if (myConf.Size() != 23)
164  {
165  MFEM_ABORT("Invalid user-defined configuration of length "
166  << myConf.Size());
167  }
168  }
169
170  // Test for a pre-programmed configuration
171  if (!user && cfg >=0 && cfg < NUM_CONFIGURATIONS)
172  {
173  myConf.SetSize(23);
174  myConf.Assign(conf[cfg]);
175  }
176
177  // Validate the configuration if it has been set
178  for (int i=0; i<myConf.Size(); i++)
179  {
180  if (myConf[i] < 0 || myConf[i] > 3)
181  {
182  MFEM_ABORT("Invalid entry \"" << myConf[i]
183  << "\"in configuration at position " << i);
184  }
185  }
186
187  if (!visualization) { anim = false; }
188
189  // Define an empty mesh
190  Mesh mesh(3, 6 * 24, 24);
191
192  // Add vertices for 24 elements
193  double c[9];
194  int v[6];
195  for (int i=0; i<12; i++)
196  {
197  // Add vertices for a pair of elements
198  // First Upward-facing wedge
199  c[0] = i-6; c[1] = 0.0; c[2] = i-6;
200  c[3] = i-6; c[4] = 0.0; c[5] = i-5;
201  c[6] = i-5; c[7] = 0.0; c[8] = i-5;
202
206
207  c[1] = 1.0; c[4] = 1.0; c[7] = 1.0;
211
212  for (int j=0; j<6; j++) { v[j] = 12 * i + j; }
214
215  // Next Downward-facing wedge
216  c[0] = i-6; c[1] = 0.0; c[2] = i-5;
217  c[3] = i-5; c[4] = 0.0; c[5] = i-4;
218  c[6] = i-5; c[7] = 0.0; c[8] = i-5;
219
223
224  c[1] = 1.0; c[4] = 1.0; c[7] = 1.0;
228
229  for (int j=0; j<6; j++) { v[j] = 12 * i + j + 6; }
231  }
232
233  mesh.FinalizeTopology();
234
235  // Paint elements with alternating colors
236  FiniteElementCollection *fec = new L2_FECollection(0, 3, 1);
237  FiniteElementSpace fespace(&mesh, fec);
238  GridFunction color(&fespace);
239
240  for (int i=0; i<24; i++) { color[i] = (i%2)?1.0:-1.0; }
241
242  // Output the initial mesh to a file
243  {
244  ostringstream oss;
245  oss << "snake-init.mesh";
246  ofstream ofs(oss.str().c_str());
247  ofs.precision(8);
248  mesh.Print(ofs);
249  ofs.close();
250  }
251
253  if (myConf.Size() > 0 && !anim) { trans(myConf.GetData(), mesh); }
254
255  // Output the resulting mesh to GLVis
256  if (visualization)
257  {
258  char vishost[] = "localhost";
259  int visport = 19916;
260  socketstream sol_sock(vishost, visport);
261  sol_sock.precision(8);
262  sol_sock << "solution\n" << mesh << color << "keys Am\n"
263  << "palette 22\n" << "valuerange -1.5 1\n"
264  << "autoscale off\n" << flush;
265
266  // Animate the twists of the selected configuration
267  if (myConf.Size() > 0 && anim)
268  {
269  sol_sock << "pause\n" << flush;
270  cout << "GLVis visualization paused."
271  << " Press space (in the GLVis window) to resume it.\n";
272
273  while (anim_step(myConf.GetData(), mesh))
274  {
275  static int turn = 0;
276  sol_sock << "solution\n" << mesh << color;
277  if (turn++ % 2 == 0)
278  {
279  sol_sock << "pause\n";
280  }
281  sol_sock << flush;
282  }
283  }
284  sol_sock << "autoscale on\n" << "valuerange -1.5 1\n" << flush;
285  }
286
287  // Join the elements together to form a connected mesh
288  MergeMeshNodes(&mesh, 1);
289
290  // Output the resulting mesh to a file
291  {
292  ostringstream oss;
293  if (user)
294  {
295  oss << "snake-user.mesh";
296  }
297  else if (cfg >= 0)
298  {
299  oss << "snake-c" << cfg << ".mesh";
300  }
301  else
302  {
303  oss << "snake-joined.mesh";
304  }
305  ofstream ofs(oss.str().c_str());
306  ofs.precision(8);
307  mesh.Print(ofs);
308  ofs.close();
309  }
310
311  // Clean up and exit
312  return 0;
313 }
314
315 void
316 rotate(double * x)
317 {
318  if (notch_ == 0) { return; }
319
320  double cent[3];
321  Vector cVec(cent,3);
322  Vector xVec(x,3);
323
324  if (joint_%2 == 0)
325  {
326  cVec[0] = -5.5 + joint_ / 2; cVec[1] = 0.5; cVec[2] = 0.0;
328  switch (notch_)
329  {
330  case 1:
331  swap(xVec[0], xVec[1]);
332  xVec[0] *= -1.0;
333  break;
334  case 2:
335  xVec[0] *= -1.0;
336  xVec[1] *= -1.0;
337  break;
338  case 3:
339  swap(xVec[0], xVec[1]);
340  xVec[1] *= -1.0;
341  break;
342  }
344  }
345  else
346  {
347  cVec[0] = 0.0; cVec[1] = 0.5; cVec[2] = -4.5 + joint_ / 2;
349  switch (notch_)
350  {
351  case 1:
352  swap(xVec[1], xVec[2]);
353  xVec[1] *= -1.0;
354  break;
355  case 2:
356  xVec[1] *= -1.0;
357  xVec[2] *= -1.0;
358  break;
359  case 3:
360  swap(xVec[1], xVec[2]);
361  xVec[2] *= -1.0;
362  break;
363  }
365  }
366 }
367
368 void
369 trans(const int * new_conf, Mesh & mesh)
370 {
371  for (int i=0; i<23; i++)
372  {
373  joint_ = i;
374  notch_ = new_conf[i];
375
376  if (notch_ != 0)
377  {
378  for (int k=0; k<6*(i+1); k++)
379  {
380  rotate(mesh.GetVertex(k));
381  }
382  }
383  }
384 }
385
386 void
387 rotate_step(double * x)
388 {
389  if (notch_ == 0) { return; }
390
391  double cent[3], y[3];
392  Vector cVec(cent,3);
393  Vector xVec(x,3);
394  Vector yVec(y,3);
395
396  if (joint_%2 == 0)
397  {
398  cVec[0] = -5.5 + joint_ / 2; cVec[1] = 0.5; cVec[2] = 0.0;
400  switch (notch_)
401  {
402  case 1:
403  yVec[0] = cosa_ * xVec[0] - sina_ * xVec[1];
404  yVec[1] = sina_ * xVec[0] + cosa_ * xVec[1];
405  yVec[2] = xVec[2];
406  break;
407  case 2:
408  yVec[0] = cosa_ * xVec[0] - sina_ * xVec[1];
409  yVec[1] = sina_ * xVec[0] + cosa_ * xVec[1];
410  yVec[2] = xVec[2];
411  break;
412  case 3:
413  yVec[0] = cosa_ * xVec[0] + sina_ * xVec[1];
414  yVec[1] = -sina_ * xVec[0] + cosa_ * xVec[1];
415  yVec[2] = xVec[2];
416  break;
417  }
419  }
420  else
421  {
422  cVec[0] = 0.0; cVec[1] = 0.5; cVec[2] = -4.5 + joint_ / 2;
424  switch (notch_)
425  {
426  case 1:
427  yVec[0] = xVec[0];
428  yVec[1] = cosa_ * xVec[1] - sina_ * xVec[2];
429  yVec[2] = sina_ * xVec[1] + cosa_ * xVec[2];
430  break;
431  case 2:
432  yVec[0] = xVec[0];
433  yVec[1] = cosa_ * xVec[1] - sina_ * xVec[2];
434  yVec[2] = sina_ * xVec[1] + cosa_ * xVec[2];
435  break;
436  case 3:
437  yVec[0] = xVec[0];
438  yVec[1] = cosa_ * xVec[1] + sina_ * xVec[2];
439  yVec[2] = -sina_ * xVec[1] + cosa_ * xVec[2];
440  break;
441  }
443  }
444 }
445
446 bool
447 anim_step(const int * new_conf, Mesh & mesh)
448 {
449  if (notch_ == 2 && step_ == 2 * nstep_) { joint_++; step_ = 0; }
450  if (notch_ != 2 && step_ == nstep_) { joint_++; step_ = 0; }
451  if (joint_ == 23) { return false; }
452  notch_ = new_conf[joint_];
453
454  if (notch_ == 0)
455  {
456  step_ = nstep_;
457  return true;
458  }
459  else
460  {
461  for (int k=0; k<6*(joint_+1); k++)
462  {
463  rotate_step(mesh.GetVertex(k));
464  }
465  }
466  step_++;
467  return true;
468 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1005
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
T * GetData()
Returns the data.
Definition: array.hpp:112
void MergeMeshNodes(Mesh *mesh, int logging)
Merges vertices which lie at the same location.
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:300
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1601
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
constexpr int visport
void rotate(double *x)
Definition: snake.cpp:316
void Assign(const T *)
Copy data from a pointer. &#39;Size()&#39; elements are copied.
Definition: array.hpp:900
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2878
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void rotate_step(char axis, int incr, double *x)
Definition: rubik.cpp:1003
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
Definition: mesh.cpp:1697
Vector data type.
Definition: vector.hpp:60
bool anim_step(char axis, int incr, Mesh &mesh)
Definition: rubik.cpp:1076
int main()
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150