MFEM  v4.6.0
Finite element discretization library
klein-bottle.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 // Klein Bottle Miniapp: Generate Klein bottle meshes
14 // ---------------------------------------------------
15 //
16 // This miniapp generates three types of Klein bottle surfaces. It is similar to
17 // the mobius-strip miniapp. The klein-bottle and klein-donut meshes in the
18 // data/ directory were generated with this miniapp.
19 //
20 // Compile with: make klein-bottle
21 //
22 // Sample runs: klein-bottle
23 // klein-bottle -o 6 -nx 8 -ny 4
24 // klein-bottle -t 0
25 // klein-bottle -t 0 -o 6 -nx 6 -ny 4
26 // klein-bottle -t 2
27 
28 #include "mfem.hpp"
29 #include <fstream>
30 #include <iostream>
31 
32 using namespace std;
33 using namespace mfem;
34 
35 void figure8_trans(const Vector &x, Vector &p);
36 void bottle_trans(const Vector &x, Vector &p);
37 void bottle2_trans(const Vector &x, Vector &p);
38 
39 int main(int argc, char *argv[])
40 {
41  const char *new_mesh_file = "klein-bottle.mesh";
42  int nx = 16;
43  int ny = 8;
44  int order = 3;
45  int trans_type = 1;
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(&trans_type, "-t", "--transformation-type",
59  "Set the transformation type: 0 - \"figure-8\","
60  " 1 - \"bottle\", 2 - \"bottle2\".");
61  args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
62  "Use discontinuous or continuous space for the mesh nodes.");
63  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
64  "--no-visualization",
65  "Enable or disable GLVis visualization.");
66  args.Parse();
67  if (!args.Good())
68  {
69  args.PrintUsage(cout);
70  return 1;
71  }
72  args.PrintOptions(cout);
73 
74  // The mesh could use quads (default) or triangles
75  Element::Type el_type = Element::QUADRILATERAL;
76  // Element::Type el_type = Element::TRIANGLE;
77  Mesh mesh = Mesh::MakeCartesian2D(nx, ny, el_type, 1, 2*M_PI, 2*M_PI);
78 
79  mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
80 
81  {
82  Array<int> v2v(mesh.GetNV());
83  for (int i = 0; i < v2v.Size(); i++)
84  {
85  v2v[i] = i;
86  }
87  // identify vertices on horizontal lines (without a flip)
88  for (int i = 0; i <= nx; i++)
89  {
90  int v_old = i + ny * (nx + 1);
91  int v_new = i;
92  v2v[v_old] = v_new;
93  }
94  // identify vertices on vertical lines (with a flip)
95  for (int j = 0; j <= ny; j++)
96  {
97  int v_old = nx + j * (nx + 1);
98  int v_new = (ny - j) * (nx + 1);
99  v2v[v_old] = v2v[v_new];
100  }
101  // renumber elements
102  for (int i = 0; i < mesh.GetNE(); i++)
103  {
104  Element *el = mesh.GetElement(i);
105  int *v = el->GetVertices();
106  int nv = el->GetNVertices();
107  for (int j = 0; j < nv; j++)
108  {
109  v[j] = v2v[v[j]];
110  }
111  }
112  // renumber boundary elements
113  for (int i = 0; i < mesh.GetNBE(); i++)
114  {
115  Element *el = mesh.GetBdrElement(i);
116  int *v = el->GetVertices();
117  int nv = el->GetNVertices();
118  for (int j = 0; j < nv; j++)
119  {
120  v[j] = v2v[v[j]];
121  }
122  }
123  mesh.RemoveUnusedVertices();
125  }
126 
127  switch (trans_type)
128  {
129  case 0: mesh.Transform(figure8_trans); break;
130  case 1: mesh.Transform(bottle_trans); break;
131  case 2: mesh.Transform(bottle2_trans); break;
132  default: mesh.Transform(bottle_trans); break;
133  }
134 
135  if (!dg_mesh)
136  {
137  mesh.SetCurvature(order, false, 3, Ordering::byVDIM);
138  }
139 
140  GridFunction &nodes = *mesh.GetNodes();
141  for (int i = 0; i < nodes.Size(); i++)
142  {
143  if (std::abs(nodes(i)) < 1e-12)
144  {
145  nodes(i) = 0.0;
146  }
147  }
148 
149  ofstream ofs(new_mesh_file);
150  ofs.precision(8);
151  mesh.Print(ofs);
152  ofs.close();
153 
154  if (visualization)
155  {
156  char vishost[] = "localhost";
157  int visport = 19916;
158  socketstream sol_sock(vishost, visport);
159  sol_sock.precision(8);
160  sol_sock << "mesh\n" << mesh << flush;
161  }
162 
163  return 0;
164 }
165 
166 void figure8_trans(const Vector &x, Vector &p)
167 {
168  const double r = 2.5;
169  double a = r + cos(x(0)/2) * sin(x(1)) - sin(x(0)/2) * sin(2*x(1));
170 
171  p.SetSize(3);
172  p(0) = a * cos(x(0));
173  p(1) = a * sin(x(0));
174  p(2) = sin(x(0)/2) * sin(x(1)) + cos(x(0)/2) * sin(2*x(1));
175 }
176 
177 void bottle_trans(const Vector &x, Vector &p)
178 {
179  double u = x(0);
180  double v = x(1) + M_PI_2;
181  double a = 6.*cos(u)*(1.+sin(u));
182  double b = 16.*sin(u);
183  double r = 4.*(1.-cos(u)/2.);
184 
185  if (u <= M_PI)
186  {
187  p(0) = a+r*cos(u)*cos(v);
188  p(1) = b+r*sin(u)*cos(v);
189  }
190  else
191  {
192  p(0) = a+r*cos(v+M_PI);
193  p(1) = b;
194  }
195  p(2) = r*sin(v);
196 }
197 
198 void bottle2_trans(const Vector &x, Vector &p)
199 {
200  double u = x(1)-M_PI_2, v = 2*x(0);
201  const double pi = M_PI;
202 
203  p(0) = (v<pi ? (2.5-1.5*cos(v))*cos(u) :
204  (v<2*pi ? (2.5-1.5*cos(v))*cos(u) :
205  (v<3*pi ? -2+(2+cos(u))*cos(v) : -2+2*cos(v)-cos(u))));
206  p(1) = (v<pi ? (2.5-1.5*cos(v))*sin(u) :
207  (v<2*pi ? (2.5-1.5*cos(v))*sin(u) :
208  (v<3*pi ? sin(u) : sin(u))));
209  p(2) = (v<pi ? -2.5*sin(v) :
210  (v<2*pi ? 3*v-3*pi :
211  (v<3*pi ? (2+cos(u))*sin(v)+3*pi : -3*v+12*pi)));
212 }
void bottle2_trans(const Vector &x, Vector &p)
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
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
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:12045
const Element * GetElement(int i) const
Return pointer to the i&#39;th element object.
Definition: mesh.hpp:1143
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
void RemoveInternalBoundaries()
Definition: mesh.cpp:12203
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
void figure8_trans(const Vector &x, Vector &p)
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1083
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
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:12096
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
int main(int argc, char *argv[])
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double a
Definition: lissajous.cpp:41
virtual int GetNVertices() const =0
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
void bottle_trans(const Vector &x, Vector &p)
Abstract data type element.
Definition: element.hpp:28