MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mobius-strip.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 // Mobius Strip Miniapp: Generate Mobius strip meshes
14 // ---------------------------------------------------
15 //
16 // This miniapp generates various Mobius strip-like surface meshes. It is a good
17 // way to generate complex surface meshes. Manipulating the mesh topology and
18 // performing mesh transformation are demonstrated. The mobius-strip mesh in the
19 // data/ directory was generated with this miniapp.
20 //
21 // Compile with: make mobius-strip
22 //
23 // Sample runs: mobius-strip
24 // mobius-strip -t 4.5 -nx 16
25 // mobius-strip -c 1 -t 1
26 // mobius-strip -c 1 -t 4 -nx 16
27 // mobius-strip -c 0 -t 0.75
28 
29 #include "mfem.hpp"
30 #include <fstream>
31 #include <iostream>
32 
33 using namespace std;
34 using namespace mfem;
35 
36 double num_twists = 0.5;
37 void mobius_trans(const Vector &x, Vector &p);
38 
39 int main(int argc, char *argv[])
40 {
41  const char *new_mesh_file = "mobius-strip.mesh";
42  int nx = 8;
43  int ny = 2;
44  int order = 3;
45  int close_strip = 2;
46  bool dg_mesh = false;
47  bool visualization = true;
48 
49  OptionsParser args(argc, argv);
50  args.AddOption(&new_mesh_file, "-m", "--mesh-out-file",
51  "Output Mesh file to write.");
52  args.AddOption(&nx, "-nx", "--num-elements-x",
53  "Number of elements in x-direction.");
54  args.AddOption(&ny, "-ny", "--num-elements-y",
55  "Number of elements in y-direction.");
56  args.AddOption(&order, "-o", "--mesh-order",
57  "Order (polynomial degree) of the mesh elements.");
58  args.AddOption(&close_strip, "-c", "--close-strip",
59  "How to close the strip: 0 - open, 1 - closed, 2 - twisted.");
60  args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
61  "Use discontinuous or continuous space for the mesh nodes.");
62  args.AddOption(&num_twists, "-t", "--num-twists",
63  "Number of twists of the strip.");
64  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
65  "--no-visualization",
66  "Enable or disable GLVis visualization.");
67  args.Parse();
68  if (!args.Good())
69  {
70  args.PrintUsage(cout);
71  return 1;
72  }
73  args.PrintOptions(cout);
74 
75  // The mesh could use quads (default) or triangles
76  Element::Type el_type = Element::QUADRILATERAL;
77  // Element::Type el_type = Element::TRIANGLE;
78  Mesh mesh = Mesh::MakeCartesian2D(nx, ny, el_type, 1, 2*M_PI, 2.0);
79 
80  mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
81 
82  if (close_strip)
83  {
84  Array<int> v2v(mesh.GetNV());
85  for (int i = 0; i < v2v.Size(); i++)
86  {
87  v2v[i] = i;
88  }
89  // identify vertices on vertical lines (with a flip)
90  for (int j = 0; j <= ny; j++)
91  {
92  int v_old = nx + j * (nx + 1);
93  int v_new = ((close_strip == 1) ? j : (ny - j)) * (nx + 1);
94  v2v[v_old] = v_new;
95  }
96  // renumber elements
97  for (int i = 0; i < mesh.GetNE(); i++)
98  {
99  Element *el = mesh.GetElement(i);
100  int *v = el->GetVertices();
101  int nv = el->GetNVertices();
102  for (int j = 0; j < nv; j++)
103  {
104  v[j] = v2v[v[j]];
105  }
106  }
107  // renumber boundary elements
108  for (int i = 0; i < mesh.GetNBE(); i++)
109  {
110  Element *el = mesh.GetBdrElement(i);
111  int *v = el->GetVertices();
112  int nv = el->GetNVertices();
113  for (int j = 0; j < nv; j++)
114  {
115  v[j] = v2v[v[j]];
116  }
117  }
118  mesh.RemoveUnusedVertices();
120  }
121 
122  mesh.Transform(mobius_trans);
123 
124  if (!dg_mesh)
125  {
126  mesh.SetCurvature(order, false, 3, Ordering::byVDIM);
127  }
128 
129  GridFunction &nodes = *mesh.GetNodes();
130  for (int i = 0; i < nodes.Size(); i++)
131  {
132  if (std::abs(nodes(i)) < 1e-12)
133  {
134  nodes(i) = 0.0;
135  }
136  }
137 
138  ofstream ofs(new_mesh_file);
139  ofs.precision(8);
140  mesh.Print(ofs);
141  ofs.close();
142 
143  if (visualization)
144  {
145  char vishost[] = "localhost";
146  int visport = 19916;
147  socketstream sol_sock(vishost, visport);
148  sol_sock.precision(8);
149  sol_sock << "mesh\n" << mesh << flush;
150  }
151 
152  return 0;
153 }
154 
155 void mobius_trans(const Vector &x, Vector &p)
156 {
157  double a = 1.0 + 0.5 * (x[1] - 1.0) * cos( num_twists * x[0] );
158 
159  p.SetSize(3);
160  p[0] = a * cos( x[0] );
161  p[1] = a * sin( x[0] );
162  p[2] = 0.5 * (x[1] - 1.0) * sin( num_twists * x[0] );
163 }
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void mobius_trans(const Vector &x, Vector &p)
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:931
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:11551
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
void RemoveInternalBoundaries()
Definition: mesh.cpp:11707
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
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5312
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
const Element * GetElement(int i) const
Definition: mesh.hpp:1033
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
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11600
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 GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:925
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
virtual int GetNVertices() const =0
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
double num_twists
Abstract data type element.
Definition: element.hpp:28
int main()
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1037
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150