MFEM  v4.1.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-2020, 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  Mesh *mesh;
76  // The mesh could use quads (default) or triangles
77  Element::Type el_type = Element::QUADRILATERAL;
78  // Element::Type el_type = Element::TRIANGLE;
79  mesh = new Mesh(nx, ny, el_type, 1, 2*M_PI, 2.0);
80 
81  mesh->SetCurvature(order, true, 3, Ordering::byVDIM);
82 
83  if (close_strip)
84  {
85  Array<int> v2v(mesh->GetNV());
86  for (int i = 0; i < v2v.Size(); i++)
87  {
88  v2v[i] = i;
89  }
90  // identify vertices on vertical lines (with a flip)
91  for (int j = 0; j <= ny; j++)
92  {
93  int v_old = nx + j * (nx + 1);
94  int v_new = ((close_strip == 1) ? j : (ny - j)) * (nx + 1);
95  v2v[v_old] = v_new;
96  }
97  // renumber elements
98  for (int i = 0; i < mesh->GetNE(); i++)
99  {
100  Element *el = mesh->GetElement(i);
101  int *v = el->GetVertices();
102  int nv = el->GetNVertices();
103  for (int j = 0; j < nv; j++)
104  {
105  v[j] = v2v[v[j]];
106  }
107  }
108  // renumber boundary elements
109  for (int i = 0; i < mesh->GetNBE(); i++)
110  {
111  Element *el = mesh->GetBdrElement(i);
112  int *v = el->GetVertices();
113  int nv = el->GetNVertices();
114  for (int j = 0; j < nv; j++)
115  {
116  v[j] = v2v[v[j]];
117  }
118  }
119  mesh->RemoveUnusedVertices();
120  mesh->RemoveInternalBoundaries();
121  }
122 
123  mesh->Transform(mobius_trans);
124 
125  if (!dg_mesh)
126  {
127  mesh->SetCurvature(order, false, 3, Ordering::byVDIM);
128  }
129 
130  GridFunction &nodes = *mesh->GetNodes();
131  for (int i = 0; i < nodes.Size(); i++)
132  {
133  if (std::abs(nodes(i)) < 1e-12)
134  {
135  nodes(i) = 0.0;
136  }
137  }
138 
139  ofstream ofs(new_mesh_file);
140  ofs.precision(8);
141  mesh->Print(ofs);
142  ofs.close();
143 
144  if (visualization)
145  {
146  char vishost[] = "localhost";
147  int visport = 19916;
148  socketstream sol_sock(vishost, visport);
149  sol_sock.precision(8);
150  sol_sock << "mesh\n" << *mesh << flush;
151  }
152 
153  delete mesh;
154  return 0;
155 }
156 
157 void mobius_trans(const Vector &x, Vector &p)
158 {
159  double a = 1.0 + 0.5 * (x[1] - 1.0) * cos( num_twists * x[0] );
160 
161  p.SetSize(3);
162  p[0] = a * cos( x[0] );
163  p[1] = a * sin( x[0] );
164  p[2] = 0.5 * (x[1] - 1.0) * sin( num_twists * x[0] );
165 }
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 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:407
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:696
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:9864
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
int main(int argc, char *argv[])
Definition: ex1.cpp:62
void RemoveInternalBoundaries()
Definition: mesh.cpp:10020
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
const Element * GetElement(int i) const
Definition: mesh.hpp:775
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:9913
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
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:690
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual int GetNVertices() const =0
Vector data type.
Definition: vector.hpp:48
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
double num_twists
Abstract data type element.
Definition: element.hpp:28
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:779
bool Good() const
Definition: optparser.hpp:125